前面我们讲了R批量下载B细胞和T细胞受体VDJ序列文件，那么如何将这些fasta序列读到R里面，方便后面处理呢？今天小编就给大家演示一下如何利用R将fasta序列转成data.frame。我们就用上次下载到的BCR的VDJ序列为例，7个fasta文件存放在BCR_seq文件夹中。
#安装Biostrings包
BiocManager::install("Biostrings")
library("Biostrings")
library(plyr)
filenames=gsub("\\.fasta","",list.files("BCR_seq"))
filepath=list.files("BCR_seq",full.names = T)
#循环读入7个fasta文件额内容
data <- llply(filepath, function(x){
fastaFile <- readDNAStringSet(x)
#获取序列名字，只取前两列
seq_name = do.call(rbind,strsplit(names(fastaFile),"\\|"))
id=seq_name[,1:2]
#获取序列信息，删掉.
sequence = gsub("\\.","",paste(fastaFile))
#生成数据框
df <- data.frame(id, sequence,stringsAsFactors = F)
names(df)=c("ID","name","seq")
df
})
names(data)=filenames
读完之后，data是一个长度为7的list。
其中每一个元素都是一个data.frame。
前面我们讲了四种获取fasta序列长度的方法，其实读到R里面之后，也能获取每条fasta序列的长度。
all_len=lapply(data,function(x){
tmp=apply(x,1,function(x){
c(x[2],x[3],nchar(x[3]))
})
tmp=t(tmp)
colnames(tmp)=c("name","seq","seq_len")
row.names(tmp)=tmp[,1]
tmp
})
最终得到的all_len也是一个长度为7的list
其中每一个元素也是一个data.frame
参考文献
