DNA配列のGC含量

GorCの塩基を多く含むDNA配列は比較的柔らかく、AorTの塩基を多く含む配列は硬いと言われている。
このことから、DNA配列のGC（AT）含量はヌクレオソーム形成などの物性の考察に有用であると考えられている。

ヌクレオソーム配列などの特定長配列でのGC含量計算

例えば、
「
DANPOS2(https://sites.google.com/site/danposdoc/ )でヌクレオソーム配置を予測
→ bedtools getfasta(http://kazumaxneo.hatenablog.com/entry/2017/07/23/213505 )でヌクレオソーム配列の取得
→ FASTA形式 147bp ヌクレオソーム配列群データ
」

この配列群の各場所（各塩基）でのATGC率を計算して、ヌクレオソーム親和的な配列特徴を探りたい。とか。

つまり

特定長（この例だと50bp）のFASTAファイルを

hoge_50bp.fa

>chr1:765821-765871
ttttctggaattctttgattgataattttttcttctcagtcttttatctt
>chr1:766223-766273
catgGTGCATCCTCCTGTTGACTCCTGACCATCTGTTCACATGgggctgc
>chr1:793081-793131
CGGCTTCCCAGCGAAGGTTTTTGAAGCACGGTTTGATTTTCTCTCTCCCC
>chr1:793465-793515
ACCTGCTGGTGGCAGTGATCTGACCGCTGTTTAAGCTTTCTTTGAACTCC
>chr1:1582152-1582202
tggttcttctctctgagcagcttgagtatgccactccactgccttctggc
・・・(配列数いくらでも)

こうしたい。（配列方向は縦　上流→下流）

A T C G
0.301 0.306 0.194 0.199 　　　　#１番目塩基
0.296 0.307 0.185 0.212 　　　　#２番目塩基
0.309 0.287 0.23 0.174 　　　　 #3番目塩基
0.329 0.289 0.188 0.194 　　　　#4番目塩基
0.317 0.31 0.183 0.19 　　　　  #5番目塩基
0.299 0.291 0.205 0.205 　　　　#6番目塩基
0.284 0.325 0.196 0.195 　　　　#7番目塩基
・・・(今回なら50行)

ワンライナー

awkで２次元リストを使う暴挙。便利です。

bash


cat hoge_50bp.fa | grep -v ">" | tr '[a-z]' '[A-Z]'| awk '{for(i=1;i<=length($0);i++){a[i,substr($0,i,1)]+=1}}END{OFS=" ";for(j=1;j<=length($0);j++){n=a[j,"A"]+a[j,"T"]+a[j,"C"]+a[j,"G"];print a[j,"A"]/n,a[j,"T"]/n,a[j,"C"]/n,a[j,"G"]/n}}'

初心者用コード解説

cat →　FASTAファイル全体表示
grep -v　→　">"のない行を検索
tr →　小文字を大文字に変換

awk内部
a[i, substr($0,i,1)]++　→　先頭から１塩基ずつ切り取って、その場所（i）のその塩基(substr($0,i,1))をカウント
n　→　ATGCの割合にするために、個数の合計を計算。print時にこれで割る。

以上。