拼接序列的编程

      今日偶然听云龙兄说他写了一个拼接序列的脚本。我一心血来潮也想将这个在R语言上实现。随即我要来了云龙兄的意（测）大（序）利（数）炮（据）。

        数据是由云龙兄的师姐预处理好了的，删掉了每个序列差异部分（差异部分很少），听云龙兄说这个些序列需要分类拼接为尽可能长的序列。由于数据太多，能够批量实现是最好。

        知道了前提条件我也有了思路，首先用MEGA将序列分类对齐，然后输出FASTA格式。为了教学需要，我进以前6个为例

         这个问题可以抽象成为在一个矩阵上做投影的问题。即将一类序列构建成一个矩阵，如果某条或者某几条序列在某个位点有氨基酸或者核酸X，其余序列在该位点为gap“-”那边这个位点就为就为X。如果所有序列在这个位点均为gap”-“那么这个位点就是gap”-“。对于序列的读取依然使用R语言的第三方包Biostrings。

library(Biostrings)

seq<-readAAStringSet(“云龙兄.fasta”)

seq1<-seq[1:6,]

        当提取了这6个序列的fasta格式后，需要将序列部分单独提取出来组成一个矩阵，要完成这个目标，由需要将提取出的序列按字符拆分，这样才能保证每个序列每个位点能输入矩阵中

AAseqarry<-NULL

for (i in 1:length(seq1)) {

  AAseq<-toString(seq1[i])

  AAseq2<-strsplit(AAseq,””)

  AAseqarry<-rbind(AAseqarry,AAseq2[[1]])

}

        我们得到了AAseqarry，这里面储存了我们的序列矩阵

        接下来需要对列进行分析，如果该列中有元素属于氨基酸，那么保留该元素至长序列的这个位点，对于每一列是一个向量，我们可以构建一个氨基酸序列的向量，如果两个向量有交集代表这个位点有氨基酸，如果没有交集则代表这个位点是gap。

LAAseq<-NULL

AA<-c(“A”,”R”,”D”,”C”,”Q”,”E”,”H”,”I”,”G”,”N”,”L”,”K”,”M”,”F”,”P”,”S”,”T”,”W”,”Y”,”V”,”*”)

for (i in 1:ncol(AAseqarry)) {

  xbp<-AAseqarry[,i]

    if(length(intersect(AA,xbp))>=1){

      LAAseq<-c(LAAseq,intersect(AA,xbp))

    }else{LAAseq<-c(LAAseq,”-“)}

}

        LAAseq储存了拼接好的序列信息，但是这里面每个位点的信息是分开的，需要用将每个位点信息无分隔地粘贴到一起，再转化为fasta格式，最后输出

LAAfasta<-NULL

for (i in 1:length(LAAseq)) {

  LAAfasta<-paste(LAAfasta,LAAseq[i],sep = “”)

}

LAAfasta<-AAStringSet(LAAfasta)

names(LAAfasta)<-names(seq1[1])

AAout<-tempfile()

writeXStringSet(LAAfasta,AAout)

       上图就是最后输出的拼接了6条序列的长序列。如果在做实验过程中拿到了较多需要拼接序列的数据就可以使用这个代码批量完成。只需要修改你的输入文件名称即可。下面附上原始代码

        云龙兄他的代码也很好地完成了这个分析，日后云龙兄也很加入进来共同分享生物信息学知识。

         下面几期内容将逐渐开展Linux及其该系统上的一些软件使用介绍。会融入在实际生物学问题中进行讲解。