fastaファイルからシーケンスを抽出する

fastaファイルからシーケンスを抽出する

次のように、さまざまな長さの数十万のDNA配列を含むfastaファイル(正しい形式ではない)があります。

>NODE_213384_length_62_cov_8686_ID_2134025ATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAACGGAATCGAATGG>NODE_213385_length_62_cov_7933_ID_2134027ATCATCATCGAATGGAATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAAC>NODE_213386_length_62_cov_7184_ID_2134029AATGATTATTCCATTCGAGTCCATTCGATGATTCCATTCGATTCCATTCGATGATGATTGCA>NODE_213387_length_62_cov_8639_ID_2134031CAGAGCAGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAG>NODE_213388_length_62_cov_6833_ID_2134033AGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAGGTCAAT

単純なLinuxコマンドを使用して1000bpより長いシーケンスを抽出し、次のように正しいfasta形式で出力したいと思います。

>NODE_213384_length_62_cov_8686_ID_2134025
ATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAACGGAATCGAATGG
>NODE_213385_length_62_cov_7933_ID_2134027
ATCATCATCGAATGGAATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAAC
>NODE_213386_length_62_cov_7184_ID_2134029
AATGATTATTCCATTCGAGTCCATTCGATGATTCCATTCGATTCCATTCGATGATGATTGCA
>NODE_213387_length_62_cov_8639_ID_2134031
CAGAGCAGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAG

助けてくれる人に感謝します。

ベストアンサー1

線が長いようですね。 sedとawkはこの操作を処理するために利用可能なメモリを使用するため、問題が発生する可能性があります(もちろん、ファイルサイズ/行の長さによって異なります)。したがって、2段階のアプローチが採用されていますが、メモリの制限により、tr上記のawk、perl、またはsedソリューションのいずれかを使用できます。

head -20 inputfile | 
tr '>' '\n'  > stage1
perl -ne 'print ">$1 $2\n" if /^(.*?)([ACGTU]+)$/ && length($2)>1000' < stage1 > output

最初の20行の外観が満足のいくものであれば、実際には単一のパイプラインでこれを行うことができます。

tr '>' '\n' inputfile | 
perl -ne 'print ">$1 $2\n" if /^(.*?)([ACGTU]+)$/ && length($2)>1000' > output

私のPerlスクリプトは他のスクリプトほど効率的ではありませんが、作業を完了する必要があります。明確にするために書きました。 1000以上の塩基対を含む行がある場合にのみラベルを印刷し、その後にスペース、関連する塩基対、および改行文字を印刷します。

おすすめ記事