はじめに
multifasta形式のファイルを投げるとGC含量を計算します。
ATGC以外の文字(Nとか) が入った場合の動作は保証していません。
動作確認にはNucleic and/or Amino Acid contentsを使用させていただきました。
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
何かいい方法があればご教授ください。