☰

リストからfastaエントリを抽出するには、読みながら使用してください。

tag-icon • tag-icon text-processing sed bioinformatics
リストからfastaエントリを抽出するには、読みながら使用してください。

私はそれぞれ約14,000個の「アイテム」を含む28個のファイルを持っています。単一の項目は、ヘッダー（>文字列で表される）、改行文字、および一連の文字列で構成されています。各項目には可変長シーケンス/文字列があります。 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

私は、gene.txtという項目のヘッダーにそれぞれ対応する14,000行で構成される別々のテキストファイルを持っています。ファイルの最初の数行は次のとおりです。

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個のファイルを生成しますが、各ファイルには1つのシーケンスしかありません。

>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

シェルは、この種の解析には実際には適していません。コードからファイル全体を一度読むことがわかります。各ファイルから読み取った遺伝子名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ときに遺伝子名から出力を生成します。 filenameをアンダーgenename.faスコアの前に現在のファイル名部分を含めるようにヘッダーを再構築します。

元のヘッダーがXXX_ref.fasta次のような場合

>genename

これは次のように変換されます。

>XXX_genename

スクリプトの最後の部分は、awkすべての行を適切な出力ファイルに送信します。

提供したデータでテストすると、次の3つのファイルが生成されます。

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

$ cat FBgn0076379.fa
>CR1_FBgn0076379
ATGCTGCGCACCCTTTTCGCCGTGCGTGGTCAGTGCCAGCAGCTGCTGAGGAGAACATTCACCCCCCATTGCAGTGGCCAACGA

関連情報