☰

使用 while read 从列表中提取 fasta 条目

tag-icon • tag-icon text-processing sed bioinformatics
使用 while read 从列表中提取 fasta 条目

我有 28 个文件，每个文件都有大约 14,000 个“条目”。单个条目由一个标头（用 >string 表示）、一个换行符和一个字符串序列组成。每个条目都有可变长度的序列/字符串。所有 28 个文件都有相同的条目标头，但每个条目的顺序是可变的。

例如，一个文件 CR1_ref.fasta 看起来像

>FBgn0080937
ATGGATAAAAGGCTCAGCGATAGTCCCGGAGATTGTCGCGTAACCAGATCCAGCATGACGCCCACCCTCCGCTTGGAGCACAGTCCCCGGCGGCAACAACAGCAACAACA
>FBgn0076379
ATGCTGCGCACCCTTTTCGCCGTGCGTGGTCAGTGCCAGCAGCTGCTGAGGAGAACATTCACCCCCCATTGCAGTGGCCAACGA
>FBgn0070974
ATGCAGACGCGTCCGAGCAGTGAACCGCAGCGCGCCAAGGAGCAACTCCTGCGGGAGCTGCCGCCGCAGAAATGCTCCAGCGCCACGCTGGCCAAGAAGGTGCTGTCGCAGAGCCCGCCGGCAGCCCCGCCGCCCACACCGGCCACAATTGTGCCGCTCACTGCGGTGCCCGTCATCCAGCTGACGCCTCCGTCGCACTCCGGCGACACGCCGCAAAAGCCAGCACCTCCGGCGCCGCCGCCGCC

总体目标是创建约 14,000 个新文件。其中每个文件都是与所有 28 个文件中的特定 ID/标头关联的条目。

要从单个文件中提取单个条目，我可以使用以下命令

sed -n '/^>FBgn0080937$/{p;n;p;}' CR1_ref.fasta

要在所有 28 个文件中提取此条目（每个文件都以 ref.fasta 结尾），我可以这样做

for i in *ref.fasta; do sed -n '/^>FBgn0080937$/{p;n;p}' $i; done > FBgn0080937.fasta

我有一个单独的文本文件，每行有 14,000 行，对应于名为gene.txt 的条目的标题。该文件的前几行看起来像

FBgn0080937
FBgn0076379
FBgn0070974
FBgn0081668
FBgn0076576
FBgn0076572
FBgn0079684
FBgn0070907
FBgn0080226
FBgn0072746

我想通读该文件，为每个标头 ID 创建一个新的文本文件。 $F 下面提取特定标头 (FBgn*) 的条目并将其存储在新文件中。我正在使用替换命令根据它们来自的 while ref.fasta 文件来重命名序列。

while read -r line;
do F=$line
for i in *ref.fasta
do sed -n "/^>$F$/{s/FB.*/$i/;p;n;p;}" $i > $line.fasta
done
done < "gene.txt"

目前该脚本创建了 14,000 个文件，但每个文件只有一个序列。

>Z9_ref.fasta
ATGCAGACGCGTCCGAGCAGTGAACCGCAGCGCGCCAAGGAGCAAC

我期望每个 *ref.fasta 文件有 28 个序列一个序列。 sed 命令正在输出最后一个条目。预期输出是

    >CR1_ref.fasta
    ATGCAGACGCGTCCGAGCAGTGAACC
    >FH2_ref.fasta
    AGCAGTGAACCGCAGCGCGCCAAGGAGCAAC
    >MSH10_ref.fasta
    CGCGTCCGAGCAGTGAACCGCAGCGCGCCAAGGAGCAAC
    >Z9_ref.fasta
    ATGCAGACGCGTCCGAGCAGTGAACCGCAGCGCGCCAAGGAGCAAC

答案1

shell 并不真正适合这种类型的解析。您可以在自己的代码中看到，您读取了整个文件一次每个从文件中读取的基因名称gene.txt。

下面的单个awk命令可以更快地完成同样的事情。

awk -F '>' '
    FNR == NR           { genes[$1]; next }
    /^>/ && $2 in genes { if (out != "") close(out);
                          out = $2 ".fa"
                          split(FILENAME, a, "_")
                          $0 = ">" a[1] "_" $2 }
    out != ""           { print >>out }' genes.txt *_ref.fasta

首先读取genes.txt文件并创建一个从中调用的关联数组，genes其中基因名称作为键。

当它到达 Fasta 文件时（代码假设这些文件都被称为类似的东西XXX_ref.fasta），当我们读取 Fasta 标头，并且标头中的基因是我们列表中的键时genes，然后我们从基因名称创建一个输出文件名asgenename.fa并重写标头以包含当前文件名下划线之前的部分。

如果原始标头XXX_ref.fasta是

>genename

那么这将被转化为

>XXX_genename

脚本的最后一部分awk将所有行发送到适当的输出文件。

使用您提供的数据进行测试会生成三个文件：

$ ls *.fa
FBgn0070974.fa FBgn0076379.fa FBgn0080937.fa

$ cat FBgn0076379.fa
>CR1_FBgn0076379
ATGCTGCGCACCCTTTTCGCCGTGCGTGGTCAGTGCCAGCAGCTGCTGAGGAGAACATTCACCCCCCATTGCAGTGGCCAACGA

相关内容