# 茨城大学集中講義2019

## 第1回 NGSのシーケンサー原理紹介

* [[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/1.pptx|講義1]]

## 第2回 RNA-seqデータ解析

最初に手元のコンピュータだけでも解析できる[[https://github.com/c2997108/OpenPortablePipeline/blob/master/README_jp.md|Portable Pipeline]]という私が開発中のソフトウェアを使って解析を行い、その後Maserも使ってみる。事前準備として、[[バイオインフォマティクス_ツールインストール方法_2019|このリンク先のページ]]に従ってWindowsならばWSL、MacであればDockerのインストールを済ませておく。

今回使用するデータ

* [[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/osanai-rna-seq-snp-100k/fastq.zip|シーケンスデータ]] [[https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4806880/|Osanai-Futahashi et al. Heredity (Edinb). 2016 Feb; 116(2): 135–145.]]より

* [[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/osanai-rna-seq-snp-100k/Bomo_genome_assembly.minus.scaf034.scaf395.plus.NC_002355.fasta.gz|カイコゲノム配列2016改変版]]

* [[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/osanai-rna-seq-snp-100k/Bomo_gene_models.withnote.plus.NC_002355.gff3.gz|カイコ遺伝子モデル2017 注釈付き]]

* ファイル名がDRRXXXXXとなっていて、サンプル情報がわからないので、[[https://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=DRA003068|DDBJ]]のSequence Read Archive (SRA)から入手する。SRAデータベースの構造は[[https://www.ddbj.nig.ac.jp/dra/submission.html|こちら]]を参照。

* [[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/osanai-rna-seq-snp-100k/dmel-all-translation-r6.29.fasta.gz|遺伝子アノテーション用にショウジョウバエの遺伝子のアミノ酸配列(FASTA)]]

* [[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/hisat-index/hisat2_index.zip|カイコゲノムのHISAT2のインデックスファイル(メモリが15GB以上あるPCなら自分で作れるから不要)]]

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では、次のスクリプトを順番に実行する。

 1. mapping-illumina~HISAT2

 (今回使用するデータはstrand specificなデータではないため、hisat2 optionの--rna-strandness RFを削除する。)

 2. RNA-seq~featureCounts

 (入力ファイルとして、1.の中の```output.unsorted_bam```の中のBAMファイルすべてと、```input_3/Bomo_gene_models.withnote.plus.NC_002355.withname.gtf```を指定する。また、strand specificではないため、featureCounts optionの-s 2を削除する。)

 3. statistics~sample_distance

 (入力ファイルとして、```count.id.forDEGanalysis.txt```を指定する。)

 4. statistics~all-sample-combinations-DESeq2-edgeR

 (入力ファイルとして、```count.id.forDEGanalysis.txt```、サンプルの状態を記述したファイル(自分で作る)を指定する。input3は空でOK)

 5. convert~extract-FASTA-from-GFF3orGTF-and-FASTA

 6. annotation~blast

 (入力ファイルとして、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には未登録。


## 第3回 SNP解析

SNP解析は基本的にはGATKのベストプラクティスワークフローを参照するのが無難。10万人規模のヒトゲノム解析にはこのツールが使われている。

https://software.broadinstitute.org/gatk/best-practices/

データの種類によって異なる3つの解析ツールがある。

*  Germline SNPs + Indels

 遺伝性疾患の原因変異を探索する際に使用される最もメジャーな使用方法であり、BWAでマッピングして、picardで重複除去を行い、GATKでSNPをコールする。Portable Pipelineで使用する下記のスクリプトではRAD-seq用に重複除去を行わない方針を採用している。

 WGS~genotyping-by-GATK

* Somatic SNVs + Indels

 主にがん細胞を対象として、変異を見つける場合に使う。BWAでマッピングして、Mutect2でSNPをコールする。ゲノムがduploidではない状態、極端に言えば変異が1%しかなくても変異を検出可能になるように設計されている。Portable Pipelineでは未登録。

* RNAseq SNPs + Indels

 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に変更する)

## 第4回 Excelファイル後のデータ解析

統計的な解析はやはりRを使ったほうが便利だけど、そういった解析が終わった後でデータを俯瞰するにはやっぱりExcelが便利。Excelで良く使う機能としては、フィルター機能、条件付き書式、グラフ作成あたりか。TXTファイルをExcelで開くときの注意点などは[[データ解析入門2018#テキストファイルをexcelで読み込む|こちら]]。

* Gene Ontologyやパスウェイ解析（エンリッチメント解析）

 これらの解析を大雑把に説明すると、DEG解析で出てきた遺伝子が幾つかあって、それらの遺伝子に偏っている特徴や機能は何か大雑把に要約したいときに使う。目的が大雑把なので、あまり細かいことを気にしてはいないツールが多いし、気にしても仕方がない（解がない）ほうが普通である。

 ツール側の内部で行われている手順としては、遺伝子名とGene Ontology、パスウェイなどを結び付けるデータベースを使って遺伝子名をそれらの情報に変換し、対象の生物種が持っている全遺伝子のGene Ontology、パスウェイと比較して与えられた遺伝子群に有意に多いカテゴリを抜き出してくる。

 遺伝子名を投げるだけという簡単な操作で結果を返してくれるサイトが多く、偶に発現量の情報も加味してくれるサイトもある。ただし、これらのサイトの多くはヒトやマウスが対象であり、ショウジョウバエやゼブラフィッシュなどがギリギリ使えるかどうかといった状態で、カイコなどはまず対象外となることが多い。ショウジョウバエであれば、下記のサイトが便利かも。

 https://amp.pharm.mssm.edu/FlyEnrichr/

* エンリッチメント解析を行うためのExcelでの下準備

 [[http://www.suikou.fs.a.u-tokyo.ac.jp/yosh/lib/exe/fetch.php?media=result.edger.isoforms.count_table.c108.p50t.txt.c108.down.p50t.up.txt.blastx.xlsx|前日までの解析でExcelファイルを作れなかった人向けのテストデータ]]

 データ解析を行うときにプログラミングスキルは必須ではなく、目的の処理を行うツール名を知っているか、ツールのインストールが出来るか、オプションを間違えずに設定できるか、必要なデータベースにアクセスできるか、ウェットの実験と同じく目的のデータだけ使わずにコントロールとして別のデータで検証する姿勢があるかといったことのほうが圧倒的に大事だと感じているが、そうは言ってもちょっとした一回限りのフォーマット変換や、集計などは頻繁に行うため、簡単なプログラミングは知っておいたほうが便利である。

 今回はエンリッチメント解析に投げる遺伝子名のリストをExcelを使って抽出してみる。

 今回Excelで使う関数は、FIND, LEFT, RIGHT, LEN関数である。

* Excelで情報の集約

 複数のファイルに別の情報があって、それを一つのファイルにまとめることが出来るとフィルター機能などで一度に選別出来て便利である。例として、二群間比較結果のファイルと遺伝子の詳細情報ファイルの情報を一つにまとめてみる。

 使用する関数は、VLOOKUPである。


## 第5回 IGV、Geneiousを使用したデータ解析 + Maserの利用方法紹介

### 私が良く使う可視化ツールの紹介。

* [[https://software.broadinstitute.org/software/igv/|Integrative Genomics Viewer (IGV)]]

 マッピングした結果を見るときに使用する。NGSでどんな解析を行ったとしても、一度は使うことになるはず。マッピングの様子を見ると、バイアスなく遺伝子領域をシーケンス出来ているかとか、参照ゲノムが間違っていそうかとか、ヘテロ・ホモ変異が想定通りかとか、色々な情報を読み取ることが出来る。RNA-seqであれば、DEG解析後に興味を持った遺伝子領域付近を、一度はIGVで開いて見るべき。

 使用方法は、IGVを起動した後、まずはゲノムファイル(FASTA形式)を上部メニューの中の「Genomes」->「Load Genome from File」から開く必要がある。それから遺伝子モデルのGFF3ファイルやマッピング結果のBAMファイルを上部メニューの中の「File」->「Load from File」から開いていく。

* [[https://www.geneious.com/|Geneious]]

 マルチプルアライメントを表示するときに便利。アンプリコンメタゲノム解析でOTUごとの違いを可視化したいときに使える。系統樹も簡単に作れるけど、プロの人たち(?)はモデル選択をバリバリ使ったりするので、そういう人はアライメントの後は[[https://www.megasoftware.net/|MEGA]]を使うみたい。

### Maserの利用方法

あらかじめアカウントを取得していると思うので、https://cell-innovation.nig.ac.jp/members/maser3/ にログインするところから説明する。基本的な流れは、プロジェクトを作成し、プロジェクトページにデータ(FASTQやFASTA等)をアップロードし、

1. 上部メニューの「Project」を選択し、「Create New Project」で新しいプロジェクトを作成する。この時、「Name」の項目は記入が必要。

2. 作成したプロジェクトを開いたら、ページ中央の「Upload My Data」からデータをアップロードする。今回はWEBブラウザだけで完結する「Upload via HTTPS」を使った方法を紹介するが、たくさんのデータをアップロードする場合は、SFTPを使ったほうが便利。

 「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」ボタンを押す。

3. Uploadが終わったらプロジェクトページをリロードして更新すると、アップロードしたデータのアイコンが見えるようになる。必要なデータセットを全てアップロードしたら、ひとまず「multi fastq (paired-end)」を選択して、「f(x)」と書かれたボタンを押してみる。そうすると、「multi fastq (paired-end)」のデータセットタイプで使用可能な解析メニューが表示される。(恐らく現在は「HISAT2->StringTie」と「Trinity->kallisto」パイプラインのみ表示される。)

4. 解析したいパイプラインを選んで「Analysis」ボタンを押し、不足しているデータセットがあれば適宜選択し、「Set option and run」をクリックしてオプションを確認し、必要があれば適宜変更してから「Run」ボタンを押す。

5. ページ上部の「Analysis Status」タブをクリックすると、ジョブの実行状態が確認できる。ジョブが終了していれば、プロジェクトページに結果が表示されるので、解析結果のデータセットのアイコンをクリックすると表示されるダウンロードボタン「↓」をクリックすると結果が表示される。

## 第6回 NGSを使うようになって分かってきた研究成果

* [[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/6.pptx|講義6]]

## 第7回 メタゲノム解析

メタゲノム解析は16S rDNAや特定の遺伝子領域のみシーケンスするアンプリコンメタゲノム解析と、採取したDNAをそのまま全部読むショットガンメタゲノム解析に分かれる。

* アンプリコンメタゲノム解析

 基本的にはリードをBlastでデータベースに当てていき、データベースに登録されている分類名を割り当てていく。単純なために色々な解析ツールや受託解析業者が存在する。Portable Pipelineで使用するスクリプトは

 metagenome~silva_SSU+LSU (16S rDNA用)

 metagenome~use-genbank-fasta-as-reference (菌類のITSといった配列をNCBIからダウンロードしてくる場合)

 また、WindowsやMacでほかにも手軽に解析する方法がある。[[2018|ナノポアで読んだ食品のメタゲノムデータからの食品当てクイズ]]参照。

 使用するデータは、[[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/foodmeta/foodmeta.zip|こちら]]

* ショットガンメタゲノム解析

 ショットガンメタゲノム解析は得られる情報量は多いが、1TB程度のメモリのサーバが必要になったりする等解析コストが非常に高く、また解析フローが共通化されてお決まりのパターンがあるとは言い難い。解析フローの一例を挙げると、MEGAHITでアセンブルを行い、MetaBATでメタゲノムの中から個別のバクテリアゲノムを分離し、Metaxaでバクテリアゲノムのアノテーション(種同定)を行う。個人的な感想としては、20年後くらいに環境の違いを議論するときのために現在はデータを取っておくフェーズなのかなという印象。既知遺伝子の新しいバリエーションを探すには有効だと思われる（例：CRISPR/Cas探索など）。

データ解析が終わっていない人は、[[http://www.suikou.fs.a.u-tokyo.ac.jp/yosh/lib/exe/fetch.php?media=result.zip]]からダウンロードしてください。

## 第8回 NGSを使った研究実例紹介

* [[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/8-1.pptx|講義8-1]]

* [[http://suikou.fs.a.u-tokyo.ac.jp/ibaragi/8-2.pptx|講義8-2]]