Shell-Script

for循環在單個命令中將多個文本從父文件grep到多個文件

  • April 5, 2019

我有 29 個 fasta 文件(.fa 作為副檔名)根據它們的基因命名和儲存序列。

(例如:核醣體蛋白 L1、核醣體蛋白 L6P/L9E、…)

在這29個fasta文件中,共有722種存在。每個序列都在第一行標記了它們的基因和物種名稱,第二行用它的序列填充。

1個物種將有超過1個基因序列。

我想將根據基因排序的 29 個 fasta 文件中的 722 個物種轉移到單獨的 722 個文件中(在物種而不是基因下對它們進行排序)。

父文件中的物種名稱用方括號括起來[ ]

如何使用 for 循環提取 722 個文件並根據其序列名稱命名文件?

範例來自Ribosomal Protein L1.fa

>gi|103486926|ref|YP_616487.1| 50S ribosomal protein L1 [Sphingopyxis alaskensis RB2256]
MAKLTKKQKALEGKVDAQKLHGVDEAIKLVRELATAKFDETLEIAMNLGVDPRHADQMVRGVVTLPAGTGKDVKVAVFAR

範例來自Ribosomal Protein L6PL9E.fa

>gi|410479108|ref|YP_006766745.1| ribosomal protein L6P/L9E [Leptospirillum ferriphilum ML-04]
MGFTHTVEFTLPSLIKASIEKQTIITLSSPDKELLGQFAADVRSIRPPEPYKGKGIKYSGEKILRKEGKTGKK

對於第一個例子,

種名:Shingopyxis alaskensis RB2256

基因序列:MAKLTKKQKALEGKVDAQKLHGVDEAIKLVRELATAKFDETLEIAMNLGVDPRHADQMVRGVVTLPAGTGKDVKVAVFA

我想將文件命名為Sphingopyxis alaskensis RB2256.fa並將具有此物種名稱的所有序列插入此文件中。

我正在使用 bash shell 來執行此操作。我可以grep用來做事:

grep -A+1 "Sphingopyxis alaskensis RB2256" *.fa >> Sphingopyxis alaskensis RB2256.fa

但是我需要做 722 次才能讓我的序列根據物種進行分類。

是否可以使用 for 循環中的 grep 來簡化工作?或者有其他方法可以做到這一點?

Fasta 格式不要求所有序列都在一行上。事實上,這並不常見,因為大多數生物序列都很長。因此grep,在 ID 有不止一行序列的任何情況下,您都會失敗。此外,您的grep命令將創建一個名為的文件Sphingopyxis,而不是一個名為Sphingopyxis alaskensis RB2256.fa.

在任何情況下,您都可以執行以下操作將每個序列放入物種之後的文件名中:

awk -F'[][]' '/>/{n=$2}; {print >> n".fa"}' *.fa 

但是,我強烈建議您不要在文件名中使用空格,因為這只會讓您的生活更加艱難。更安全的方法是:

awk -F'[][]' '/>/{n=$2; gsub(/ /,"_",n)}; {print >> n".fa"}' *.fa 

gsub替換物種名稱中的所有空格_,生成以下文件:

Leptospirillum_ferriphilum_ML-04.fa  Sphingopyxis_alaskensis_RB2256.fa

請注意,上述兩種方法都可以處理多行序列。

引用自:https://unix.stackexchange.com/questions/510674