GRCh37とHg19

突然ですがヒトのゲノムにマッピングをするとき、NCBIの配列を使いますか？　それともUCSCの配列を使いますか？

GRCh37 （Genome Reference Consortium human build 37） ？　
Hg19 (UCSC human genome 19) ？

リファレンス配列、1番から22番+X、Y染色体に関して言えば、配列はどちらも同じです。
そう言われても信じない方もいると思うし、現にヒューマンエラーで多少違うんでは？と思っている方もいるでしょう。

どちらも全く同じ配列だ、ということはUCSCのゲノムブラウザなどで、GRCh37/hg19　と記述されていることから何となく信じている方がほとんどだと思います。
でも、たまにお客さんから、どっちがいいの？　本当に同じ配列なの？　と聞かれることもありました。

そういう疑問を持つ方もいるらしく、私はとりあえず、A、T、C、GそれとNの塩基の数がGRCh37とhg19で同じかどうか、染色体ごとに調べてみました。
「とりあえず」というのは、本当は塩基の並び順も調べるべきなのでしょうが、4つの塩基+Nの数が全く同じなら「合っている」としてしまおう、という大雑把な試みだからです。


GRCh37の配列は
ftp://ftp.ncbi.nih.gov/genomes/H_sapiens/Assembled_chromosomes/seq/
から、hs_ref_GRCh37.p5_chr1.fa.gzなどを1番から22番+XとY
Hg19の配列は
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/　
から、chr1.fa.gz などを1番から22番+XとY

解凍してfastaファイルにしたら、以下のようなコマンドを実行して塩基の種類と数を得ました。
参考にしたのはこのサイトです。

# UCSCのchr1の場合

tail -n+2 chr1.fa | awk '{ for ( i=1; i<=length; i++ ) arr[substr($0, i, 1)]++ }END{ for ( i in arr ) { print i, arr[i] } }' > hg19_01.txt
A 32485284
N 23970000
C 23064132
T 32559153
G 23070958
a 33085607
c 23960280
g 23945604
t 33109603

# GRCh37のhs_ref_GRCh37.p2_chr1の場合（GRCh37.p5ではなくて.p2でしたね、私が落とした時は。でも同じです。）

tail -n+2 hs_ref_GRCh37.p2_chr1.fa | awk '{ for ( i=1; i<=length; i++ ) arr[substr($0, i, 1)]++ }END{ for ( i in arr ) { print i, arr[i] } }' > GRCh37_01.txt
A 65570891
N 23970000
C 47024412
G 47016562
T 65668756

UCSCのは、塩基の大文字と小文字に分かれていますが、足すとGRCの数と一致することがわかります。

私はこれを、1番から22番+X、Yまで順番に実行して足し算して比べました。
（暇人ではないですよ）
結果、A、T、C、G、Nの数は、3番染色体を除く全ての染色体でUCSCのとGRCのとが一致しました。　
唯一つ、3番染色体は、Nの数が違ったのです。
正確には、
60,830,534番目の塩基が、UCSCではN、GRCではM
60,830,763番目と次の塩基が、UCSCではNN、GRCではRR
でした。
以下、ViewerはCLCのGenomics Workbenchを使用

上段がUCSC hg19、下段がGRCh37の3番染色体
NとMが違います。

こちらはNNとRR

ちなみに、これら塩基のある場所は、FHITという遺伝子のイントロン領域でした。

別に、どうってことは無いんでしょうけど、全く同じだと思っていた予想は外れてしまいました。

30億分の3塩基ですからね。

あと個人的に、Nという塩基がどれくらいあるのか気になったので、染色体ごとにグラフにしてみました。　これはGRCh37を基準にしています。

縦軸は塩基の数です。

全体を揃えて%で見るとこんな感じ。

へえー、Yって半分以上がまだNなんだー。





メチレーションをダイレクトに検出するシーケンサー

シーケンサーは塩基配列を読むもの

というのはSanger型から始まってSolexa型、454型、SOLiD型、Ion Torrent型等、今までのシーケンサーの常識です。
でも、配列だけではなく、メチレーションの有無も同時に読んでくれたら・・・
これがPacBioではできる（らしい）のです。　

Flusberg et al. Direct detection of DNA methylation during single-molecule, real-time sequencing.
Nat Methods. 2010 Jun;7(6):461-5.
昨年のthe Advances in Genome Biology and Technology (AGBT)のポスター（実際には行っていないんですけど）から注目していたのですが、PacBioのSingle Molecule Real Time （SMRT）シーケンスを使えば、塩基の修飾も配列と同時に検出することができるようなのです。

SMRTテクノロジーについては以前、
PacBio　データの特徴　1　と　2　に書きました。

SMRTからの塩基配列データは、パルス　（4種類の異なる光）として出力されます。

（AGBTポスターより）

パルスから得られる情報には、塩基そのものの他にも、一塩基が読みとられる時のパルスの幅 Pulse width (PW) と、次の塩基が読みとられるまでのパルスの間隔 Interpulse duration (IPD) もあります。
このIPDが、塩基が修飾されていると変化するのです。
具体的には、Aがメチル化されている　N6-methyladenosine (mA)　の場合、Tと結合する部分のNにメチル化が起こっているため、Tと結合するのが遅くなるそうです。
そうすると、メチル化されたAのパルスの間隔　IPDが、メチル化されていないAと比べて大きくなります。　（下図の上がメチル化A、下がコントロールA。　メチル化Aの後、次のTまでの間隔が、コントロールのAからTまでの間隔と比べて大きいのがわかります）
この原理を用いて、メチル化されているアデノシン（に限らずシトシンも）を見つけよう！　ということです。

でもこの検出方法は、うまくいくのでしょうか？　

論文ではまず、人工的に作ったメチレーション配列とコントロール配列199塩基を、両方とも何百回も同じところを読んで、平均のメチル配列のIPDと、非メチル配列のIPDを比較しています。

これを見ると、必ずしもメチル化された箇所のIPDが大きく変化するわけではなく、その周辺の塩基のIPDがより影響を受けることがわかりました。

縦軸が塩基ごとのIPDの比（メチルサンプル/ 非メチルサンプル）の平均で、横軸が塩基のポジション、▲印がメチル化の場所です。

上からN6-methyladenosine、5-methylcytosine、5-hydroxymethylcytosine　の順で、それぞれN=346、504、393　だそうです。

どうでしょう？　メチル化された塩基の場所のIPDが変化しているものもありますが、その前後で大きく変化を示しているパターンのものもあるのがわかります。

N6-methyladenosineでは、メチル化塩基とその5番めの塩基

5-methylcytosineでは、2番、3番、6番めの塩基

5-hydroxymethylcytosineでは、2番、6番めの塩基

のIPDが大きく変化しています。　（縦軸の値が異なるのに注意）

このようなパターンがあることは、メチル化の種類も見分けることができる、ということです。

既存のバイサルファイト法のような形では不可能でしょう。

もちろん何百回も読まなくても、数回、しかし繰り返し同じ所を読むことで、メチル塩基と非メチル塩基を区別できるそうです。　

SMRTのCircular topology型DNAテンプレートなら可能なのです。

赤線がメチルAを含むテンプレート、青線がコントロールです。

5回読めばはっきりメチルテンプレートとコントロールの、メチルAのパルスが区別できることを示しています。

人工配列はもういい、という方のために、彼らはDNA adenosine methyltransferase positive (dam+)の大腸菌から取り出したC. elegans fosmid　配列をも読んでいます。

比較対象は、Whole Genome Amplification（WGA）したコントロール配列です。

これは3.7kbの配列にフォーカスして解析しています。

dam+ のIPDをコントロール（WGA）のIPDで割った、IPD Ratio　が縦軸です。

Aメチル化はGATC配列に起こりやすいので、non-GATCの塩基（○）とGATC塩基（■）を分けています。

N=106で、平均IPDをプロットしています。

全体的に、non-GATC塩基はIPD ratio が2以下なのに対し、GATC塩基（メチル化塩基）は、IPD Ratioが3から7の間にいるので、区別できる！　ということは言えるかな？

ちょっと2500塩基付近の下のGATCが微妙ですが。

実際、これからの技術ですね。

IPDの変化のパターンは、メチル化の種類にだけ依存しているのか、周りの塩基によっても違うのか、生物種によっても違うのか、光学的な条件によっても違ってくるのか。

でも、希望は持てます。　一度パターンがわかってしまえば、シーケンスするだけでメチル化の種類もわかってしまうのですから。

Optical Mapping って

RNA-Seqの続きを書く予定でしたが、ちょっと寄り道して、別なことを。

只今、CBI学会（Chem－Bio Informatics学会）に来ています。
日立ソリューションズのブース＆プレゼンで、面白い技術を見つけました。
Optical Mapping　というものです。

アメリカでは数年前から導入されていたようですが、このたびマシーンが完成して日本では初めて発表されました。
その名も「全ゲノムソリューション」　
外挿遺伝子やリピート領域の同定、ゲノム構造の解析、高精度な菌株同一性に威力を発揮するといいます。

OpGEN 社のテクノロジーは、簡単に言うと全ゲノムレベルで制限酵素のマップを作ることにあります。
測定原理は、

1.  菌体から長い（200kb以上必要だそうです）ゲノムDNAを抽出し
2. ガラス基板上に並べて固定し、 
3. 制限酵素（1種類）でDNAを切断し
4. 蛍光により切断断片の長さを測定
5. 制限酵素マップをもとにアセンブルして、全ゲノムの制限酵素マップを作る

というものです。
詳しくはこちらのページ
DNAを基板に固定してから制限酵素で切っているので、ばらばらにはならず、Cuttingの位置情報が得られる、というのが重要な点です。

制限酵素のマップを作って、これが次世代シーケンサーの解析にどう役立つかというと、今までのショートリードによるde novo アセンブルは、どうしてもうまく読めないところや、リピートや外からの挿入遺伝子配列などがちゃんとつながらない、という問題がありました。
コンティグはできてもそれがつながらないということです。

それを補正する形で、あらかじめ制限酵素マップという物理的な地図を作り、これをもとにしてアセンブルしたときのコンティグをアライメントしていく、というのがこの技術の良いところです。

Optical Mappingでは、制限酵素のマップはできますが、これには塩基情報はありません。
あるのは一つの制限酵素が切った場所、の情報です。
では、どうやってこの制限酵素マップと、NGSからの配列情報、de novo Assemblyで出来たContigをアラインさせるのでしょうか？

答えは単純でした。
ある一定の長さを持つContigの中の、制限酵素Cutting部位を検索し、Contigひとつひとつに制限酵素マップをつくります。　これは付属のソフトで出来るそうです。
このContig内の制限酵素マップを、Optical Mappingで求めた全ゲノムマップに対してアラインしていくと、Contigに含まれる配列情報が全ゲノムマップに反映される、ことになります。

普通のNGSの感覚では塩基配列を基にアラインすると思ってしまうのですが、これは制限酵素のマップ（Cuttingの場所情報）を基にアラインするのです。　へぇー

実例として、あの大腸菌O-114のゲノムを紹介していました。
PacBio、Ion Torrent、Illumina　等で読まれたデータのContigは、すべてOptical Mappingで作った制限酵素マップがカバーしていました。

現在はバクテリアサイズのゲノムに対応しているそうですが、本国ではヒトゲノムサイズのマップにも成功しているそうです。
最近のモントリオールの人類遺伝学会でも発表したそうですが、普通の（遺伝病などを持たないという意味でしょうが）白人男性の全ゲノムをこの制限酵素マッピング技術で読んで、既知のゲノム配列の制限酵素マップと比較した結果、なんと

• 22番長鎖に240kbのInsertion
• 8番短鎖に4.5MbのInversion

など、新しいゲノム変異が見つかったそうです。
こんなに大きな変異が普通にあったとしたら、今のリシークエンスに使っているUCSCのHG19とかの既知ゲノム配列って、何なんだろう、って思いました。
我々が「既知だ」、と思っているヒトゲノムのリファレンス配列が、人によって数メガ塩基の単位で異なるとしたら、 このリファレンス配列にReadをMappingしているだけでは、染色体の構造変化を見つけることはできません。
ドラフト配列だらけの非モデル生物ゲノムについては、言うまでもありません。

改めて、染色体構造を物理的に知ることの重要性を確認しました。

なお、この技術やマシンについては日立ソリューションズさんにお問い合わせください。
きっと親身に説明してくれると思います。

RNA-Seqのカウントって....どうする？

“ There are known knowns; there are things we know we know.
We also know there are known unknowns; that is to say we know there are some things we do not know.
But there are also unknown unknowns -- the ones we don't know we don't know."

これは2002年 2月12日、記者から、イラクにあるとされる大量破壊兵器がテロリストの手に渡ったという確かな証拠が無いことを指摘されて、これに答えたドナルド・ラムズフェルド国防長官（当時）の言葉です。

当時、あまりに意味難解な返答に、たくさんのコメディアンがジョークのネタにしていたことを覚えています。

既知なこと、知っていること、わかっていること。　これらが存在することを知っていること、わかっていること。

わからないことがあるということを知っていること。　あるいは知らないこと。

・・・哲学的ですね。　

これとRNA-Seqとが何の関係があるかというと、ノーマライズの方法について、知っていることと知っていると仮定していること、知らないこと、をちゃんと知っておく必要があるからです。

タグの数から発現量を算出する時の、RPKMという方法を以前書きました。

RPKMでは、タグカウントを、サンプル内の全カウントで割って、さらに転写産物の長さで割って、定数を乗じて、ノーマライズするのでした。
くわしくはWikiなどを参照。

ここで、
タグ（リード）数はKnownですか？
サンプル内総タグ数はKnownですか？
転写産物の長さはKnownですか？

Isoformの存在、未知Exonの存在などを加味すると、転写産物の長さは、厳密にはUnknownですね。　しかしこれをKnownと仮定して、計算します。
タグ数も、この仮定の転写産物の長さにマップされた数としてカウントされるのですから、厳密にはUnknownです。
こうなってくると、サンプル内総タグ数は、転写産物にMappableなタグという意味では、Knownとするべきでしょうか？

とは言っても、転写産物の長さを厳密に認識するのは、ショートリードでは不可能ですから、みなさんRefSeqなどで公開されている既知のmRNAの長さを用いるのです。　
でも頭のどこかに、これは仮定の長さかもしれない、と覚えておくのは意味があるでしょう。

さて、
複数サンプル間で発現量を比較したいとき、遺伝子の長さがとても影響するということがわかっています。
Oshlack et al., Transcript length bias in RNA-seq data confounds systems biology. Biology Direct 2009, 4:14.
 
転写産物が長い遺伝子ほど、異なるサンプル間で発現量の変動が大きく計算される傾向があるのです。
この論文では遺伝子の長さで補正していません。　発現変動の大きさと、遺伝子の長さをそれぞれ軸にプロットすると、発現変動の大きい遺伝子ほどmRNAが長い遺伝子だということがわかったそうです。

aとcがRNA-Seqのデータで、縦軸があるしきい値で発現変動ありとされた転写産物の割合（％）、横軸がその転写産物の長さ（bp）です。　bは同じくマイクロアレイのデータです。

RNA-Seqの「発現変動ありとする」やりかたに問題があるのは、長い遺伝子ほど変動ありとされている遺伝子の割合が高いことから、明らかです。

そこで、RPKMのように転写産物の長さで補正することが必要になるのですが、これまた、RPKMも万能ではありません。

つづく