最初に手元のコンピュータだけでも解析できるPortable Pipelineという私が開発中のソフトウェアを使って解析を行い、その後Maserも使ってみる。事前準備として、このリンク先のページに従ってWindowsならばWSL、MacであればDockerのインストールを済ませておく。
今回使用するデータ
RNA-seqの解析の流れは、参照ゲノム・遺伝子モデルがあるかどうかで大きく手法が分かれる。
参照ゲノム・遺伝子モデルがある場合でも、さらに既存の遺伝子モデルのみを使うのか、新規の遺伝子モデルを使うのかでも手法が異なる。
Trinityでトランスクリプトームをde novoアセンブルし、kallistoで発現量算出、DESeq2, edgeR, sleuthでDEG解析を行う。Trinityの公式HPに下流の解析方法なども例があり参考になる。https://github.com/trinityrnaseq/trinityrnaseq/wiki/Post-Transcriptome-Assembly-Downstream-Analyses
Portable Pipelineで使用するスクリプトは、
RNA-seq~Trinity-kallisto-sleuth
(Maserでも使用可能)
HISAT2で参照ゲノムにRNA-seqデータをスプライス部位を考慮してマッピングし、StringTieで遺伝子モデルを構築し、DESeq2, edgeR, ballgown, cuffdiffでDEG解析を行う。少し前だとTopHat2でマッピングして、Cufflinksで遺伝子モデルを構築し、cuffdiffでDEG解析を行うのがスタンダードだった。https://www.nature.com/articles/nprot.2012.016
確かにHISAT2→StringTieまでであれば、TopHat2→Cufflinksよりも速くて正確そうではあるが、最後のDEG解析まで含めると正直どっちが良いのかなと思う。
Portable Pipelineで使用するスクリプトは、
RNA-seq~HISAT2-StringTie-DEGanalysis
(Maserでも使用可能)
HISAT2で参照ゲノムにRNA-seqデータをスプライス部位を考慮してマッピングし、featureCountsで既知遺伝子の発現量を算出し、DESeq2, edgeRでDEG解析を行う。シンプルで速い解析手法。シングルセル解析等の3'や5'の一部しか読まないような解析では、基本的にこのパターンの解析を行う。マッピングにはHISAT2やSTARが良く用いられ、発現量算出にはfeatureCountsやStringTie、その他諸々のカウントツールが存在する。今回はこの手法で解析を行ってみる。
Portable Pipelineでは、次のスクリプトを順番に実行する。
(今回使用するデータはstrand specificなデータではないため、hisat2 optionの–rna-strandness RFを削除する。)
(入力ファイルとして、1.の中のoutput.unsorted_bamの中のBAMファイルすべてと、input_3/Bomo_gene_models.withnote.plus.NC_002355.withname.gtfを指定する。また、strand specificではないため、featureCounts optionの-s 2を削除する。)
(入力ファイルとして、count.id.forDEGanalysis.txtを指定する。)
(入力ファイルとして、count.id.forDEGanalysis.txt、サンプルの状態を記述したファイル(自分で作る)を指定する。input3は空でOK)
(入力ファイルとして、queryには5.の出力ファイルのBomo_gene_models.withnote.plus.NC_002355.gff3.with-geneid.genes.fastaを、databaseにはショウジョウバエの遺伝子のタンパク質配列を、Gene expression fileとして例えば4.の出力ファイルのresult.edgeR.isoforms.count_table.C108.p50T.txt.C108.down.p50T.up.txtを指定する。DNA→タンパク質での相同検索なので、blastのprogramとしてblastxを指定する。)
遺伝子、もしくはスプライスバリアントの配列ファイルだけがあってゲノム配列を使わない場合が最もシンプルであり、kallisto等のツールで発現量を算出し、sleuth等のツールでDEG解析を行えばよい。最も単純だけど、同一遺伝子内のスプライスバリアント間の関係を見たくなったりして、ゲノムにマッピングして確かめたくなることが多い。Portable Pipelineには未登録。
SNP解析は基本的にはGATKのベストプラクティスワークフローを参照するのが無難。10万人規模のヒトゲノム解析にはこのツールが使われている。
https://software.broadinstitute.org/gatk/best-practices/
データの種類によって異なる3つの解析ツールがある。
遺伝性疾患の原因変異を探索する際に使用される最もメジャーな使用方法であり、BWAでマッピングして、picardで重複除去を行い、GATKでSNPをコールする。Portable Pipelineで使用する下記のスクリプトではRAD-seq用に重複除去を行わない方針を採用している。
WGS~genotyping-by-GATK
主にがん細胞を対象として、変異を見つける場合に使う。BWAでマッピングして、Mutect2でSNPをコールする。ゲノムがduploidではない状態、極端に言えば変異が1%しかなくても変異を検出可能になるように設計されている。Portable Pipelineでは未登録。
RNAseqデータからSNPを探索する場合に推奨されているパイプラインとしては、STARでマッピングして、picardで重複除去を行い、GATKでSNPをコールする。Portable Pipelineで使用するスクリプトは
RNA-seq~SNPcall-STAR-GATK
今回はこのRNAseqからの解析を行ってみたいところだが、メモリが20GB以上必要そうなので、簡略化して次のように行う。
1．mapping-illumina~HISAT2 (第2回RNA-seq解析の結果をそのまま使うので実行する必要はない）
2．SNPcall~bcftools-mpileup (時間がかかるので１番染色体だけのSNPを解析するため-r Bomo_Chr1をbcftools mpileup optionに追加する)
3．graph~VCFtoTree (データ数が少ないので、閾値を下げて数を増やすため、GQ thresholdを0に変更する)
統計的な解析はやはりRを使ったほうが便利だけど、そういった解析が終わった後でデータを俯瞰するにはやっぱりExcelが便利。Excelで良く使う機能としては、フィルター機能、条件付き書式、グラフ作成あたりか。TXTファイルをExcelで開くときの注意点などはこちら。
これらの解析を大雑把に説明すると、DEG解析で出てきた遺伝子が幾つかあって、それらの遺伝子に偏っている特徴や機能は何か大雑把に要約したいときに使う。目的が大雑把なので、あまり細かいことを気にしてはいないツールが多いし、気にしても仕方がない（解がない）ほうが普通である。
ツール側の内部で行われている手順としては、遺伝子名とGene Ontology、パスウェイなどを結び付けるデータベースを使って遺伝子名をそれらの情報に変換し、対象の生物種が持っている全遺伝子のGene Ontology、パスウェイと比較して与えられた遺伝子群に有意に多いカテゴリを抜き出してくる。
遺伝子名を投げるだけという簡単な操作で結果を返してくれるサイトが多く、偶に発現量の情報も加味してくれるサイトもある。ただし、これらのサイトの多くはヒトやマウスが対象であり、ショウジョウバエやゼブラフィッシュなどがギリギリ使えるかどうかといった状態で、カイコなどはまず対象外となることが多い。ショウジョウバエであれば、下記のサイトが便利かも。
前日までの解析でExcelファイルを作れなかった人向けのテストデータ
データ解析を行うときにプログラミングスキルは必須ではなく、目的の処理を行うツール名を知っているか、ツールのインストールが出来るか、オプションを間違えずに設定できるか、必要なデータベースにアクセスできるか、ウェットの実験と同じく目的のデータだけ使わずにコントロールとして別のデータで検証する姿勢があるかといったことのほうが圧倒的に大事だと感じているが、そうは言ってもちょっとした一回限りのフォーマット変換や、集計などは頻繁に行うため、簡単なプログラミングは知っておいたほうが便利である。
今回はエンリッチメント解析に投げる遺伝子名のリストをExcelを使って抽出してみる。
今回Excelで使う関数は、FIND, LEFT, RIGHT, LEN関数である。
複数のファイルに別の情報があって、それを一つのファイルにまとめることが出来るとフィルター機能などで一度に選別出来て便利である。例として、二群間比較結果のファイルと遺伝子の詳細情報ファイルの情報を一つにまとめてみる。
使用する関数は、VLOOKUPである。
マッピングした結果を見るときに使用する。NGSでどんな解析を行ったとしても、一度は使うことになるはず。マッピングの様子を見ると、バイアスなく遺伝子領域をシーケンス出来ているかとか、参照ゲノムが間違っていそうかとか、ヘテロ・ホモ変異が想定通りかとか、色々な情報を読み取ることが出来る。RNA-seqであれば、DEG解析後に興味を持った遺伝子領域付近を、一度はIGVで開いて見るべき。
使用方法は、IGVを起動した後、まずはゲノムファイル(FASTA形式)を上部メニューの中の「Genomes」→「Load Genome from File」から開く必要がある。それから遺伝子モデルのGFF3ファイルやマッピング結果のBAMファイルを上部メニューの中の「File」→「Load from File」から開いていく。
あらかじめアカウントを取得していると思うので、https://cell-innovation.nig.ac.jp/members/maser3/ にログインするところから説明する。基本的な流れは、プロジェクトを作成し、プロジェクトページにデータ(FASTQやFASTA等)をアップロードし、
「Data Label」にデータを識別するための名前を入力する。
「Data Type」にデータの種類を指定するが、注意事項として、Portable Pipelineにも登録されている比較的最近のデータ解析パイプライン用にデータをアップロードする場合、「Data Type」にはかなり下のほうまでスクロールして「multi fastq (paired-end) : This dataset type can contain 1 - 1000 fastq files.」を指定する。これは一つのデータセットに複数のFASTQファイルを含むデータセット。古めのパイプライン(数は多い)を利用する場合は、「fastq (paired-end)」や「fastq (single-end)」を選ぶ。こちらは1つのデータセットに1組のペアエンドのFASTQか、1つのシングルエンドのFASTQのみアップロード可能。また、FASTA形式は、塩基配列なのか、アミノ酸配列なのかでデータセットが異なり、DNAの場合は「fasta (nucleotide)」を選択する。サンプル情報(サンプル毎の実験条件が書かれたファイルで統計解析の際に使用される)を記入したテキストファイルは「SampleList」形式を選択する。目的の形式を探すときは、Ctrl-Fを押してブラウザのページ内検索機能を使うと便利。
「Data File:」には適宜「Add file」ボタンを押しながら、「ファイルを選択」ボタンから必要なファイルをすべて選択する。選択し終えたら「Upload」ボタンを押す。
メタゲノム解析は16S rDNAや特定の遺伝子領域のみシーケンスするアンプリコンメタゲノム解析と、採取したDNAをそのまま全部読むショットガンメタゲノム解析に分かれる。
基本的にはリードをBlastでデータベースに当てていき、データベースに登録されている分類名を割り当てていく。単純なために色々な解析ツールや受託解析業者が存在する。Portable Pipelineで使用するスクリプトは
metagenome~silva_SSU+LSU (16S rDNA用)
metagenome~use-genbank-fasta-as-reference (菌類のITSといった配列をNCBIからダウンロードしてくる場合)
また、WindowsやMacでほかにも手軽に解析する方法がある。ナノポアで読んだ食品のメタゲノムデータからの食品当てクイズ参照。
使用するデータは、こちら
ショットガンメタゲノム解析は得られる情報量は多いが、1TB程度のメモリのサーバが必要になったりする等解析コストが非常に高く、また解析フローが共通化されてお決まりのパターンがあるとは言い難い。解析フローの一例を挙げると、MEGAHITでアセンブルを行い、MetaBATでメタゲノムの中から個別のバクテリアゲノムを分離し、Metaxaでバクテリアゲノムのアノテーション(種同定)を行う。個人的な感想としては、20年後くらいに環境の違いを議論するときのために現在はデータを取っておくフェーズなのかなという印象。既知遺伝子の新しいバリエーションを探すには有効だと思われる（例：CRISPR/Cas探索など）。
データ解析が終わっていない人は、http://www.suikou.fs.a.u-tokyo.ac.jp/yosh/lib/exe/fetch.php?media=result.zipからダウンロードしてください。