ヒメツリガネゴケゲノムの制限酵素認識配列を数える

西山智明 tomoakin@staff.kanazawa-u.ac.jp

2014年11月4日

要点:

unfoldした配列に対して制限酵素の認識配列でgrepして-o オプションを使い認識配列を出力

前提

標準的unix shell tools (grep, wc)
ruby
fatt (https://github.com/mkasa/klab)が使える状態になっている

参照配列を途中に改行の入らないテキストにする。

fatt unfoldを利用して一行。

$ fatt unfold Physcomitrella_patens.1_1.fasta  >  Pp.ufa

これで行をまたがるサイトの問題を気にせずgrepで探せるようになる。

grepしてみる

$ grep -c GACGTC  Pp.ufa
924

行を書き出さず数だけを示す-cオプションをつけて実行してみると、すぐ終わって大体scaffold総数より少し少ないくらいだからおそらく大丈夫だろう。でもこのままではサイトの数は分からない。

-o オプションとwcを組み合わせる

grepのマニュアルをさらに読んで-oオプションを使うとマッチした部分を一行１マッチで書き出せる事がわかる。その行数を数えればマッチした数だろう。

$ grep -o GACGTC  Pp.ufa | wc
  27527   27527  192689

大体良さそうである。無駄な出力があるけど-lオプションをつけてやれば行数だけになるはずなのでそこは簡単。

参照配列の読み込み

上記がコマンドラインレベルで確認できれば、認識配列を表から読み込んで上記のコマンドの出力を付け足すrubyプログラムは以下のように書ける。タブ区切りテキストの２カラム目に認識配列が書いてあるデータを予定している。

#!/bin/env ruby
#
ARGF.each_line do |l|
  ar = l.chomp.split("\t")
  ar << `grep -o #{ar[1]} Pp.ufa | wc -l `
  puts ar.join("\t")
end

検索実行

上記プログラムが REcountgrep.rbというファイルに書かれていて制限酵素と認識配列のペアがrepairというファイルに格納されているとして

$ ruby REcountgrep.rb repair > REpaircounts

とすると結果がREpaircountsに書き出される。ここはそれなりに時間がかかる(約10分)。その間に、平均長を計算する方法を検討する。

scaffoldのgapでない部分の長さの和を求める。

fatt statでcontig total lengthを見れば概数としては妥当な値が得られるであろう。contig total lengthを制限酵素サイトの数で割れば、これまた端の数の誤差を除き概ね、平均長になるはずであろう。

$ fatt stat Pp.ufa
Scaffold (w/gap) statistics
Total # of bases = 479,985,347
# of scaffolds = 2,106
Max size = 5,387,533 (# = 1)
N50 scaffold size = 1,316,933 (# = 111)
N70 scaffold size = 837,602 (# = 203)
N80 scaffold size = 614,044 (# = 270)
N90 scaffold size = 369,512 (# = 369)
Min size = 1,000
Total scaffold # = 2,106
Avg size = 227,913

Scaffold (wo/gap) statistics
Total # of bases = 453,927,444
# of scaffolds = 2,106
Max size = 5,108,762 (# = 1)
N50 scaffold size = 1,281,184 (# = 108)
N70 scaffold size = 830,544 (# = 196)
N80 scaffold size = 623,949 (# = 258)
N90 scaffold size = 389,436 (# = 349)
Min size = 1,000
Total scaffold # = 2,106
Avg size = 215,540

Contig statistics
Total # of bases = 453,927,444
# of contigs = 21,041
Max size = 593,018 (# = 1)
N50 contig size = 72,478 (# = 1,846)
N70 contig size = 44,147 (# = 3,452)
N80 contig size = 31,680 (# = 4,660)
N90 contig size = 16,724 (# = 6,575)
Min size = 1
Total contig # = 21,041
Avg size = 21,573

total scaffoldの479 Mではなく、Contig statisticsのTotal # of bases の453,927,444の値を採用する。

平均長算出

この数を表の３カラム目の数で除すると平均長の概数が算出できる。rubyの簡単なプログラム

#!/usr/bin/env ruby
#
ARGF.each_line do |l|
  ar = l.chomp.split("\t")
  ar <<  453927444.0 / ar[2].to_i
  puts ar.join("\t")
end

を通して終わりでありましょう。

$ ruby REavlength.rb  REpaircounts 
Aat II	GACGTC	27527	16490.262069967666
Afl II	CTTAAG	115030	3946.1657306789534
Age I	ACCGGT	13002	34912.12459621597
Apa I	GGGCCC	12403	36598.19753285495
ApaL I	GTGCAC	51191	8867.329100818502
BamH I	GGATCC	66103	6866.971907477724
Bci VI	GGATAC	55246	8216.476197371756
Bgl II	AGATCT	196098	2314.798947465043
BspH I	TCATGA	175768	2582.5374584679807
BssH II	GCGCGC	8763	51800.46148579253
Bst1107 I	GTATAC	64606	7026.088041358388
Cla I	ATCGAT	73215	6199.924113911084
Dra I	TTTAAA	679157	668.3689397296943
Eci I	TCCGCC	24629	18430.607982459704
EcoR I	GAATTC	147014	3087.6477342293933
EcoR V	GATATC	132635	3422.380548120783
EcoT14 I	CCWWGG	0	Infinity
EcoT22 I	ATGCAT	182986	2480.6676139158185
Hinc II	GTYRAC	0	Infinity
Hind III	AAGCTT	299590	1515.1622016756235
Hpa I	GTTAAC	50481	8992.045403221015
Kpn I	GGTACC	35529	12776.251625432746
Mlu I	ACGCGT	14525	31251.4591394148
Nae I	GCCGGC	11316	40113.77200424178
Nco I	CCATGG	74919	6058.909542305691
Nde I	CATATG	135050	3361.1806293965196
Nhe I	GCTAGC	30871	14704.00842214376
Nru I	TCGCGA	18069	25121.890752116884
Pst I	CTGCAG	69531	6528.4181731889375
Pvu I	CGATCG	22279	20374.677678531352
Pvu II	CAGCTG	59532	7624.931868574884
Sac I	GAGCTC	65039	6979.31155153062
Sac II	CCGCGG	11001	41262.38014725934
Sal I	GTCGAC	33011	13750.793493078065
Sca I	AGTACT	93262	4867.228281615235
Sma I	CCCGGG	8232	55141.81778425656
Sml I	CTYRAG	0	Infinity
SnaB I	TACGTA	36834	12323.598957484932
Spe I	ACTAGT	84724	5357.719701619376
Sph I	GCATGC	60382	7517.595376105462
Spl I	CGTACG	9844	46112.093051605036
Ssp I	AATATT	738106	614.9895055723704
Stu I	AGGCCT	53884	8424.160121743003
Xba I	TCTAGA	194032	2339.4462975179354
Xho I	CTCGAG	55011	8251.575939357583
Xma I	CCCGGG	8232	55141.81778425656

これにて計算終了。なお、degenerateなパターンを持つサイトは手抜きによって対応していない。今回は、目標が5k-10kの断片サイズを作る酵素だったので切れ過ぎるはずのものを計算する必要はない。なお、degenerateなのに合わせて正規表現パターンを置き換えれば十分なはずである。メチル化されているサイトは切れない事があることも考慮していない。