如何根据序列id组合FASTA序列？

Question

我有9个FASTA文件，代表9个基因的DNA测序。

每个FASTA文件包含121个序列，代表121个菌株。每个序列的名称是每个菌株的 ID。

但是在每个文件中，id是没有排序的，比如在gene1.fasta:

>1
AAA
>16
TTT
>2
GGG
...

在gene2.fasta中：

>2
CCC
>34
AAA
>1
GGG
...

我想把这9个基因的FASTA文件改成121个菌株的FASTA文件，在每个文件中，简单地把一个菌株的9个基因组合起来。例如，在 strain1.fasta:

AAAGGG

在strain2.fasta中：

GGGCCC

我如何在 R 中执行此操作？

Answer 1

这里有一个 R 中的解决方案，使用 Biostrings 包读取 fasta 文件。

它有效，但我不得不说这是我很长时间以来编写的最丑陋的代码。我只是想看看我能否以某种方式完成这项工作 - 这 100% 不是最佳解决方案。

library("Biostrings")
library("tidyverse")

convertStringSet = function(seq){
  return(df = data.frame("names" = names(seq), "seq" = paste(seq)))
}

# change the path accordingly
filenames = list.files("/home/x/y/z", pattern="gene*", full.names=TRUE)%>%
  lapply(readDNAStringSet)

fastaDF = filenames %>% lapply(convertStringSet) %>% 
  reduce(full_join, by = "names") %>% 
  unite("seq", -1,  sep="")

writeOutput = function(x){

  header = paste(">",x[1],sep="")
  fileName = paste("strain",x[1],".fasta",sep="")

  writeLines(c(header,x[2]), fileName)
}

apply(fastaDF, 1, writeOutput)

作为替代方案，如果您使用的是 unix 系统，此 awk 行应该给出相同的结果：

awk '[=11=] ~ /^>/ {i=substr([=11=],2); next;} i != -1 {printf "%s", [=11=] >> "file"i; i=-1;}' gene*

如何根据序列id组合FASTA序列？

How to combine FASTA sequence according to sequence id?

r

bioinformatics

fasta