読者です 読者をやめる 読者になる 読者になる

CORDEA blog

Programming及びFedora21等のLinux OSのことが多めです。

GC含量を計算する

Bash

はじめに

過去の記事にて、ほぼawkGC含量の計算をしていたのですが、久々に見たらあまりにひどかったので書き直しました。

multifasta形式のファイルを投げると前回とは違ってGC含量だけ計算します。
もちろん一行書き換えるだけでもGC以外の計算は出来ます。
あとATGC以外の文字(Nとか)を予期していないので予期しない文字が入った場合の動作は保証しません。*1

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

 
 
Python, awk, Bash版があります。
Pythonawkが早いかと思います。

 

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

動作確認はDebianCentOSしかしていません。
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

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

*1:修正は簡単かと思います。