SOAPdenovo2になって，アセンブリで得られた配列の統計情報がscafStatisticsというファイルに出力されるようになった．そのため，自分でアセンブリ結果を解析しなくとも，コンティグ数やN50などの基本情報をチェックすることができる……のだが，scafStatisticsは項目を並べただけのテキストファイルなので，複数のアセンブリ結果を比較しようと思った場合，いちいちファイルを開いて統計情報を集計するのが非常に面倒になる．
ということで，SOAPdenovo2のアセンブルで出力される情報を一つにまとめるスクリプトを簡単に書いてみた．タブ区切りテキストで出力されるので，Excelなどでも見やすいようにしてある．
インストール
上記のgistに置いてあるソースコードをダウンロードするか，wgetでソースコードを落としてくる．
1
$ wget --no-check-certificate https://raw.github.com/gist/3870629/6547d69e768e3d87a140d6405fdf19102ff525cb/soapdenovo2_stats.rb
動作確認はRuby 1.8.7および1.9.3で行っているが，まあ変なことはしてないのでどんな環境でもだいたい動くと思う．RubyのバージョンよりかはSOAPdenovo2の未知の部分や仕様変更などが怖いわけで，そもそも以下のデータセットで出てきたscafStatisticsを参考にしているため，SOAPdenovo2のパラメータなどによっては期待通りの動きをしないかもしれない．あと，例外処理などをあまり真面目に書いていないので，入力ファイルが見つからない場合にエラーを吐かなかったりするが，そこらへんは後々修正するかも．
使い方
基本的な使い方
ruby soapdenovo2_stats.rb
1
2
3
4
$ ruby soapdenovo2_stats.rb path/to/k23.scafStatistics path/to/k25.scafStatistics
Filename Size_includeN Size_withoutN Scaffold_Num Mean_Size Median_Size Longest_Seq Shortest_Seq Singleton_Num N50
k23.scafStatistics 4533843 4533843 845 5365 2919 55290 100 845 11457
k25.scafStatistics 4536449 4536449 654 6936 3609 71235 100 654 14798
その他の詳細な統計量を出力する
-a--all-c--contig
1
2
3
$ ruby soapdenovo2_stats.rb -a -c path/to/k23.scafStatistics
Size_includeN Size_withoutN Contig_Num Mean_Size Median_Size Longest_Seq Shortest_Seq Contig>100 Contig>100 ( %) Contig>500 Contig>500 ( %) Contig>1K Contig>1K ( %) Contig>10K Contig>10K ( %) Contig>100K Contig>100K ( %) Contig>1M Contig>1M ( %) Nucleotide_A Nucleotide_A ( %) Nucleotide_C Nucleotide_C ( %) Nucleotide_G Nucleotide_G ( %) Nucleotide_T Nucleotide_T ( %) GapContent_N GapContent_N ( %) Non_ACGTN Non_ACGTN ( %) GC_Content N10 Contigs >0 in N10 N20 Contigs >0 in N20 N30 Contigs >0 in N30 N40 Contigs >0 in N40 N50 Contigs >0 in N50 N60 Contigs >0 in N60 N70 Contigs >0 in N70 N80 Contigs >0 in N80 N90 Contigs >0 in N90 NG50 Contigs >0 in NG50 N50_contig-NG50_contig_length_difference Number_of_contigs_in_scaffolds Number_of_contigs_not_in_scaffolds( Singleton) Average_number_of_contigs_per_scaffold
4533843 4533843 845 5365 2919 55290 100 840 99.41 646 76.45 589 69.7 147 17.4 0 0.0 0 0.0 1122311 24.75 1154274 25.46 1146494 25.29 1110764 24.5 0 0.0 0 0.0 50.75 29554 13 21527 31 17094 55 13748 85 11457 121 9180 165 7136 2225257 296 3315 403 NaN NaN NaN 0 845 0
ワイルドカードを使って入力ファイルを複数指定する
Zshなどのシェルコマンドと同様に，入力ファイルの指定には*?
1
2
3
4
5
6
7
8
9
$ ruby soapdenovo2_stats.rb path/to/*.scafStatistics
Filename Size_includeN Size_withoutN Scaffold_Num Mean_Size Median_Size Longest_Seq Shortest_Seq Singleton_Num N50
k23.scafStatistics 4533843 4533843 845 5365 2919 55290 100 845 11457
k25.scafStatistics 4536449 4536449 654 6936 3609 71235 100 654 14798
k27.scafStatistics 4539473 4539473 578 7853 3781 79498 100 578 18157
k29.scafStatistics 4544164 4544164 570 7972 3672 103369 100 570 17945
k31.scafStatistics 4550461 4550461 621 7327 3306 77302 100 621 17052
k33.scafStatistics 4582773 4582773 1655 2769 1688 22953 100 1655 5435
k35.scafStatistics 2393501 2393501 16300 146 132 1345 105 16300 142
特定のカラムでソートする
-s
1
2
3
4
5
6
7
8
9
$ ruby soapdenovo2_stats.rb -s N50 path/to/*.scafStatistics
Filename Size_includeN Size_withoutN Scaffold_Num Mean_Size Median_Size Longest_Seq Shortest_Seq Singleton_Num N50
k27.scafStatistics 4539473 4539473 578 7853 3781 79498 100 578 18157
k29.scafStatistics 4544164 4544164 570 7972 3672 103369 100 570 17945
k31.scafStatistics 4550461 4550461 621 7327 3306 77302 100 621 17052
k25.scafStatistics 4536449 4536449 654 6936 3609 71235 100 654 14798
k23.scafStatistics 4533843 4533843 845 5365 2919 55290 100 845 11457
k33.scafStatistics 4582773 4582773 1655 2769 1688 22953 100 1655 5435
k35.scafStatistics 2393501 2393501 16300 146 132 1345 105 16300 142
今回用いたサンプルデータ
なお，今回上記で示したデータは，E.coliのSRR001665のデータを用いてSOAPdenovo2_revision210でアセンブルしたものを使用している．アセンブルのコマンドとconfigファイルは以下のとおり．SEQanswersのWikiにあったExampleをそのまま使わせてもらっている（ http://seqanswers.com/wiki/How-to/de_novo_assembly#SOAP_denovo ）．
1
$ SOAPdenovo-63mer all -K 25 -R -s cont.config -o output/k25
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#maximal read length
max_rd_len = 36
[ LIB]
#average insert size
avg_ins = 200
#if sequence needs to be reversed
reverse_seq = 0
#use for contig building only
asm_flags = 1
#in which order the reads are used while scaffolding
rank = 1
#fastq files
q1 = /home/yag_ays/tmp/soapdenovo2_testrun/SRR001665_1.fastq
q2 = /home/yag_ays/tmp/soapdenovo2_testrun/SRR001665_2.fastq
参考
ソースコード
#/usr/bin/env ruby
require "optparse"
def v(s)
if s == "NaN"
return s
elsif s.include?("%")
return s.gsub("%","").to_f
else
return s.to_i
end
end
def parse_scafstatistics(file)
stats = []
header = []
open(file) { |f|
header << "Filename"
stats << File.basename(file)
f.each_line do |line|
if !line.include?("<--") && line != "\n"
l = line.chomp.split("\t")
if l[0] == "GC_Content"
header << l[0]
stats << v(l[1])
elsif l.length == 3
if l[2].include?("%")
header << l[0].strip
header << l[0].strip + " (%)"
stats << v(l[1])
stats << v(l[2])
else
header << l[0]
header << "Contigs >0 in " + l[0]
stats << v(l[1])
stats << v(l[2])
end
else
header << l[0]
stats << v(l[1])
end
end
end
}
return header, stats
end
def print_tsv(header, stats, print_all, print_scaffold)
h = []
s = []
if print_scaffold
if print_all
h = header[0..57]
s = stats.map{|a| a[0..57]}
else
h = [header[0..8],header[45]]
s = stats.map{|a| [a[0..8],a[45]]}
end
else
if print_all
h = header[58..header.length]
s = stats.map{|a| a[58..a.length]}
else
h = [header[58..65].flatten,header[98]]
s = stats.map{|a| [a[58..65].flatten,a[98]]}
end
end
puts h.join("\t")
s.each do |f|
puts f.join("\t")
end
end
if __FILE__ == $PROGRAM_NAME
header = []
stats = []
sort_column = nil
print_all = false
print_scaffold = true
ARGV.options do |opt|
opt.on( "-a","--all") { print_all = true }
opt.on( "-s VAL","--sort") { |a| sort_column = a }
opt.on( "-c", "--contig") { print_scaffold = false }
opt.on( "-h","--help") { puts opt ;exit }
opt.parse!
end
Dir.glob(ARGV).each do |f|
h, s = parse_scafstatistics(f)
header = h
stats << s
end
if sort_column
i = header.index(sort_column)
if i == nil
puts "ERROR : unknown column name '#{sort_column}' "
exit 1
end
stats = stats.sort{|a,b| b[i] <=> a[i] }
end
print_tsv(header, stats, print_all, print_scaffold)
end