Exome と SNVについてのメモ 1

最近ブログを更新していませんでした。
先月末から何ていうのでしょう、一種のスランプみたいな、そんな気分でしたので。
もう脱しつつあるので、大丈夫です。

昔からGWASなどでSNPを解析しているひとには常識かもしれませんが、NGSが出た最近やっとSNP（SNV）を始めてみようかな、というひとには意外と「目からうろこ」なことがあります。

私もSNPを昔からやっていたわけではないので、たまに論文を読んでいて気づかされることがあります。
で、そういうことを他の研究者と話すと、意外とその人も知らなかったり。

具体的には、NGSを使ってWhole Exomeをやっている解析で、SNVを見つける、そのときのワークフローです。
当たり前のように行っていたフィルタリングのプロセスで、「なぜそのフィルタリングを行うのか」「基本となる考え方・仮定は何か」という基本的なところがおろそかになっていたと。
そこについてメモ書きのようですが、まとめました。

Exome解析の最終目的は、おおよそ、「疾患原因となるSNVの候補を絞り込む」ことになると思います。
そこで、

1. ヒトならヒトのExon領域をカバーしたターゲットシーケンスを行い、その領域で十分な厚みを持ってマップされた場所から、SNVを見つけます。　（例：SureSelectなどでターゲットキャプチャーしたあとIlluminaマシンで大量に読み、BWAでヒトゲノムにマッピング、キャプチャー領域にてSamtoolsでSNVを検出）
2. 検出されたSNVの場所とリストをもとに、dbSNPなどのデータベースに無いものを抽出
3. さらにその中からNon-synonymousのSNV（アミノ酸置換を伴うSNV）のみを抽出
4. そして残ったnsSNVのうち、患者複数のサンプルで共通するものを選び出す

ここから先はその疾患原因SNV候補のバリデーションになるわけですが、上記の１～４の処理はそれぞれ理由、というか前提や仮定があります。
それを理解していないでなんとなくパイプライン的に解析を処理していると、後から他の研究者に突っ込まれてタジタジ・・・となってしまう！

以下はSNPをやってたひとには当たり前すぎることを書いていると思います。　どうぞお許しを。

1. Exon領域をキャプチャーする大前提は、テーマとしている疾患の変異を、遺伝子のコーディング領域から探そうとしているわけです。　つまりコーディングされない場所にいくら疾患原因SNVがあったとしても、それは見ないことにしよう、としています。　
   もちろんコーディング領域は重要です。　しかしそれ以外の転写制御領域のSNVも大変重要ですが、これらは一般的なExomeでは解析対象から外されます。　
2. dbSNPに無い、SNVのみを探すというのは、一見、新規性を探索するのに合理的な方法です。　この疾患原因SNVは非常にレア、珍しいという前提があります。　しかし、データベースに登録されるSNVは、現在爆発的に増え続けていますので、既知であっても自分の研究疾患では新規ということは十分あり得ます。
3. Non-synonymousのSNVに絞るということは、アミノ酸置換、あるいはフレームシフト、ミッセンス等を伴うような変異は特に重要だという前提に立っています。　「この疾患の原因は、アミノ酸変異を伴うSNVだ」　という仮定がそこにあります。　そうでない仮定の場合にNon-synonymous SNVに絞るのは適当でないでしょう。
   また気をつけるべきは、Alternative Exon（スプライスバリアント）の存在です。　これを無視してSNVを見つけているケースが多いのではないでしょうか。　ゲノムのSNVではなく遺伝子コーディング領域中のSNVとなると、スプライスバリアントが異なればそのSNVの有る無しの意味は大きいはずですよね。　（言葉足らず）
4. 複数の患者のサンプルで同じSNVを見つけるというのは、そのフェノタイプに共通するSNVを見つけ出すということです。　SNV（原因）があれば必ずそのフェノタイプ（結果）を引き起こすことをComplete Penetranceと言います。　一方、そのフェノタイプ（結果）には必ずそのSNV（原因）が見られるときをComplete Detectance　と言います。
   つまりこのフィルタリングは、全ての患者に共通するSNVがあるはずだ、という前提に立っているわけです。

まとめると、上記のワークフローで解析するための前提・仮定は、
「全ての患者に共通するSNVがあって、かつそのSNVを持っていればかなりの確率でその疾患にかかり、そのSNVはとても珍しく、そしてタンパク質のアミノ酸コードを変えることで疾患フェノタイプを引き起こす」　
ということになりますでしょうか。

今日は字ばっかりですみませんでした。

SIFT SNP機能予測アルゴリズム

今日もSNVの話です。
１塩基の変異が、タンパク質にどんな影響を与え得るのか、の予測アルゴリズムです。
今までPolyPhenとかGrantham Scoreとかについて少し触れましたが、今日はSIFTというアルゴリズム＆ツールについて。

"SIFT"　ってググると、上位には画像検出アルゴリズムがたくさん出てきますが、これは違うので必ず、"SIFT　SNP"　とかで検索して下さいね。　（画像検出の方も興味がありますが）
http://sift.jcvi.org/

SIFTとは、Sorting Intolerant From Tolerant　の略です。
？？？
そもそもtolerant とは何でしょう？
辞書で引くと寛容な、とか耐性、とかがあります。　
ここではタンパク質の機能に変化が無い、ということを意味します。


SIFTの考え方の基本は、配列保存性です。
同じタンパクファミリーの中で、保存性の高いアミノ酸配列があったとします。　
そのような配列の中のアミノ酸変異・置換は、タンパク質の機能に致命的な影響を与える可能性が高いと思います。　
機能を失ったタンパク質は、進化の過程で生き残れませんね。

反対に、タンパク質機能に重要な影響を与えない部分のアミノ酸ならば、進化の途中でどんどん置換が起こり得るだろうし、またそのような機能に影響の無い場所のアミノ酸置換は、生き残っているタンパク質にも見られるはずでしょう。

つまり、配列の保存性が高い場所のアミノ酸置換は、タンパク質機能の変化という点でIntolerant（影響あり）。
配列保存性の低い場所のアミノ酸置換は、タンパク質機能変化にTolerant（影響なし）。

この方法は、正しいタンパク質ファミリー（オーソログタンパク）を集めてアライメントすることが重要です。その中で配列保存性を見ていくからです。　
UniprotやNCBIのnrなどから配列データを借りています。

さて、このように配列相同性のみを使用した予測ツールでは、タンパク質の立体構造までは見ていません。　
ループ構造部分に変異があるのか、疎水性が失われるのか、といった情報までは見ていないのです。
これらは別の予測ツール（Polyphenなど）が必要です。

NGSでSNVを見つけた。
Non-synonymousのSNVを絞り出した。
次に考えるのは、そのnsSNVの重要性です。
そのアミノ酸が変わったらタンパク質として生き残れないよ、という情報は、重要性の指標になります。
そのアミノ酸が変わっても、タンパク質は生き残っていたのなら、その置換はさほど重要ではないでしょう。
最初にSIFTなどで、配列保存性を基準にnsSNVの重要性をフィルターにかけ、残った重要そうなやつを、次にタンパク質立体構造上で確認していく方が近道のような気がします。

さて、実際にNGSでSNVを見つけてからこのツールを使って解析するにはどうしたら良いか？
先にURLを示した場所よりも、以下から入る方がお勧めです。
VCFフォーマットをSIFT専用のフォーマットにコンバートしてくれるからです。

http://sift.bii.a-star.edu.sg/www/SIFT_intersect_coding_submit.html

ここで、SAMtoolsなどで求めたVCFファイルをアップします。
そうするとしばらくして、コンバートされたファイルができますので、これをPCに保存します。
SNP用のファイルと、InDel用のファイルの2種類できます。
ここではSNP用ファイルを例に、落とします。

保存したら、upload here のリンクをクリックして、そのファイルをアップロードします。
その時、自分のE-mail アドレスと、オプションでどんなアノテーションを付けたいかを指定します。

結果は、HTML、またはテキストで参照できます。
これがSNPの結果

リンクできるところはそれぞれのDBにリンクしています。
HTMLだと、ちょっと扱いにくいので、私はテキスト（tsv）で落として、Excelで見ることもお勧めします。
例えば、こんな風に絞り込めば、新規っぽいSNPの機能を予測できます。

この絞り込みの条件は、

1. dbSNPに登録の無い "novel" SNPで
2. アミノ酸置換をもたらす Non-synonymous SNP

です。　ここで、SNP-Typeの右の、Predictionカラムでソートして、DAMAGING（=Intolerant）を選ぶのです。
その時、ScoreとMedian Infoにも注目します。
Scoreは0から1の値を取り、0.05以下のとき、タンパク質機能に影響あり（=Intolerant = Damaging）とされます。
Median Infoは保存性の信頼度で、0から4.32の値をとり、大きいほど信頼度が低くなります。　3.25以上のときはLow confidenceとして警告が出ます。

どうでしょう？　SNPの機能を調べるにはとても使いやすいツールだと思います。


 もっと詳しく知りたい方は、

 http://sift.jcvi.org/www/SIFT_help.html

Kumar P et al. Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm. Nature Protocols 4, 1073-82 (2009).

Ensembl - Variant Effect Predictor でSNPアノテーション

先日、VCFフォーマットを読みこんでSNPのアノテーションを付けるツールとして、SeattleSeq Annotationを紹介しました。ここ

ちょっと調べると、ほかにもいろいろありますね。
今日は定番？のEnsembl　です。
その中の、Variant Effect Predictor というウェブツールです。

Species　は、ここではヒト。　ゲノムバージョンと一緒に選びます。
ファイルの場所を指定して、フォーマット（VCF）を選びます。　
もう気づいたかもしれませんが、私はVCFを読みこめるツールが好きです。
VCFは変異レポートファイルのデファクトスタンダード。
書き込める自由度も高く、ほぼ満足です。
VCF自体を人間が読むことはあまり無いでしょうが。

追加するアノテーションを選びます。
ここからはヒトのデータのみです。

Non-Synonymous SNP（アミノ酸コーディング領域で起こる一塩基変異で、コードされるアミノ酸を置換するもの）について、その影響を予測するために、SIFT、PolyPhen、の2種類のスコアをアノテーションします。

その下のフィルタリングオプションは、マイナー・アレル・フリークエンシー（MAF）で、データをカットできます。
Nextを押して実行すると、検索が始まります。
こんな画面がでたら終了。
HTMLを「別のタブまたは別画面」で開いてみましょう。　そのまま開くと前に戻れないので注意！

HTMLのレポート画面です。　
Ensembl-IDとか、SNVの遺伝子上の場所、変異のアノテーション、アミノ酸変異があるものはその機能への影響度合い、などがリストで表示されます。

すぐ上のスクリーンショットには、Coding領域にあるNon-Synonymous SNVがありますね。

SIFTとPolyPhenの値も右にあります。

SIFT=tolerated(0.07); PolyPhen=benign(0.16) 

このうちPolyPhenの値とDescriptionについては以前SeattleSeqのところでちょっと書きましたね。

SIFTについてはまた次回。

さて、HTMLのほかに、テキストでダウンロードもできます。
これは後で編集するのに便利でしょう。

もうひとつは、
Ensemblのゲノムビューワーへのリンク。私は使ったことはことがないですが。　こんな感じでSNVの位置がわかります。

このEnsembl Variant Effect Predictorですが、残念ながらウェブ上で行うにはデータ量MAX７８０変異と、制限があるそうです。
それ以上の変異を解析するには、Perl APIという手があります。

制御領域はSNP検出に不向き？

遺伝子発現のオン・オフに重要な働きをする、CpG island と その多くが存在する5'-UTR
この部分はゲノムの他の部分よりGCコンテンツが高いことが知られています。
制御領域を次世代のシーケンサーで読んで、その変異を調べようという研究がありますが、それは注意が必要だ、というレポートを見つけました。

その名も、「Next generation sequencing has lower sequence coverage and poorer SNP-detection capability in the regulatory regions」　

このレポートでは、Cancer Genome AtlasからIlluminaとSOLiDのデータをダウンロードして、まずマッピングをたくさんのツールで行って比較しています。

比較しているのは、マッピングツールは

• Bowtie
• BWA
• SOAP2
• RMAP
• ZOOM
• Maq
• Novoalign
• SHRiMP

 SNP検出には

• MAQ
• SOAPsnp
• SNVmix

です。
マッピングツールのサマリーは、サプリメント資料に詳しくあります。
Smith-Waterman alignment系、Burrows-Wheeler transform系、Hash Table系、とそれぞれ。

マッピングプログラムの比較をまとめた論文は良く見かけるのですが、このレポートでは単にどれが速い・低メモリー・高カバレージだ、というのにとどまらず、タイトルにもありますが、

GC％が高い制御領域は、ゲノムの他の領域と比べてマッピングカバレージが極端に低い
よってCpG island などの制御領域でSNPを見つけたいときは気をつけよ！

と言っています。

私も今まで、カバレージと言っても、ゲノム平均カバレージ、あるいは遺伝子平均カバレージにのみ注意が行って、遺伝子の中の5'-UTR, Exon, Intron, 3'-UTR と分けて考えたことは無かったのでなるほどと思いました。

ちなみにこのレポートでは、

• BWT系のマッパー（Bowtie, BWA）がイルミナデータのマッピングには、ペアエンド、シングルリードのどちらでも最も成績が良かった。
• SOLiDデータのマッピングには、NovoalignCSが最も良かった。
• SNP検出はMAQが一番良かった。

と結論つけています。　ここは異論がある方もいるでしょうが。
ただ、なんとなくですが、無難な結論ではあります。

SeattleSeq Annotation : WebベースのHuman SNP機能予測システム

SeattleSeq　というサイトをご存じでしょうか？
ワシントン大学が運営する、ヒトSNV機能予測データベース＆システムです。

Web上でデータをやり取りします。
ユーザはSNVのリストファイルをアップロードして、欲しいアノテーションにチェックを入れ、メールアドレスを入力して、ボタンを押して、後は待つだけ。
インターネット上にデータを投げるわけですから、機密性が低い、公共データやデモデータ、自分で責任の負えるデータ、でやるのが無難です。

ファイルのフォーマットは割合と自由です。ここに詳しく
公共データのReadファイルから、マッピング、SNV検出してきた結果ファイルVCFフォーマットを入れてみましょうか。
VCFフォーマットの場合はちょっとコツがいりまして、SNVのみの結果にしなければいけません。
もしSNVとInDelが混ざって出力されているVCFの場合、InDelのデータは除きます。
VCFからSNVとInDelを分けて別々のファイルに保存するには、
こんなawkファイルを作って、

----------------Separate_SNV_and_InDel.awk     ここから------------------------------
/^#/    {
    print $0 > "snv_only.vcf";
    print $0 > "indel_only.vcf";   
    next;   
    }
   
/^[^\t]+\t[0-9]+\t[^\t]*\t[atgcATGC]\t[a-zA-Z]\t/   {
    print $0 > "snv_only.vcf";
    next;   
    }   
   
    {   
    print $0 > "indel_only.vcf";   
    next;   
    }
------------------------ここまで---------------------------------

awk -f separate_snp_indel.awk [あなたのファイル.vcf]

なんてコマンドを打つと、2つのファイルができて、snv_only.vcf ファイルにはSNVのみのデータが、
indel_only.vcf　にはInDelのみのデータが作られます。
このコマンドはあるBio-info仲間が教えてくれたものですが、簡単で重宝しています。

さて、もうひとつ、VCF 4.0 ファイルであることを示すため、できた snv_only.vcf ファイルの先頭行に
##fileformat=VCFv4.0
という一行を入れてあげます。

sed -e "1i ##fileformat=VCFv4.0" [あなたのファイル.vcf]  > [あなたのファイル2.vcf]
これで準備OK。

SeattleSeqの画面を見てみましょう。　http://snp.gs.washington.edu/SeattleSeqAnnotation131/
ヒトゲノムのバージョンが、hg18 と　hg19 の２つあるので、hg19を選択。

メールアドレスを入力し、ファイルをUpします。
あとは欲しいアノテーションにチェックを入れて、Submitボタンを押せばOK！
画面が切り替わり、ジョブの進行状況が確認できます。
しばらくすると、メールアドレス宛てに、メールで結果ファイルのURLが届きます。

そこをクリックして、ファイルをダウンロードします。

これがアップロードしたVCFファイル

そしてこれが、SeattleSeqの結果

New　で示したところが、SeattleSeqで新しくつけられたアノテーションですね。
VCFフォーマットを保ったまま、INFO　のところにrs番号をはじめ、SNVの分類、アミノ酸置換の種類、配列、GranthamScore、PolyPhen分類、などが新たについています。

以下の例で言うとこれらは、rs番号はrs3748597、SNVの分類(FD)はNon-synonymousでmissense、アミノ酸置換（AC）はイソロイシンからバリン（ILEからVAL）、GranthamScore（GS）は29、PolyPhen分類（PH）はpossibly-damaging、というようなことが書かれています。

1 888659 rs3748597 T C 149 . DP=11;AF1=1;CI95=1,1;DP4=0,0,3,8;MQ=60;
FQ=-60;DN;DT;DA=C/T;GM=NM_015658;GL=NOC2L;FG=missense;FD=missense;

AC=ILE/VAL;PP=300/750;GS=29;PH=possibly-damaging;CP=0.5990;CG=2.300;AA=C;
CN=2294,3274,30362;HA=12.1;HE=6.7;HC=6.2;DG;DV=by-frequency,by-cluster;
PS=MAAAGSRKRRLAELTVDEFLASGFDSESESESENSPQAETREAREAARSPDK
PGGSPSASRRKGRASEHKDQLSRLKDRDPEFYKFLQENDQSLLNFSDSDSSEEEE
GPFHSLPDVLEEASEEEDGAEEGEDGDRVPRGLKGKKNSVPVTVAMVERWKQA
AKQRLTPKLFHEVVQAFRAAVATTRGDQESAEANKFQVTDSAAFNALVTFCIRD
LIGCLQKLLFGKVAKDSSRMLQPSSSPLWGKLRVDIKAYLGSAIQLVSCLSETTV
LAAVLRHISVLVPCFLTFPKQCRMLLKRMVVVWSTGEESLRVLAFLVLSRVCRH
KKDTFLGPVLKQMYITYVRNCKFTSPGALPFISFMQWTLTELLALEPGVAYQHA
FLYIRQLAIHLRNAMTTRKKETYQSVYNWQYVHCLFLWCRVLSTAGPSEALQPL
VYPLAQVIIGCIKLIPTARFYPLRMHCIRALTLLSGSSGAFIPVLPFILEMFQQVDFN

RKPGRMSSKPINFSVILKLSNVNLQEKAYRDGLVEQLYDLTLEYLHSQAHCIGFP

ELVLPVVLQLKSFLRECKVANYCRQVQQLLGKVQENSAYICSRRQRVSFGVSEQ

QAVEAWEKLTREEGTPLTLYYSHWRKLRDREIQLEISGKERLEDLNFPEIKRRKM

ADRKDEDRKQFKDLFDLNSSEEDDTEGFSERGILRPLSTRHGVEDDEEDEEEGEE

DSSNSEDGDPDAEAGLAPGELQQLAQGPEDELEDLQLSEDD* GT:PL:GQ 1/1:182,33,0:63



ここで登場した、granthamScoreとpolyPhen、ほかにもscorePhastCons（CP）、consScoreGERP（CG）などは、塩基の変異によりアミノ酸が置換されるとき、これがタンパク質にどれくらい影響を及ぼすかの予測された指標です。

Grantham Scoreは、アミノ酸が変わるときに隣近所のアミノ酸との間にどれだけ化学的な影響を及ぼすか、分子量や極性などをもとに数値化しています。

conservative (0-50)、moderately conservative (51-100)、moderately radical (101-150)、radical (>151) というような分類もあるようです（Li et al. J Mol evol. 21, 58-71(1984).）

PolyPhen は、他の生物種とのホモログの配列情報から、変異があった箇所の配列保存性、化学的な特徴、タンパク質構造やドメイン情報をもとに、この変異がタンパク質の機能に及ぼす影響を数値化します。

probably damaging (>2.00)、possibly damaging (1.50-1.99)、potentially damaging (1.25-1.49)、borderline (1.00-1.24)、benign (0.00-0.99) という分類があるようで、SeattleSeqではこの分類表記になっていますね。

ちなみに、PolyPhenは現在、PolyPhen-2というのもありますが、こちらはアカデミック・非営利機関のみ利用できます。

SNVを探してから、missense のSNVのみを抽出し、そのアミノ酸変異がタンパク質にどれくらい影響を及ぼすのか、という研究は昔からありました。

現在はNGSのおかげでどんどん新しいSNVが見つかっています。

その塩基の変異が、タンパク質レベルでどのように働くのか、というテーマは今後も続いていくでしょう。

まずはPolyPhenなどの予測アルゴリズムを用いて、分子レベルの機能変化をスコア化し、数個に絞った後で、（今はここで終わっている論文がほとんどですが）今後はタンパク質の3次元構造モデリングなどを行ったり、アッセイに持ち込んで本当にタンパク質の活性が変化するのかを確かめたり、する研究も増えていくかもしれません。

Exome 解析 non-synonymous SNVを見つけた後は・・・

Agilent社のSureSelect、Illumina社のTruSeq、NimbleGen社のSeqCap といえば、最近盛んに宣伝されているExomeの実験用キットです。
Exomeというのは、ゲノム全体ではなく、その中のExon、実際はExon＋α　の配列のみを選択してシーケンスする解析方法です。
原理はまず、Exonの部分配列DNA（プローブ）が用意されていて、断片化されたゲノムと、溶液中またはアレイ上でハイブリし、Exonプローブにハイブリされた（キャプチャーされた）ゲノム断片だけを抽出します。
そのExon部分配列のみをシーケンスするというわけです。
実際Exonの上を何塩基ごとにプローブが設定されているのか、という情報は非公開ですが、Exon全てをキャプチャーできるようには設計されているそうです。

このキットが普及されるにつれ、Exome、特にヒトExomeの研究が盛んになってきました。
それまでは1000GenomeやCancer Genome Instituteのようなゲノムセンターでの研究が一般的でしたが、昨年末ごろから一般の大学の研究室レベルでも、ヒトExome実験が行われるようになってきたようです。

実験プロトコールもいろいろ大変なのでしょうが、私はデータ解析屋なので、データが出てからの料理に興味があります。

世の中に、フリーツールはたくさんありますが、Exome解析は、それらをいくつも組み合わせていく、一見面倒くさいけれどもエレガントなワークフローです。

ちょっと検索すると、
NGS Surfer's Wiki のリシークエンス
BioStar
などにパイプライン（ワークフロー）が出てきます。

まとめると、

1. リードをゲノムにマッピングし　（BWA）
2. 冗長性のあるリード・Duplicateを除去し　（SamtoolsやPicard）
3. キャプチャーしたExon領域だけを取り出し　（Bedtools）
4. SamtoolsでSNVを抽出し
5. SNVにアノテーションをつける non-synonymous SNV, missense, frame-shift など

というところまでを行っています。
SNVというのはSNPとほぼ同意語で、変異か多型かの違いです。厳密には異なりますが、ここでは同じ意味とします。

さて、Exomeの論文はたくさんありますが、その中でいいな、と思ったものを3つ

1. Vissers et al. A de novo paradigm for mental retardation. Nature Genetics 42, 1109-12 (2010).
2. Timmermann et al. Somatic mutation profiles of MSI and MSS colorectal cancer identified by whole exome next generation sequencing and bioinformatics analysis. PLoS One 22, e15661 (2010).
3. Ng et al. Exome sequencing identifies MLL2 mutations as a cause of Kabuki syndrome. Nature Genetics 42, 790-3 (2010).

 3のKabukiシンドロームは、10人の患者データセットがあるのですが、10人全部に共通する新規レアSNVを見つけようとしたらうまくいかなかったのですね。　そこで10人全員ではなく9人に共通、8人に共通、7人に共通というふうにレベルを下げていったのです。
そうするとNonsense置換またはフレームシフトを起こしていたMLL2遺伝子上の変異が、10人中7人に見つかったのです。
後はちゃんとCGHアレイとサンガーシーケンスで確認しています。

2のMSI/MSSは、直腸癌の種類（microsatellite instable / stable）が違う6人の患者サンプルで、MSIとMSSとで見つかった変異にどう違いがあるか、を見ています。
それぞれのサンプルで55,000程のSNVをリストした後、遺伝子の中にあるか、タンパクコード領域にあるか、dbSNP/1000Genomeに登録されていないか、カバレージは十分にあるか、Somaticな変異かどうか、とフィルタリングしていき、最終的に「タンパク質の機能を変化させる変異であるか」というところまで行き、数十程のSNVまで絞り込んでいます。
せっかくタンパク質の3D立体構造まで表示しているので、もう少しここから先が欲しかったです。

1、2、3、共通するのが de novo SNVを見つけた後にその変異が引き起こすであろうアミノ酸置換またはフレームシフトが、タンパク質の機能にどのような影響を及ぼすか、をスコア化していることです。

1はPhyloPとGrantham、2はPolyPhen（Polymorphism Phenotyping）とMutation Taster、3はGERP　というプログラムを使ってスコアを計算しています。
このような潜在的なnon-synonymous SNVの機能予測プログラムは他にも、SIFT（Sorting Intolerant from Tolerant）、PolyPhen‐2などがあります。
私も全部を知っているわけではないので（知っていなくてはいけないのでしょうが）、これらのうちの
 いくつか、多分PolyPhenとGranthamを、SeattleSeq というWebツールと一緒に紹介します。

発現解析パイプラインを作るぞ！ その2： TopHat の使い方 2

さて、無事TopHatが終了し、下のような画面が出て、
accepted_hits.bam, junctions.bed, insertions.bed / deletions.bed　というファイルができたとしましょう。

この3ファイルが何かについては前回説明しました。

TopHatは、コマンド実行時に指定したgtfファイルを使って、gtfに書かれている既知の転写配列情報から、Exon Junction　を認識してマッピングしてくれるのです。

Exon Junctionを認識するということは、Fusion Geneも特定してくれそうです。

それは、Tophat-fusion というツールがあるのですが、また今度別のときに書きます。

さて、うまくマッピングされ、BAMファイルができました。

では次に、Samtools を使って、BAMファイルを染色体番号順にソートします。

これは後で必要になります。

samtools sort ./tophat_SRP003186_Control/accepted_hits.bam ./tophat_SRP003186_Control/sort_accepted_hits

samtools sort ./tophat_SRP003186_MCF7/accepted_hits.bam ./tophat_SRP003186_MCF7/sort_accepted_hits

samtools sort ./tophat_SRP003186_SKBR3/accepted_hits.bam ./tophat_SRP003186_SKBR3/sort_accepted_hits
最後のアウトプットに.bam と書くと、出力が.bam.bam と2重拡張子になるので注意！

次に、後でUCSCとかIGVとかのGenome Viewerで見るために、ソートしたBAMファイルにインデックスをつけておきます。
samtools index ./tophat_SRP003186_Control/sort_accepted_hits.bam
samtools index ./tophat_SRP003186_MCF7/sort_accepted_hits.bam
samtools index ./tophat_SRP003186_SKBR3/sort_accepted_hits.bam
拡張子.bai というファイルができるはずです。　これがインデックスです。

ついでにソートしたBAMファイルを、SAMフォーマットにしておきましょうか？　別に必要無いですけど、samtools に慣れていないひとは練習がてら。　

ちなみにIGVは、SAMファイルも取り込めます。その際には、IGV-toolsでインデックスをつけるひつようがあります。

また、以前紹介した、SAM-Mate　というWindows 対応のRNA-SEQ解析ツールは、SAMでないと受け付けてくれません、今のところ。　

samtools view ./tophat_SRP003186_Control/sort_accepted_hits.bam -o ./tophat_SRP003186_Control/Control_sort_accepted_hits.sam
samtools view ./tophat_SRP003186_MCF7/sort_accepted_hits.bam -o ./tophat_SRP003186_MCF7/MCF7_sort_accepted_hits.sam
samtools view ./tophat_SRP003186_SKBR3/sort_accepted_hits.bam -o ./tophat_SRP003186_SKBR3/SKBR3_sort_accepted_hits.sam

 

発現変動を見るには、Cufflinksのコマンドのひとつ、Cuffdiff　というものを実行します。

ソートしたBAMファイル、既知トランスクリプトのGTFファイルの2つを使います。

cuffdiff -p 4 -v /Path_to_gtf/refGene.gtf ./tophat_SRP003186_Control/sort_accepted_hits.bam ./tophat_SRP003186_MCF7/sort_accepted_hits.bam -o ./cuffdiff/MCF7

cuffdiff -p 4 -v /Path_to_gtf/refGene.gtf ./tophat_SRP003186_Control/sort_accepted_hits.bam ./tophat_SRP003186_SKBR3/sort_accepted_hits.bam -o ./cuffdiff/SKBR3

コントロール v.s. MCF-7、コントロール v.s. SK BR-3 の2つの比較ですね。

本当はリプリケートがあると良かったのですが、これはN=1です。

パラメータは単純にデフォルトで行っています。

唯一、CUPを　-p 4 にしたことくらいでしょうか。

Cuffdiffのアウトプットファイルは地味で、タブ区切りテキストファイルができるだけです。

もちろん加工すればGenome Viewer上でも表示できますが。

FPKM tracking ファイルと、Differential expression 結果ファイルの2種類があります。

1) FPKM tracking ファイル：　トランスクリプトごとのFPKMを出力し、各cuffdiffランごとに4ファイルできる
isoforms.fpkm_tracking：　Transcript FPKMs
genes.fpkm_tracking：　Gene FPKMs.  同じgene_idを持つトランスクリプトのFPKMの積算
（cds.fpkm_tracking：　Coding sequence FPKMs.）
（tss_groups.fpkm_tracking：　Primary transcript FPKMs）

2) Differential expression 結果ファイル：　　各cuffdiffランごとに4ファイル

isoform_exp.diff：　Transcript differential FPKM.
gene_exp.diff：　Gene differential FPKM.
（tss_group_exp.diff：　Primary transcript differential FPKM）
（cds_exp.diff：　Coding sequence differential FPKM）

このうち、発現変動の結果ファイル、isoform_exp.diff　をExcelで開いてみるとこんな感じです。

トランスクリプトのIDや名前、ポジション情報が先ずあり、そのあとHカラムくらいから、実際のカウンティングデータがきます。

• FPKMx：　FPKM of the gene in sample x
• FPKMy：　FPKM of the gene in sample y
• ln(FPKMy/FPKMx)：The natural log (自然対数) of the fold change y/x
• test stat：　FPKMの変動が有意であるかどうかの統計量
• p value：　補正無しの p-value
• q value：　FDR-adjusted p-value
• significant："yes" or "no"　マルチサンプル比較の場合はBenjamini-Hochberg補正後、FDR-adjusted p-valueが、cuffdiffランのパラメータ --FDR （デフォルト：0.05）より小さい場合は、"Yes"　有意だとされる

Ratioは、自然対数のログ表記という点にご注意ください。

また、Excelはsignificant = Yes で絞り込んでいます。

こんな風にして、RNA-Seq　の実験データから、発現変動の有意な遺伝子を解析することができるのです。　

実際は、パラメータの調整などが必要ですが、感じはわかって頂けたと思います。

このあと考えられる操作は、

1. 発現変動の有意だった遺伝子のポジション情報からBedファイルを作り、
2. 先ほどソート済みのBAMファイルをGenome Viewerと一緒にとりこみ、
3. ゲノム上で発現の変動と、マッピングの様子を同時に見る

ことでしょうか。　これはまた、Viewerに関する書き込みのときに詳しく説明します。

つづく

SNP検出のアルゴリズム2

SNP検出で前回、NQSを書いた。
クオリティスコアについては、様々な閾値で設定することができるが、文献で見たことがあるのは、
“11-base NQS 20/15 threshold”
というもの。
これはSNP検査対象の塩基を中心に前後5ベース、の計11ベースのクオリティを見て、
中心は20、前後5ベースは15、以上あるとき、の中心塩基を検査する、というもの。
先のCLC-Bio Genomics Workbench を例にしたときのデフォルト。
NQSは次世代シーケンサーの前から良く使われていたらしい。

Shen Y, Wan Z, Coarfa C, Drabek R, Chen L, Ostrowski EA, Liu Y, Weinstock GM,
Wheeler DA, Gibbs RA, Yu F.
A SNP discovery method to assess variant allele probability from next-generation resequencing data. Genome Res. 2010 Feb;20(2):273-80. Epub 2009 Dec 17. PubMed PMID: 20019143

僕もこの閾値が絶対だとは思わないので、いろいろ変える必要はありそうだ。

実はSNP検出に関して、この閾値よりももっと大事なのは、マッピングするときのアルゴリズム。
つまりギャップありマッピングか、ギャップなしマッピングかで、結果は大きく変わってくるのだ。
ギャップありの場合、1塩基単位のInsertion、Deletionも検出可能になる。
ギャップありマッピングは、BWAなどで可能だが、精度が高い分ランに時間がかかる。
しかし小さいInDelを検出できるとともに、SNPの擬陽性確率が低いという報告もある。
昨年末のBMBでもあるポスター発表者と話したのだが、BWAなどのギャップありマッピングは、SNP検出時に行うマッピングとしては、Bowtieなどのギャップ無しマッピングと比べて優れているそうだ。

ちなみに有償ソフトのNextGENeも、ユニークなSNP検出アルゴリズムを用いている。
計算式が複雑すぎて僕は理解できないが、興味のあるひとはのぞいてみると良い。
http://www.softgenetics.com/NextGene_UsersManual_version_2-0.pdf
114ページ、Overall Mutation Score の章から先

SNPの検出アルゴリズム

Single Nucleotide Polymorphism（SNP：1塩基多型）は、遺伝学では疾患マーカーとして良く研究される分野である。　たった1塩基の差が、その病気を決定づけたり、ある薬の効果を全く無くしたり、あるいは重篤な副作用をもたらしたりするからだ。

実際、その塩基の差が、病気や副作用に関係するかどうかは、遺伝統計学の力でなんとか証明するのだけれど、それはここではひとまず置いておいて、その「塩基の違い」はどうやって見つけるのかについて考えてみよう。　もちろんシーケンスからだ。
そう、始まりはいつもマッピングファイル。

SNPの検出は、マッピングされた後のデータから始める。
有名、というか一般的な方法に、NQSという手法がある。これはNeighborhood Quality Standardの略で、リファレンス配列にマッピングしたときのリードのクオリティスコアをもとにSNPを判断する方法だ。

参考になるのはこの文献
Brockman et al. Quality scores and SNP detection in sequencing-by-synthesis systems. (2008) Genome Res. 18, 763-70. (PMID: 18212088)
ここで提唱されている手法は、多くの論文でも採用されているし、いくつかの有償ソフトでも基本的なところではそのままデフォルトになっている。　ここから先は、そんな有償ソフトの一つ、CLC-Bio社のGenomics Workbenchを例に説明しよう。

先ず、リファレンスにマッピングされている複数のリードを想像して欲しい。　そのリード配列のあるところの塩基がAで、同じ場所のリファレンス配列の塩基もA、ところが複数あるリードのなかには、その場所がCのものもある。さて、この塩基が、SNPであるかどうか？

一番上はリファレンスで、その下にたくさんあるのがマップされたリードの配列だ。
真ん中の縦の列に注目して欲しい。このAとCがSNPであるかどうかを判断するために、リードの各真ん中の塩基の、右と左の塩基のクオリティスコアを見る。

ここから先、クオリティスコアなどの閾値の例を示す。
上の絵の横のハイライトを見て欲しい。
CATAT A AATAA と、真ん中のAを中心に右に5つ、左に5つの塩基がハイライトされている。
この左右10個の塩基のクオリティスコアの平均が、最低15以上でなければならない。
また、中心の全ての塩基のクオリティスコアは、最低20でなければいけない。
これら5＋5＋1、計11塩基のことを、Windowサイズ11と呼ぶ。

さらに条件を加えるならば、同Windowの中にはミスマッチ・ギャップが2つ未満でなければいけない。
これらの条件から外れた場合、そのリードは、この箇所の塩基において（真ん中の縦の列）、SNP検出の候補にはならない。
また、縦の列、中心の塩基すべてのクオリティスコアが最低20に満たない場合、この塩基はSNP判定の対象にならない。

ここで、SNP判定の対象になる塩基と、その塩基の判定に用いられるリードが残るわけだ。
SNP判定には、
1．何本のリードがカバーしているか
2．その変異は、全体の何%必要か、あるいは数はいくつ必要か
などを基準に決められる。

このようにして検出されたSNPは、リストになって出力され、どのSNPがアミノ酸置換を起こしたか、などが参照できる。

フリーのソフトは、僕はまだ試していないが、やっているひとはいる。
先の文献、一度PubMedで調べて、その文献を引用している論文を見てみると、フリーツールが見つかるのではないか。
NQSは結構有名なのだ。