塩基配列のAT含量（CGまで）をちゃちゃっと計算させるシェルコマンド

#AT含量計算に無駄な機能はいらない
AT含量計算なんて、そこら中にツールがあるので探した方が早い。

が、自分の場合、コピペした塩基配列とかFASTAからとかbedtools getfasta出力からAT含量をちゃちゃっと計算したい瞬間が多く、

１．ターミナルで動く
２．カスタムできる
３．後の描画用にATCGそれぞれで出力してほしい

を条件に加えると、作ったほうが早い。

atcg.bash

#!/bin/bash
set -e
if [ -p /dev/stdin ]; then
    allseq=$(cat -)
else
    allseq=$(cat $1)
fi

echo -e "$allseq" | tr '[a-z]' '[A-Z]' | awk '/^>/{print $0}!/^>/{for(i=1;i<=length;i++){a[substr($0,i,1)]++};n=a["A"]+a["T"]+a["C"]+a["G"];OFS="\t";print a["A"]/n,a["T"]/n,a["C"]/n,a["G"]/n}'

最初のif-elseは、標準入力受け取りとファイル受け取りに対応するため。
これをコピペしてファイル名「atcg」で保存して

bash

sudo ln -s /path/to/atcg /usr/local/bin/

とでもしてPATHつなげれば「atcg」コマンドとして動く。

##実行
出力はコマンド名通り、A T C Gの順でタブ区切り出力。
コピペからは

bash

echo tctgttcttcagaatggataatttctac | atcg
#0.25	0.428571	0.178571	0.142857

FASTAからは

bash

atcg FASTA.fa

#>CHR9:4001761-4001911
#0.306667	0.326667	0.22	0.146667
#>CHR21:45231589-45231739
#0.273333	0.233333	0.243333	0.25
#>CHR3:116878420-116878570
#0.268889	0.217778	0.248889	0.264444
#>CHR16:26967230-26967380
#0.268333	0.205	0.266667	0.26
#>CHR5:155754502-155754652
#0.258667	0.233333	0.253333	0.254667

配列名まで大文字になるのはご愛嬌。
bedtools getfasta出力はデフォルトでFASTA形式なので、パイプで同様に行ける。

bash

bedtools getfasta -fi /path/to/refgenome -bed /path/to/bedfile | atcg

全部つなげたければ

bash

grep -v "^>" FASTA.fa | tr -d '\n' | atcg
#0.258667	0.233333	0.253333	0.254667

各場所でそれぞれ出したいとか、Gnuplotでそのまま書きたいとかなら、以下参照。

【ATGC率計算】特定長DNA配列データから各塩基のATGC率を計算する。ヌクレオソーム配列解析など。（ https://qiita.com/Alzt/items/2da327030234c9971520 ）

以上。