使用 awk 从 fasta 文件中选择序列组的问题

Question

我想对我的 fasta 文件进行子集化，以检索属于给定群体的序列。以下是我的文件示例。

>CLocus_12706_Sample_44_Locus_36326_Allele_0 [JoJo_s113.fq; groupI, 125578, +]
TGCAGCATGCTGGTGAACGCGTCATCATAAGCCTGTTGGCGAGCCAGCAGAAGGCGGCATGGGCAGCACTTAATAGGACGCACGTCCTCTGTGTCA
>CLocus_12706_Sample_46_Locus_34641_Allele_0 [JoJo_s115.fq; groupI, 125578, +]
>CLocus_12706_Sample_69_Locus_37751_Allele_0 [LakeCamp_s033.fq; groupI, 125578, +]
TGCAGCATGCTGGTGAACGCGTCATCATAAGCCTGTTGGCGAGCCAGCAGAAGGCGGCATGGGCAGCACTTAATAGGACGCACGTCCTCTGTGTCA
>CLocus_12706_Sample_70_Locus_33595_Allele_0 [LakeCamp_s034.fq; groupI, 125578, +]
TGCAGCATGCTGGTGAACGCGTCATCATAAGCCTGTTGGCGAGCCAGCAGAAGGCGGCATGGGCAGCACTTAATAGGACGCACGTCCTCTGTGTCA
>CLocus_72879_Sample_136_Locus_80036_Allele_0 [NaknekRiver_s148.fq; groupV, 11333693, -]
TGCAGAACGAGATGAGGACAAACACACTCACCACTCTGTGGACATGTAGACGGCTGGCCTGTCCTACCAAGGACAAATACTCCCACAACAGTCCAA

人口是 id 的一部分，例如“LakeCamp”或“JoJo”或“NaknekRiver”。

我试图按照这篇文章来弄清楚如何提取序列。 https://unix.stackexchange.com/questions/253499/extracting-subset-from-fasta-file

为此，我执行了以下操作，“JoJo”是这里选择的人口，我的输入文件是“fasta8c18subset.fa”。

awk -vrs=">" 'BEGIN{t["JoJo"]=1}{if($1 in t){printf ">%s",$0}}' fasta8c18subset.fa

我运行此程序时没有收到错误，但也没有输出。

作为输出，我想获得与该群体相关的整个标题和序列。因此，例如，如果我尝试提取“LakeCamp”样本，我希望输出文件包含以下内容

>CLocus_12706_Sample_69_Locus_37751_Allele_0 [LakeCamp_s033.fq; groupI, 125578, +]
TGCAGCATGCTGGTGAACGCGTCATCATAAGCCTGTTGGCGAGCCAGCAGAAGGCGGCATGGGCAGCACTTAATAGGACGCACGTCCTCTGTGTCA
>CLocus_12706_Sample_70_Locus_33595_Allele_0 [LakeCamp_s034.fq; groupI, 125578, +]
TGCAGCATGCTGGTGAACGCGTCATCATAAGCCTGTTGGCGAGCCAGCAGAAGGCGGCATGGGCAGCACTTAATAGGACGCACGTCCTCTGTGTCA

想法？

Answer 1

我建议您使用可用的 fasta 格式解析器。

例如，在 python3 中，您可以使用来自pyGATB 的非常高效的解析器。

你可以这样使用它：

#!/usr/bin/env python3

import sys
from gatb import Bank

fasta_file = sys.argv[1]
pop_name = sys.argv[2]

def get_pop(header):
    """Extracts the population name from the fasta header."""
    return header.decode("utf-8").split(" ")[1].split("_")[0][1:]

for seq in Bank(fasta_file):
    if get_pop(seq.comment) == pop_name:
        print(">%s\n%s" % (
            seq.comment.decode("utf-8"),
            seq.sequence.decode("utf-8")))

sys.exit(0)

使用您的示例文件运行此程序（奇怪的是，CLocus_12706_Sample_46_Locus_34641_Allele_0 [JoJo_s115.fq; groupI, 125578, +] 的序列为空）：

./extract_pop.py test.fa "JoJo"
>CLocus_12706_Sample_44_Locus_36326_Allele_0 [JoJo_s113.fq; groupI, 125578, +]
TGCAGCATGCTGGTGAACGCGTCATCATAAGCCTGTTGGCGAGCCAGCAGAAGGCGGCATGGGCAGCACTTAATAGGACGCACGTCCTCTGTGTCA
>CLocus_12706_Sample_46_Locus_34641_Allele_0 [JoJo_s115.fq; groupI, 125578, +]

如果没有python3，可以使用Biopython的SeqIO模块：

#!/usr/bin/env python

import sys
from Bio import SeqIO

fasta_file = sys.argv[1]
pop_name = sys.argv[2]

def get_pop(header):
    """Extracts the population name from the fasta header."""
    return header.split(" ")[1].split("_")[0][1:]

for seq in SeqIO.parse(fasta_file, format="fasta"):
    if get_pop(seq.description) == pop_name:
        print(">%s\n%s" % (seq.description, seq.seq))

sys.exit(0)