如何使用序列 ID 提取 FASTA 序列（shell 脚本）答案

【问题标题】：How to extract FASTA sequence using sequence ID (shell script)如何使用序列 ID 提取 FASTA 序列（shell 脚本）
【发布时间】：2016-04-12 18:31:05
【问题描述】：

我有以下序列，它们是 fasta 格式，带有序列头及其核苷酸。如何比较两个文件（Kcompare.pep 和 clade1i.txt）并提取具有相同序列头的序列？谁能帮帮我？

Kcompare.pep

>ztr:MYCGRDRAFT_45998
MAAPLHAEGPIRTPYTGVELLNTPYLNKGTAFPADERRVLGLTALLPTSVHTLDQQLQRA
WHQYQSRDNDLARNTFLTSLKEQNEVLYYRLVLDHLSEVFSIIYTPTEGEAIQRYSSLFR
>kal:KALB_5042
MTAEVAVVSDGSAIPGASPPATLPLLQDYAELVREHAGLSAVPLAVDSARLAAELCALPK
RFRAVFLTHTDPERAFQVQRAVAKAGGPLVITDQDTTAISLTASTLTTLARRGRSPSDSR

clade1i.txt

cpo:COPRO5265_0583
ble:BleG1_3845
kal:KALB_5042

预期输出

>kal:KALB_5042
MTAEVAVVSDGSAIPGASPPATLPLLQDYAELVREHAGLSAVPLAVDSARLAAELCALPK
RFRAVFLTHTDPERAFQVQRAVAKAGGPLVITDQDTTAISLTASTLTTLARRGRSPSDSR

我尝试运行它，但没有出现错误或结果。

for i  in K*
   do
     echo $i

      awk -F ' ' '{print $1}' $i/$i.pep > Kcompare.pep
      mv Kcompare.pep $i

      awk -F '_' '{print $2":"$3"_"$4}' $i/firstClade.txt > $i/clade1i.txt
       awk 'NR==1{printf $0"\t";next}{printf /^>/ ? "\n"$0"\t" : $0}' $i/Kcompare.pep | awk -F"\t"  'BEGIN{while((getline k <"$i/clade1i.txt")>0)i[k]=1}{gsub("^>","",$0);if(i[$1]){print ">"$1"\n"$2}}' > $i/firsti.pep
    done

【问题讨论】：

biostars.org 上的许多解决方案：biostars.org/p/84983/#84996
您也可以尝试使用SAMtools。

标签： shell fasta

【解决方案1】：

使用这个：

while read l; do
  sed -n '/^>'"$l"'/,/^>|$/p' Kcompare.pep
done <clade1i.txt

while 循环在clade1i.txt 文件中逐行循环。
sed -n 禁止自动打印。
/regex/,/regex/ 匹配从第一个正则表达式到第二个正则表达式的所有内容。
p 打印匹配的行。

【讨论】：

循环遍历文件并重新读取每一行的输入文件是一种常见的反模式。您可以轻松地将输入转换为单个 sed 脚本，然后在输入文件上运行该脚本。 sed 's%.*%/^>&$/,/^>|$/p' clade1i.txt | sed -n -f - Kcompare.pep

【解决方案2】：

使用 awk：

awk 'NR==FNR{a[">"$0];next}/^>/{f=0;}($0 in a)||f{print;f=1}' clade1i.txt Kcompare.pep

读取 clade1i.txt 文件并作为键存储在数组中。阅读 Kcompare.pep。对于以 '>' 开头的每一行，设置一个标志，并继续打印这些行，直到遇到以 '>' 开头的下一行。

【讨论】：