2013年 06月 24日

awkとsedで塩基配列をアミノ酸配列にする

手元のマルチファスタな塩基配列をアミノ酸配列に変えたい。
でも、色んなパッケージを取ってくるのはよくわかんないし依存関係とか考えてるより培養したい...

awk とかsed とか誰のmac にでも入ってるものだけでなんとかしてみました。
どんな環境でも動きますけど死ぬほど遅いです。
ひとつの配列を処理するのに1秒ぐらいかかります。

使う時は、#########に囲われてる部分をコピーして"nucl2prot.sh" という名前で保存。
手元のマルチファスタな塩基配列を"transcripts.fasta" という名前にして

sh nucl2prot.sh transcripts.fasta proteins.fasta

と打てば、"proteins.fasta" という名前のマルチファスタなアミノ酸配列が手に入ります。

################################################################################################

#nucl2prot.sh

#塩基配列をアミノ酸に。コドンテーブルを目的に合わせて書き直せ
#数字の3桁"," 区切りを利用しているので合わない時は考えて
#すごく遅い
#$1=[input.file] $2=[output.file]

more $1 | grep -v ">" | tr "ATCGN" "12345" | sed -e :a -e 's/\(.*[0-9]\)\([0-9]\{3\}\)/\1,\2/;ta' | tr "12345" "ATCGN" | awk '{print $1","}' | sed -e 's/TCA,/S/g' -e 's/TCC,/S/g' -e 's/TCG,/S/g' -e 's/TCT,/S/g' -e 's/TTC,/F/g' -e 's/TTT,/F/g' -e 's/TTA,/L/g' -e 's/TTG,/L/g' -e 's/TAC,/Y/g' -e 's/TAT,/Y/g' -e 's/TAA,/*/g' -e 's/TAG,/*/g' -e 's/TGC,/C/g' -e 's/TGT,/C/g' -e 's/TGA,/*/g' -e 's/TGG,/W/g' -e 's/CTA,/L/g' -e 's/CTC,/L/g' -e 's/CTG,/L/g' -e 's/CTT,/L/g' -e 's/CCA,/P/g' -e 's/CCC,/P/g' -e 's/CCG,/P/g' -e 's/CCT,/P/g' -e 's/CAC,/H/g' -e 's/CAT,/H/g' -e 's/CAA,/Q/g' -e 's/CAG,/Q/g' -e 's/CGA,/R/g' -e 's/CGC,/R/g' -e 's/CGG,/R/g' -e 's/CGT,/R/g' -e 's/ATA,/I/g' -e 's/ATC,/I/g' -e 's/ATT,/I/g' -e 's/ATG,/M/g' -e 's/ACA,/T/g' -e 's/ACC,/T/g' -e 's/ACG,/T/g' -e 's/ACT,/T/g' -e 's/AAC,/N/g' -e 's/AAT,/N/g' -e 's/AAA,/K/g' -e 's/AAG,/K/g' -e 's/AGC,/S/g' -e 's/AGT,/S/g' -e 's/AGA,/R/g' -e 's/AGG,/R/g' -e 's/GTA,/V/g' -e 's/GTC,/V/g' -e 's/GTG,/V/g' -e 's/GTT,/V/g' -e 's/GCA,/A/g' -e 's/GCC,/A/g' -e 's/GCG,/A/g' -e 's/GCT,/A/g' -e 's/GAC,/D/g' -e 's/GAT,/D/g' -e 's/GAA,/E/g' -e 's/GAG,/E/g' -e 's/GGA,/G/g' -e 's/GGC,/G/g' -e 's/GGG,/G/g' -e 's/GGT,/G/g' -e 's/---,/-/g' -e 's/NNN,/N/g' -e 's/N..,/N/g' -e 's/.N.,/N/g' -e 's/..N,/N/g' > temp_nucl2prot_seq.txt

more $1 | grep ">" > temp_nucl2prot_head.txt

paste temp_nucl2prot_head.txt temp_nucl2prot_seq.txt | awk '{print $1"\n"$2}' > $2

rm temp_nucl2prot_head.txt
rm temp_nucl2prot_seq.txt

echo "(^q^)/"


################################################################################################

以上です。
ほんとに遅いから心配になったら

wc -l temp_nucl2prot_seq.txt

とかやって進行状況を見ても安心して下さい。

#20140109追記
塩基配列の文字数が3の倍数になってないとダメです。
その場合は生成されたアミノ酸配列に「,」が入ってきます。
数字の桁区切りの「,」を利用して3文字を区切っているため、例えば4文字であれば
ATCG→1234→1,234→A,S
と残念な結果になります。
5文字だと、
AATCG→11234→11,234→AA,S

人気ブログランキングに参加しています。
応援よろしくお願いします。
FC2 Blog Ranking

[PR]

by koretoki | 2013-06-24 21:31


<< ホルマリン固定パラフィン包埋な...      ゲノムアノテーションファイルが... >>