一个R语言的实用脚本

        昨天师姐想把本地Blast输出的结果文件回溯到megahit输出的拼接数据中抓取目标序列。先通过一波图文解释一下这个过程。首先，实验室拿到了120个X病毒的测序样本(原始的reads数据)，师姐在服务器上通过运行megahit软件将这120个样本分别组装成了长序列文件。

         之后师姐又将这些长序列文件进行本地Blast判断这些序列能与X病毒的哪个片段匹配。本地Blast输出的m6格式文件如下

        可以看出第一列是自己的长序列id，第二列是匹配到的X病毒参考序列id，第三列是自己的长序列长度，第五列是匹配上的X病毒参考序列的长度。由于此文件没有序列信息。要想拿回自己序列的信息只能从第一幅图片中的fasta序列中拿回序列信息。另外，X病毒的最短片段长度也有800bp，所以师姐想抓取出尽可能是全长的序列。由于有120个样本，每个fasta文件中有至少2万条长序列。如果一条一条人工找将会很麻烦。虽然实验室的海舟老师自己写了满足这个需求的perl脚本。但是本着探索的精神，我想独立将这个需求在R语言上实现。

        首先，由于每个fasta文件对应一个blast文件，而且两个文件名是一一对应的。可以通过批量循环读取两个文件。本地Blast输出的m6格式文件可以由R语言的read.table函数读取，fasta文件需要通过加载R语言第三方包Biostrings的readDNAStringSet函数读取。对于R语言的第三方包有两个下载途径，一个是通过install.package(‘yourpackage’)函数从R的源下载，另外一个可以通过Bioconductor官网下载。Bioconductor下载包方法如下：

如果你的R版本在3.5.0以上你通过

BiocManager::install(‘yourpackage’)

如果你的R版本在3.4或者以下你通过

source(“https://bioconductor.org/biocLite.R”)

options(BioC_mirror=”https://mirrors.ustc.edu.cn/bioc/”)

BiocInstaller::biocLite(‘yourpackage’)

       在读取完文件后，首先根据myqueryid即blast文件中自己的长序列id先判断这条序列匹配上了多少bp，如果大于800(因为X病毒最短片段长度为800)，就将其在对应样本编号的长序列fasta文件中的fasta序列注释信息逐一寻找，找完所有对应的id，用DNAStringSet函数输出fasta序列，将这些序列逐条放入输出变量outseqfile中，最后批量将outseqfile逐一用writeXStringSet函数输出成fasta格式的文件。以下是详细的代码已经代码对应的图片

library(stringr)

library(Biostrings)

setwd(“D:/test”)

dir1.name<-dir(‘D:/test/megahit_res’)

dir2.name<-dir(‘D:/test/blastn’)

a<-NULL

b<-NULL

for (i in 1:length(dir1.name)) {

  setwd(“D:/test/megahit_res”)

  myqueryseq<-readDNAStringSet(dir1.name[i])

  setwd(“D:/test/blastn”)

  myqueryid<-read.table(dir2.name[i],sep = “\t”)

  outseq<-DNAStringSet()

  Outseq<-DNAStringSet()

  for (j in 1:nrow(myqueryid)) {

    if(myqueryid[j,5]>=800){

      for (k in 1:length(myqueryseq)) {

        idname<-strsplit(names(myqueryseq[k]),” “)

        if(myqueryid[j,1]==idname[[1]][1]){

          outseq<-c(outseq,DNAStringSet(myqueryseq[k]))

        }else{b<-b+1}

      }

    }else{a<-a+1}

  }

  outseqfile<-tempfile()

  writeXStringSet(outseq,outseqfile)

  outseqfile

}

       最后根据输出文件的时间排序，将文件名改为对应的样本的名称及fasta后缀即可。

       因为writeXStringSet函数没有添加文件名功能，所以在这一步需要手动改文件名，看来R在处理文本方面弱于perl啊。另外，如果确定文件中的id都是以非冗余的形式存在可以将for循环改为while循环，可以减少程序运行的时空复杂度。

      对于R语言小白或者没有生信背景的人，如果在实验中遇到这样批量处理序列的情况，可以完全使用这个代码，只需要将dir1.name<-dir(‘D:/test/megahit_res’)，dir2.name<-dir(‘D:/test/blastn’)和setwd()设置为你自己的工作路径和文件夹名称即可。如果想学R语言可以从Thr Art of R Programming这本书开始学。