[生信媛]练习题：如何根据FASTA的ID拆分数据
[生信媛]公众号
给定一个FASTA文件，他以多行形式存放数据，如下所示
请根据每个序列的ID中的染色体编号拆分数据，例如上面两个数据，需要新建两个文件，分别是chr7和chr8，将这两条序列分别输出到这两个文件中。
测试数据http://share.weiyun.com/5ecURT0 密码：pzsiab
=========================
R语言版的解决方案：
方案1：
library(readr,stringr) ## str_extract函数需要用到
gene_read <- function(x){
g <- readLines(x) #文件按行读入便于按行操作
num_l <- length(g) # 行数
star <- substr(g,1,1) # 每1行的第1个字符向量
gene <- "" ###用来存放单独的fasta格式单元
for(i in 1:num_l){
line <- paste0(g[i],"\r\n") #在每一行数据后加换行符并赋值给中间变量line
if (star[i]==">") { #通过“>”判断是否fasta格式首行
chrN <- stringr::str_extract(gene,"chr([0-9]{1,2})") #提取染色体编号
readr::write_file(gene,paste0(chrN,".txt"),append = T)# 在遇到下一个">"表示上一个fasta格式单元已经加载完毕。将其写入文件，使用了`append=T`参数，采用了追加模式。默认`append=FALSE` 使用的是替换模式，只会保留最后的书写内容。
gene <- line ##数据写到文件后，gene中间变量才可以存储新的fasta格式单元。
} else{gene <- paste0(gene,line)} # 如果没有“>”则不是fasta格式首行，将这一行内容添加在上一行的后面
}
file.remove("NA.txt") ##删除第1个循环创建的空文件
}
方案2：
# 这一版的方案与上一个方案基本相同，只是为了避免使用'tidyverse' 包，尽量使用`base`包看是否可以提高运行速度。
gene_read1 <- function(x){
g <- readLines(x)
num_l <- length(g)
star <- substr(g,1,1)
gene <- ""
for(i in 1:num_l){
line <- paste0(g[i],"\n")
if (star[i]==">") {
position <- regexec("chr([0-9]{1,2})",gene) ###获得染色体编号的匹配项的位置
chrN <- substr(gene,position[[1]][1],position[[1]][2])##通过位置信息，提取染色体编号
write.table(gene,paste0(chrN,".txt"),row.names=F, col.names=F,quote=F,eol="", append = T)
gene <- line
} else{gene <- paste0(gene,line)}
}
file.remove("NA.txt")
}
system.time(gene_read1("test.faa"))
第三种方案：
library(readr) ###使用wrIte_file函数来写文件
gene_read2 <- function(x){
g <- readLines(x)
num_l <- length(g)
star <- substr(g,1,1)
gene <- ""
for(i in 1:num_l){
line <- paste0(g[i],"\n")
if (star[i]==">") {
position <- regexec("chr([0-9]{1,2})",gene) ###获得染色体编号的匹配项的位置
chrN <- substr(gene,position[[1]][1],position[[1]][2])##通过位置信息，提取染色体编号
readr::write_file(gene,paste0(chrN,".txt"), append = T)
gene <- line
} else{gene <- paste0(gene,line)}
}
file.remove("NA.txt")
}
# 比较一下两种方式的运行速度。
system.time(gene_read1("test.faa"))
用户 系统 流逝
0.95 1.17 1.26
system.time(gene_read("test.faa"))
用户 系统 流逝
0.47 0.65 0.88
system.time(gene_read2("test.faa"))
用户 系统 流逝
0.45 0.75 0.66
通过比较，第3种方案最好，第二种方案比第一种运行速度慢好多，主要问题出在写文件环节，使用readr::write_file函数来写文件速度就提高很多。
经过查看文件发现上述代码有误，生成的文件中常规染色体chr*，而没有考虑scaffold。
需要将代码修改为如下：
gene_read3 <- function(x){
g <- readLines(x)
num_l <- length(g)
star <- substr(g,1,1)
gene <- ""
for(i in 1:num_l){
line <- paste0(g[i],"\r\n")
if (star[i]==">") {
chrN <- str_split(gene,"[>:]")[[1]][2] ### 从fasta格式的首行进行分列，以'>' 或‘:' 分隔，第1列空，第2列为染色体名称。
readr::write_file(gene,paste0(chrN,".txt"), append = T)
gene <- line
} else{gene <- paste0(gene,line)}
}
file.remove("NA.txt")
}
system.time(gene_read3("test.faa"))
用户 系统 流逝
0.45 0.78 0.68
由于对base的文本处理函数不熟悉。欢迎小伙伴来拍砖。
