如何连接相同 multiFASTA 文件中的序列,然后将结果打印到新的 FASTA 文件?
Posted
技术标签:
【中文标题】如何连接相同 multiFASTA 文件中的序列,然后将结果打印到新的 FASTA 文件?【英文标题】:How to concatenate sequences in the same multiFASTA files and then print result to a new FASTA file? 【发布时间】:2021-11-26 23:45:38 【问题描述】:我有一个包含 50 多个 FASTA 文件的文件夹,每个文件中包含 2-8 个 FASTA 序列,这是一个示例:
testFOR.id_AH004930.fasta
>AH004930|2:1-128_Miopithecus_talapoin
ATGA
>AH004930|2:237-401_Miopithecus_talapoin
GGGT
>AH004930|2:502-580_Miopithecus_talapoin
CTTTGCT
>AH004930|2:681-747_Miopithecus_talapoin
GGTG
testFOR.id_M95099.fasta
>M95099|1:1-90_Homo_sapien
TCTTTGC
>M95099|1:100-243_Homo_sapien
ATGGTCTTTGAA
它们都是根据 ID 号(在本例中为 AH004930 和 M95099)分组的,我已经使用 HERE 找到的非常方便的 seqkit 代码从原始原始 multiFASTA 文件中提取出来。
我的目标是:
-
使用
cat
将这些序列放在文件中,如下所示:
>AH004930|2:1-128_Miopithecus_talapoin
ATGAGGGTCTTTGCTGGTG
>M95099|1:1-90_Homo_sapien
TCTTTGCATGGTCTTTGAA
(我不关心核苷酸的位置,我关心的是ID和物种名称!)
-
将此结果打印到新的 FASTA 文件中。
理想情况下,我真的很想将所有这 50 个文件压缩成 1 个 FASTA,然后我可以继续过滤/对齐:
GENE_L.fasta
>AH004930|2:1-128_Miopithecus_talapoin
ATGAGGGTCTTTGCTGGTG
>M95099|1:1-90_Homo_sapien
TCTTTGCATGGTCTTTGAA
....
到目前为止,我已经找到了一种方法来实现我想要的,但一次只有一个文件(使用此代码:cat myfile.fasta | sed -e '1!/^>.*/d;' | sed ':a;N;$!ba;s/\n//2g' > output.fasta
,我很遗憾失去了信用的链接)但是很多这些文件名非常相似,所以如果我手动完成,我会不可避免地错过一些/它会太慢。
我试图把它放到一个循环中,它有点在那里!但它的作用是对每个 FASTA 文件进行分类,将其放入一个新文件中,但只保留第一个标题,给我留下大量拼接在一起的序列;
for FILE in *; do cat *.fasta| sed -e '1!/^>.*/d;'| sed ':a;N;$!ba;s/\n//2g' > output.fasta; done
output.fasta
>AH004930|2:1-128_Miopithecus_talapoin
ATGAGGGTCTTTGCTGGTGTCTTTGCATGGTCTTTGAAGGTCTTTGAAATGAGTGGT...
我想知道是否制作一个类似于 HERE 的循环会不会有什么好处,但我真的不确定如何让它在打开新文件后打印每个标题。
我怎样才能对这些序列进行分类,将它们打印到一个新文件中并仍然保留这些标题? 我非常感谢任何关于我在循环中出错的地方以及任何适合 zsh shell 的解决方案的建议!我对任何 python 或 linux 解决方案持开放态度。提前谢谢你
【问题讨论】:
如果您能更详细地解释这些输入是如何形成预期输出的,那么对于我们这些了解 Unix 文本工具但不太了解生物信息学的人来说,您的问题会更清楚。 @tripleee 这是一个很好的观点,我想我现在已经找到了解决方案(由用户@potong 提供),但最好多解释一下。我会为遇到它的人相应地更新我的问题 【参考方案1】:这可能对你有用(GNU sed):
sed -s '1h;/>/d;H;$!d;x;s/\n//2g' file1 file2 file3 ...
设置-s
分别处理每个文件。
复制第一行。
删除任何其他包含>
的行。
将所有其他行附加到第一行。
删除除最后一行之外的这些行。
在文件末尾,交换到副本并删除除第一个以外的所有换行符。
对所有文件重复。
非 GNU seds 的替代方案:
for file in *.fasta; do sed '1h;/>/d;H;$!d;x;s/\n//2g' "$file"; done
注意MacOS sed 可能需要放入脚本并使用-f
选项调用或使用-e
选项(减去;
命令)分成几部分,您的运气可能会有所不同。
或许:
for file in file?; do sed $'1h;/>/d;H;$!d;x;s/\\n/@/;s/\\n//g;s/@/\\n/' "$file"; done
【讨论】:
谢谢你,但我在 Mac OS 上,不能使用 GNU sed - 我会看看我是否能找到 Mac 的选项,然后试一试。我可以检查一下,将*.fasta
放在代码末尾而不是单独放置每个文件是否可行?
@Fitzy yes *
是任何名称的通配符。至于 -s
选项,您可以在 for 循环中调用 sed,请参阅编辑。
我不得不从非 GNU sed 选项中删除 g (它提出了 sed: 1: "1h;/>/d;H;$!d;x;s/\n//2g": more than one number or 'g' in substitute flags
)但是一旦我这样做了,它看起来确实完成了工作!我会确认并赞成这个答案,但我想先与你核实一下 global 命令的问题。非常感谢您的帮助,您是明星!
@Fitzy g
标志表示全局,因此删除它只会替换第二个换行符。我忘记了这是 GNU sed 独有的功能。另一种解决方法是用记录中不存在的东西(可能是@
)替换第一个换行符,然后全局替换所有换行符,最后用换行符替换@
。但是我似乎记得\n
也不被 MacOS sed 接受。因此,您需要逐字插入换行符。 HTH
@Fizy 查看最终编辑【参考方案2】:
不确定我是否完全理解您的问题,但如果您只是想将多个文件中的内容连接到一个文件中,我相信下面的(Python)代码应该可以工作:
import os
input_folder = 'path/to/your/folder/with/fasta/files'
output_file = 'output.fasta'
with open(output_file, 'w') as outfile:
for file_name in os.listdir(input_folder):
if not file_name.endswith('.fasta'): # ignore this
continue
file_path = os.path.join(input_folder, file_name)
with open(file_path, 'r') as inpfile:
outfile.write(inpfile.read())
【讨论】:
非常感谢您提供此解决方案,但它并不完全符合我的要求。它确实保持所有标题完好无损并将它们放入一个新的 fasta 文件中,这很棒,但它不会将每个文件中的序列放在一起。以上是关于如何连接相同 multiFASTA 文件中的序列,然后将结果打印到新的 FASTA 文件?的主要内容,如果未能解决你的问题,请参考以下文章
Java对象序列化,与ObjectOutPutStream一起使用后如何使用相同的文件