GC含量を計算する

はじめに

multifasta形式のファイルを投げるとGC含量を計算します。
ATGC以外の文字(Nとか) が入った場合の動作は保証していません。

　
動作確認にはNucleic and/or Amino Acid contentsを使用させていただきました。

2018/7/16 追記

今でもたまに参照されているようなので、少し修正を加えたものを以下に置いておきました

github.com

Mathematica (Wolfram)

fasta := URLDownload["http://example.com/hoge.fasta"]
lines = Import[fasta, {"FASTA", "LabeledData"}]
gcContent = Map[Module[{map = Counts[Characters[Last[#]]]}, First[#] -> Interpreter["PercentFraction"][N[(map["G"] + map["C"]) / Total[map]]]] &, lines]
Do[Echo[Last[x], First[x] <> ": "], {x, prc}]

Python

使いかた

python gc.py sequence.ffn

コード

#!/bin/env python
# encoding:utf-8
#
# Author:   Yoshihiro Tanaka
# Created:  2014-10-09
#
import sys
infile = open(sys.argv[1], "r")
lines = infile.readlines()
infile.close()
atgc = ['A', 'T', 'G', 'C']
fflag = False
lst = []
for l in range(len(lines)):
    line = lines[l]
    if ">" in line:
        if fflag:
            sys.stdout.write(seq)
            print("gc%: " + str( (lst[2]+lst[3]) / float(sum(lst)) ))
            lst = []
        seq = line
        fflag = True
    else:
        for i in range(len(atgc)):
            try:
                lst[i] += line.count(atgc[i])
            except:
                lst.append(line.count(atgc[i]))
sys.stdout.write(seq)
print("gc%: " + str( (lst[2]+lst[3]) / float(sum(lst)) ))

awk

awk version 20070501
　
動作確認はMac OS Xのみ

使い方

awk -f gc.awk sequence.ffn

コード

#!/bin/awk -f
# encoding:utf-8
#
# Author:   Yoshihiro Tanaka
# Created:  2014-10-26
#

function calc_gc(atgc) {
    sum = 0;
    for (i in atgc) {
        sum += atgc[i]
    }
    gc = ((atgc["G"] + atgc["C"]) / sum)*100
    return gc
}
function plus_atgc(atgc) {
    if (substr($0, i, 1) in atgc) {
        return substr($0, i, 1)
    }
}
{
    if ($0~/^>/) {
        if (length(atgc) > 0) {
            printf("%s%0.2f\n","gc%: ",calc_gc(atgc))
        }
        print $0
        for (i=1; i<=4; ++i) {
            atgc[substr("ATGC", i, 1)] = 0
        }
    }
    else {
        for (i=1; i<=length($0); ++i) {
            ++atgc[plus_atgc(atgc)]
        }
    }
}
END {
    printf("%s%0.2f\n","gc%: ",calc_gc(atgc))
}

Bash

動作確認はDebianとCentOSしかしていません。
GNU sed しか動かないようです。BSD sed(Mac OSとか)の方はgsed使うか書き換えて下さい。
　

使い方

chmod u+x gc.sh
./gc.sh sequence.ffn

コード

#!/bin/bash
#
# Author: Yoshihiro Tanaka
# date: 2014-10-08

flag=false
lines=$(wc -l $1 | sed -e 's/^  *//g' | cut -d ' '  -f 1)
j=1
for line in `cat $1`;
    do
        if [ `echo $line | grep ">"` ]; then
            if $flag; then
                echo $seq
                for val in ${array[@]}; do
                    let sum+=$val
                done
                echo "gc%: $(echo "scale=4; ((${array[2]}+${array[3]}) / $sum) * 100" | bc | cut -c 1-5)"
                sum=0
                array=()
            fi
            seq=$line
            flag=true
        else
            tmp=($(echo $line | sed -e 's/\(.\)/\1\n/g' | grep "[ATGC]" | sort | uniq -c | sed -e 's/^  *//g' | cut -d ' ' -f 1))
            i=1
            for val in ${tmp[@]}; do
                let result=${array[$i]}+$val
                array[$i]=$result
                let i++
            done
            if [ $lines -eq $j ]; then
                for val in ${array[@]}; do
                    let sum+=$val
                done
                echo $seq
                echo "gc%: $(echo "scale=4; ((${array[2]}+${array[3]}) / $sum) * 100" | bc | cut -c 1-5)"
            fi
        fi
        let j++
    done

軽く説明

主に

tmp=($(echo $line | sed -e 's/\(.\)/\1\n/g' | grep "[ATGC]" | sort | uniq -c | sed -e 's/^  *//g' | cut -d ' ' -f 1))

このあたりについて。
　

ATGCで並んだ文字列をsedで一文字ずつ改行してgrepでATGCだけ取り出します。
これは放置していると改行コードまで入ってしまうためです。

次のsortとuniqによってATGCの並び順はA, C, G, Tになります。
その後のsedでは前半のスペースを全て消しています。
ちなみにawk '{print $1}' によって以下の部分を省略できます。

sed -e 's/^  *//g' | cut -d ' ' -f 1

今回はawkを使わずに書きたかったのであえてこの書き方にしてあります。

　
そのため、以下の部分では　array[2] (C)　と　array[3] (G)　を取り出して計算しています。

echo "gc%: $(echo "scale=4; ((${array[2]}+${array[3]}) / $sum) * 100" | bc | cut -c 1-5)"

それと、以下の部分でcutしているのは、bcを行った際にどうしても5.5500のように00が付いてしまうのが紛らわしいと感じたためです。

"scale=4; ((${array[2]}+${array[3]}) / $sum) * 100" | bc | cut -c 1-5

何かいい方法があればご教授ください。