如何根据序列id组合FASTA序列？-Java 学习之路

我有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中这样做？

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 '$0 ~ /^>/ {i=substr($0,2); next;} i != -1 {printf "%s", $0 >> "file"i; i=-1;}' gene*

回复于 2024-05-03T02:27:31+08:00

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

1 回答

相关问题