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 

集団は、例えば、「LakeCamp」または「ジョジョ」または「NaknekRiver」を含むIDの一部です。

私はこの記事の後に、どのように配列を抽出するかを調べようとしました。 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)