高通量CAPS标记开发--实战篇（一）

还是按照惯例回顾一下往期内容，上期链接：

大文章下脚料利用---高通量CAPS标记开发 (一)

上次发文吐槽了一下海量重测序数据利用率低下的问题，其实很多实验室也在努力构建自己的数据库，相信不久的将来有点儿实力的实验室会陆续推出自己的数据库吧。我们这里继续讨论如何进行高通量CAPS标记开发。
上次说到一共分七步可以完成开发。大家不要想着这个工作可以一两天内完成，每一步都是要耗点儿时间的，尤其是基因组很大的情况下。今天我们学习如何进行前两步分析。

一、用所有可用的限制性内切酶对要研究物种的全基因组进行酶切，统计酶切位点。

这里我们用EMBOSS软件的restrict插件进行基因组酶切位点统计。具体如下：

restrict -sequence genome.fa -sitelen 4 -enzymes all -noambiguity -outfile genome.restrict

-sequence: input fasta format genome file

-sitelen: 设置最短酶切特异识别序列长度。如果特异序列短于这个长度，该限制性内切酶将被忽略。

-enzymes: 限制性内切酶符号。这里我们用所有可用的限制性内切酶，所以设置all。

-outfile: 输出文件名

这一步的速度取决于你基因组大小，从半小时到几小时不等吧，输出文件如下：

我们可以看到输出文件信息有酶切位点，限制性内切酶，内切酶识别序列等信息。直接输出文件这种格式用起来不太方便，我们整理一步，将output整理成表格形式，加入染色体信息一列。
格式如下：

python代码如下：

#!/usr/bin/python
#coding=utf-8
#python mergeChrSite.py [genome.restrict] [genome.restrict1]

import sys,re

inFile = open(sys.argv[1],'rU')
outFile = open(sys.argv[2],'w')

tag = 0for line in inFile:
    _line = " ".join(line.rstrip().split())    
    if re.search(r"chr",_line):
        tag = _line.rstrip().split(" ")[2]    
    if re.search(r"^\d",_line):
        outline = "\t".join(_line.split(" ")[0:7])+"\t"+tag+"\n"
        outFile.write(outline)
outFile.close()
inFile.close()

整理好格式后目前的格式后，我们将总文件按酶来分隔开:

python script 如下：