☰

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

tag-icon • tag-icon text-processing
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以上の塩基対を含む行がある場合にのみラベルを印刷し、その後にスペース、関連する塩基対、および改行文字を印刷します。

答え2

「簡単なLinuxコマンド」は次のとおりです。

sed 's/\(>[^ACGT]*_[0-9]\+\)\([ACGT]\+\)/\1\n\2\n/g' yourdnafile |egrep -B1 '^[ACGT]{1000}'

sed 部分は各グループを 2 つの行に分け、grep は 1000 を超える行とその前の行 (-B1) を表示して一致させます。

それとももっと簡単かもしれません：

sed 's/\(>[^>ACGT]*\)\([ACGT]\{1000\}[ACGT]*\)/\1\n\2\n/g;s/>[^>ACGT\n]*[ACGT]\+//g' yourdnafile

答え3

最初にターゲットフォーマットに分割し、DNA配列の十分に長い線対を印刷する2つのステップでこれを行うことができます。たとえば、1000bpDNA配列の文字の長さ（「length」という単語の後の値ではない）を参照すると、次のようになります。

cat inputfile | sed -e 's/>/\n>/g' -e 's/\(ID_[0-9]*\)/\n\1/g' |\
awk -F_ '$1=="NODE" { n=$0; next; } length($0)>1000 { print n; print ">" $0;}' \
> outputfile

関連情報