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

生信百科 · 公众号 · 医学 · 2017-07-26 09:27

正文

上期链接：

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

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

各位看官实在抱歉，最近忙于私事，未能及时更新这个标记开发相关主题。今天总算可以继续为大家继续更新。

上期我们讨论了如何进行前两步开发，这期我们将展示如何进行第三步操作。

三、对候选位点进行高通量引物设计。

进行这步分析需要两步走，首先将每个候选位点上下250bp的序列提取出来存为fasta文件，然后对每个fasta序列进行引物设计。

用R进行上下250bp的序列提取工作，代码如下：（运行前请自行安装seqinr这个包）

fa=read.fasta("genome.fasta")

files = dir("availableSite/")

for(j in 1:length(files))
{
    name = strsplit(files[j],"_")[[1]][1]
    final = read.table(sprintf("availableSite/%s",files[j]),head=T)
    posForMarker = final[which(final$pos!="N"),]
    posForMarker$pos = as.numeric(as.character(posForMarker$pos))
    posForMarker 0){
        dir.create(sprintf("fastaFile/%s",name))
        for(i in 1:nrow(posForMarker)){
            chr = posForMarker$V8[i]
            seq = toupper(fa[[chr]][(posForMarker$pos[i]-250):(posForMarker$pos[i]+250)])
            write.fasta(sequence=seq,names=paste(posForMarker$V8[i],posForMarker$pos[i],sep="_"),file.out=sprintf("fastaFile/%s/%s_%s_%s.fa",name,name,posForMarker$V8[i],posForMarker$pos[i]))
        }    
    }
    print(j)
}

用primer3进行引物设计。（由于需要设计的引物实在是太多，这一步一定要用多CPU服务器，CPU越多越好。。。。）

这步运行我们没有直接写多线程程序，而是对应每个fasta序列做一个bat运行文件，所以有多少个位点就有多少bat运行程序，然后写个程序将其分配给每个CPU进行运行。所以CPU越多，意味着你最后所用的时间越少。

bat文件制作过程如下：

#make bat file for runing primer3
#before running you need make two folders named 'bat' and 'keyfiles'
#input folder include all fasta file which is generated last step

import sys, glob, os

folder = sys.argv[1]

Lst = glob.glob(folder+"*")

for enzyme in Lst:
    name = enzyme.split("/")[-1]
    Lst2 = glob.glob(enzyme+"/*")
    for fa in Lst2:
        Fopen = open(fa,'rU')
        Fopen.readline()
        seq=Fopen.read().replace("\n","")
        seqName = fa.split("/")[-1].split(".")[0]
        Fopen.close()
        keyfile="""SEQUENCE_ID=%s
SEQUENCE_TEMPLATE=%s
SEQUENCE_TARGET=240,20
PRIMER_TASK=generic
PRIMER_PICK_LEFT_PRIMER=1
PRIMER_PICK_INTERNAL_OLIGO=0
PRIMER_PICK_RIGHT_PRIMER=1
PRIMER_OPT_SIZE=20
PRIMER_MIN_SIZE=18
PRIMER_MAX_SIZE=25
PRIMER_DNTP_CONC=0.8
PRIMER_DNA_CONC=500
PRIMER_SALT_CORRECTIONS=1
PRIMER_SALT_MONOVALENT=50
PRIMER_SALT_DIVALENT=1.5
PRIMER_MAX_END_GC=2
PRIMER_MIN_GC=35
PRIMER_MAX_GC=65
PRIMER_MIN_TM=58
PRIMER_MAX_TM=64
PRIMER_OPT_TM=60
PRIMER_PAIR_MAX_DIFF_TM=5
PRIMER_MAX_NS_ACCEPTED=1
PRIMER_PRODUCT_SIZE_RANGE=200-300
P3_FILE_FLAG=1
SEQUENCE_INTERNAL_EXCLUDED_REGION=240,20
PRIMER_EXPLAIN_FLAG=1
=\n"""%(seqName,seq)
        out=open("keyfiles/"+seqName,'w')
        bat = open("bat/"+seqName+".bat",'w')
        bat.write("~/programs/Primer3/primer3-2.3.7/src/primer3_core < %s > %s"%("keyfiles/"+seqName,seqName+".primers"))
        cmd = "chmod 775 %s"%("bat/"+seqName+".bat")
        os.system(cmd)
        out.write(keyfile)
        out.close()
        bat.close()

曾经设计过引物的童鞋，上面参数应该都可以看懂，这里不做过多注释，不明白的自行查阅一下引物设计规则吧。代码后部分加了一个chmod 775 xxx,这是为了将普通bat文件转为可运行文件。请注意修改你安装primer3的软件路径。

到这里你可以随便挑一个bat文件试运行一下，如果可以运行，说明是没问题的。

下一步是将所有bat文件分配到各个CPU进行单独运算。每个实验室应该都有自己分配CPU任务的方法，我贴一个供参考：

library(seqinr)
fa=read.fasta("genome.fasta")

files = dir("availableSite/")

for(j in 1:length(files))
{
    name = strsplit(files[j],"_")[[1]][1]
    final = read.table(sprintf("availableSite/%s",files[j]),head=T)
    posForMarker = final[which(final$pos!="N"),]
    posForMarker$pos = as.numeric(as.character(posForMarker$pos))
    posForMarker 0){
        dir.create(sprintf("fastaFile/%s",name))
        for(i in 1:nrow(posForMarker)){
            chr = posForMarker$V8[i]
            seq = toupper(fa[[chr]][(posForMarker$pos[i]-250):(posForMarker$pos[i]+250)])
            write.fasta(sequence=seq,names=paste(posForMarker$V8[i],posForMarker$pos[i],sep="_"),file.out=sprintf("fastaFile/%s/%s_%s_%s.fa",name,name,posForMarker$V8[i],posForMarker$pos[i]))
        }    
    }
    print(j)
}

运行起来大家可能会发现有些序列花几十分钟甚至几个小时都设计不出引物。这可能是primer3的bug，其实对于一个500bp的序列，我觉得1分钟以上还不能设计出引物的话可能就有问题，这样的序列可以直接扔掉了。。。

根据序列多少及服务器情况这个工作可能需要几天至一两周完成，所以这期就先到这。下期我们继续展示如何对引物进行本地blast及数据展示。

欢迎分享!