バイオインフォマティクスの基礎
CLCなどの有償のソフトウェアを使用しないで，自分でMacや大型計算機（スーパーコンピューター：スパコン）を使用してNGSのデータ解析を行うには，UNIXが使えることが必須です。はじめに，知っておくと便利なコマンドやよく使うコマンドをご紹介します。また，大型計算機（スパコン）で使用するコマンドなどをご紹介致します。
1-1) PATHの概念
UNIXコマンドを覚える前にPATHについて覚える必要があります。ファイルやデータはどこにあって，どこに移動したいのかなど，すべての操作に対してPATH(ファイルの場所)を指定する必要があります。
このページの先頭へ
1-2) UNIXコマンドを覚えよう
みなさんが，いつもPC上でデータをコピーしたり，新しいディレクトリ（ファイル）を作ったりしていますが，実際裏では，UNIX(mac)やMS-DOS(windows)が動いてくれています。macのPCを使用している方は，UNIXをインストールする必要はありませんが，windowsをお使いの方はcgwin等をインストールして，UNIX環境下にする必要があります。
このページの先頭へ
1-3) 大型計算機とジョブ管理システム用のコマンド
NGSのデータは非常に膨大な容量で，それを解析するには多くのCPUやメモリが必要になります。しかしながら，Mac Proを使用してたとしても，最大CPUやメモリには限界があります。そこで，理化学研究所の京や東工大のTSUBAMEに代表するようなスーパーコンピューターが必要になります。スーパーコンピューターには，ジョブ管理システムが使用されており，それを使うために少し特有なコマンドを覚える必要があります。
このページの先頭へ
1-4) Docker
このページの先頭へ
1-5）UNIXとDOSのコマンドの違い
UNIXコマンドを覚えるとデータの編集等がとても楽になります。windowsを使用していると，UNIXコマンドが使用したくなる状況が時々あります。ここでは，UNIXとMS-DOSのコマンドの違いをご紹介します。
このページの先頭へ
1-6) FASTA形式とFASTQ形式の違いは？
塩基配列は、一般的にFASTA形式やFASTQ形式で書かれるとこが多いです。ここで、しっかりその形式の違いを学びましょう。
このページの先頭へ
ゲノム解析
とある生物のゲノム配列を決定して，興味ある形質がどの遺伝子に起因しているのか？もしくは，今後モデル生物としての研究の基盤のためゲノム配列を決定するのか？さまざまな理由からゲノム配列を決定したい研究者が存在していると思います。ここでは，私がとあるカビゲノムを決定した時の解析方法や印象を含めて書き留めておきます。
一般に，真正細菌レベル（～数十Ｍb）であれば，PacBio RS IIを用いたロングリードを使用すれば，完全に配列を決定することができます。さらに予算次第では，ゲノムサイズが比較的小さい真核生物（～600
Mbくらい？）まではPacBio RS IIを用いた方が長いゲノム配列が得られる印象です。しかしながら，1Gbを越えるようなゲノムサイズをもつ生物になると，お金もかかりますし，いろいろな問題も出てきます。こうした場合は，ショートリード
( mate-pairやpaired-end)を一緒に読んで，hybrid assemblyをした方が，戦略的に良いかもしれません。
Pacbio RS IIを使用できない方でも，illuminaのショートリードでもゲノム配列は繋がりますが，repeat regionやGCやATリッチの領域があると，やはり繋がりにくい印象を受けます。
2-1）ゲノムサイズの推定
ゲノムサイズを推定するために，これまでフローサイトメーターや蛍光顕微鏡などを用いて推定してきました。さらに最近では，NGSで得られた総塩基数とある配列の出現回数を計算することで，より正確に推定することが可能になりました。
このページの先頭へ
2-2）得られたリードの品質の確認
NGSからリードが得られたら，まずそのクオリティを確認します。ライブラリーの品質や機器の問題などでクオリティが低くなってしまっている場合があります。
このページの先頭へ
2-3）低品質の塩基の除去
低品質の塩基は，アセンブルするのに邪魔になります。そこで，アセンブルを行う前に低品質の塩基は除去する必要があります。
このページの先頭へ
2-4）De novo assembly
A.ショートリード用
使用したNGSの機器によって，アセンブルに使用するソフトウェアが異なります。イルミナのようなショートリードをアセンブルする場合は，以下のようなアセンブラーを用います。
B.ロングリード用
PacBio RS IIや454から得られるロングリードをアセンブルする場合は，以下のようなアセンブラーを用います。
- Celera WGA Assembler
サポートが終了し、canuに移行した。
- canu
これ一択 、PacBio RSIIとOxford Nanopore MinIONに対応、ショートリードの場合は、Celeraを使用する。
- Newbler
454専用。
- Mira
あまり使えないという噂
- Sprai
PacBio RS II専用。笠原さん＠東大が開発したアセンブラ，というかパイプライン。
ゲノム増幅したサンプルなどにも対応できる。Spariの中身でCeleraが動いている。
C. ショートリードとロングリード混在する場合
現在のところ，ショートリードとロングリード混在した状態でアセンブルようなアセンブラーは開発されていないと思います。以下のソフトウェアは，あくまでもアセンブルができるというだけで，品質の保証はできません。
このページの先頭へ
2-5）ゲノム配列の完全性の評価
アセンブルが完了したら，ゲノム配列の完成度を評価します。真核生物の場合，必ず保存している遺伝子が248個あります。それをアセンブルされた配列の中に存在するか確認することで，ゲノム配列の完成度を評価します。
このページの先頭へ
2-6）リピート配列のマスク
ゲノム配列が得られたら，次は遺伝子の予測です。まずはじめに，リピート配列や既知のトランスポゾンなどが，遺伝子を予測する上で邪魔になってしまうため，そういった領域をマスクします。
このページの先頭へ
2-7）遺伝子予測
ここでも，バクテリアと真核生物で解析手法が大きく異なります。バクテリアの場合は，とても良いソフトウェア（パイプライン）ができているので，それをご紹介いたします。一方で，真核生物は，さまざまな情報（近縁種の遺伝子情報，RNAseqの結果，リピート配列やトランスポゾン情報など）を駆使して，遺伝子を予測しなければなりません。
A. バクテリア
B. 真核生物
真核生物の遺伝子予測にも，いくつかのストラテジーがあります。
1. RNAseqの情報のみで，遺伝子予測する場合
RNA-seq解析でリードがマップされた領域を遺伝子領域として抽出します。利点としては，cufflinksを使用するだけなので解析にとても単純です。しかしながら，欠点としてサンプリング条件や組織によって発現していない遺伝子は予測できません。
2. 近縁種の遺伝子情報を用いて，遺伝子予測する場合（evidence-based 遺伝子予測）
- EXONERATE
既知の生物のタンパク質やCDS配列を用いて，得られたゲノム配列から似ている領域を探します。これも欠点として，生物特異的に持つタンパク質は検出できません。
- Spada
peptide proteinの検出に特化している。ナズナとメディカゴのデータはデフォルトで入っているが，その他に興味のあるpeptide sequenceを加えて検出することも可能。
3. 遺伝子モデルを使用して，遺伝子予測する場合（ab initio 遺伝子予測）
- AUGUSTUS
登録されている遺伝子モデルを参考にして，遺伝子を予測します。登録されていなくても，その後traningすることで，改善することができます。
4. さまざまな情報を使用して，遺伝子予測する場合
これが，最も精度よく遺伝子を予測することができます。まずは，1.のRNAseq解析からイントロン領域を得ます。2.から近縁種でタンパク質をコードしているエキソン領域を得ます。最後にリピート配列やトランプゾンなどの情報をすでに得ていると思います。これら情報すべてをAUGUSTUSにインプットして，遺伝子を予測させます。こうすることで，
ab initio 遺伝子予測とevidence-based 遺伝子予測から得られた遺伝子モデルを結合することができます。
このページの先頭へ
2-8）アノテーション付け
遺伝子が予測できたら，その遺伝子をさまざまなデータベースを用いて注釈付けしていきます。
DNAやタンパク質の配列比較
オントロジー解析
オーソログ解析
パスウェイ解析
タンパク質のドメイン検索
tRNA配列の検出
- tRNAscan
もし完全にゲノムが読まれているのであれば，20アミノ酸に対応するtRNAが検出されてくるはずです。
移行シグナル配列の検出
このページの先頭へ
RNA-Seq解析 - genome-based analysis -
-
RNA-seqは，遺伝子の発現量解析だけでなく，(新規)遺伝子の同定や遺伝子カタログ作成，isoformの検出など様々な解析が可能であり，非常に強力ツールである。RNAseq解析には大きく分けて，ゲノム配列の有無によって2つのストラテジーがある。ここでは，レファレンス配列がある場合の解析法の一例を紹介する。レファレンスがない場合の解析法については，次章で紹介する。
3-1）低品質の塩基の除去
はじめに，低品質の塩基領域を取り除きます。
このページの先頭へ
3-2）リードのマッピング
ゲノム配列が明らかになっている生物を解析する場合は，レファレンスゲノムにショートリードを張り付けます。対象生物が，原核生物の場合だとエキソン-イントロン構造を考慮する必要がないため，BWAやBowtieを使用します。一方で，真核生物の場合は，エキソン-イントロン構造を考慮する必要があるため，Tophatを使用します。
このページの先頭へ
3-3）マッピング結果の可視化
どのようにリードが，レファレンスにマッピングされたのか，3-2)の結果を可視化します。IGV以外使用したことがありませんが、使い方はどれも同じだと思います。
このページの先頭へ
3-4）リードのカウントから発現変動遺伝子の同定
サンプル間で遺伝子発現の違いのある遺伝子を同定するためには、初めに遺伝子領域にマップされたリードをカウントします。その後、サンプル間で得られたリード数や遺伝子の長さ・発現量の違いを正規化する必要があります。これらの操作を1つのプログラムで行う方法が、以降でご紹介するプランAです。また、これらの操作を行うのに一番良いソフトウェアを組み合わせて、実行する方法がプランBです。
A.カウントから正規化，発現変動遺伝子の検出まで解析するパイプライン
※正規化をFPKMで行っているせいか，発現変動遺伝子の検出感度が悪い印象を受ける。
B.カウント，正規化，発現変動遺伝子の同定をそれぞれベストのソフトウェアを使用する場合
B-1）カウント
B-2)正規化から発現変動遺伝子の検出
このページの先頭へ
3-5) 遺伝子群のクラスタリング
階層型クラスタリング
非階層型クラスタリング
3-6) 融合遺伝子の解析
細胞ががん化すると遺伝子変異により、遺伝子が融合することがあることが明らかになっている。そうした遺伝子を同定するため、いくつかソフトウェアが開発されているので、紹介します。
このページの先頭へ
RNA-seq解析 - De novo transcript-based -
ここでは，レファレンス配列が明らかになっていないような非モデル生物のRNAseq解析法についてご紹介する。
4-1）低品質の塩基の除去
はじめにアセンブルするため，低品質の塩基領域を取り除きます。低品質の塩基が含まれているとアセンブルの邪魔になってしまいます。
このページの先頭へ
4-2）De novo transcriptome assembly
次に高品質のリード配列を
de novo transcriptoome assembleすることによって，contig(塩基配列の断片)を得ます。
このページの先頭へ
4-3）マッピングとリードカウント, 発現変動遺伝子の同定
De novo transcriptの場合は，trinityに付属されているRSEMを用いて，contigにマッピングします。
このページの先頭へ
4-4）マッピング結果の可視化
どのようにリードが，レファレンスにマッピングされたのか，4-3)の結果を可視化します。IGV以外使用したことがありませんが、使い方はどれも同じだと思います。
このページの先頭へ
4-5）その他の解析
- CD-HITによる類似した配列の除去
- タンパク質配列の予測
De novo transcriptoome assemblyによって得られたcontigが，どんなタンパク質をコードしているのかORFを抽出してくれるソフトウェアです。tirinityに付属されています。
| transdecoder -t [fasta_file] -m [minimum protein_length] |
このページの先頭へ
リシーケンシング解析，エキソーム解析
5-1）低品質の塩基の除去
このページの先頭へ
5-2）マッピング
このページの先頭へ
5-3）SNPs/INDELsの検出
手っ取り早くSNPやIndelを検出するなら，samtoolsとその姉妹ツールのbcftoolsで良いかと思います。きちんと解析する場合は、GATKの方が性能が良いとされています。ただ、GATKは、商業ベースで用いる場合、ライセンスの購入が必要です。
このページの先頭へ
5-4) 変異の影響を確認（同義的置換？非同義的置換？）
検出されたSNPsやIDELsによって、タンパク質配列にどのような影響が現れるのか確認します。変異が同義的置換を起こすのか？もしくは非同義的置換なのか？
- SnpEff
SNPが同義的置換か非同義的置換かチェックするソフトウェア
- Annovar
このページの先頭へ
メタゲノム解析
6-1) はじめに
ちょっと書きたいことを書いておきます。
このページの先頭へ
6-2) DNAバーコーディング
お店でバーコードを読み込んで商品や価格を認識するように，ある特定の遺伝子領域のDNA配列を解析することで生物種を同定する方法です。解析対象の生物種によって用いる遺伝子領域が異なるので、目的に合わせて解析する必要があります。
このページの先頭へ
(6-3) 合成核酸の利用
真菌類の菌叢を調べる際、ITS領域を用いてアンプリコン解析する場合が多い。しかしながら、植物の根に共生する糸状菌を解析する場合、ITS領域を用いてしまうと、糸状菌だけでなく植物由来のDNAも増幅されてしまい、糸状菌由来の解析が困難になる場合がある。同様なことが、16S
rRNA領域でも起こる。葉に付着する細菌の菌叢を解析するため、16S rRNAのV4領域を増幅させると、細菌だけでなく、葉緑体やミトコンドリア由来のDNAも同時に増幅され、細菌の解析の感度が著しく低下する場合がある。そういった問題を解決するため、PNAやLNAなどの合成核酸が利用されてきている。
このページの先頭へ
6-4) アンプリコン解析（リードの重ね合わせ）
一般にDNAバーコーディングで用いられるDNA配列は、150-400bp程度です。Illumina社のNGSから得られる片リードは現在最長300bpであるため、片側のリードだけでは足りません。そこで、paired-endリード(リード1とリード2)をマージすることで、1本の長いリードにする必要があります。その後、この長いリードを用いて、菌叢解析やOTU解析を行います。
このページの先頭へ
6-5) アンプリコン解析（菌叢解析, 群集解析，OTU解析）パイプライン
16S rRNAやITS領域、18S rRNAについては、Qiime用のデータベースが整っているが、それ以外のCOIやMiFishのターゲット領域である12S
rRNAはデータベースが整っていない。そのため、ここではUSEARCHを用いた方法をご紹介する。
このページの先頭へ
6-6) 環境DNA
最近話題の環境DNAについて、まとめました。
このページの先頭へ
6-7) メタゲノム解析
菌叢解析は，さまざまな微生物の存在を明らかにすることはできるが，機能に関する情報を得ることができない。そこで，メタゲノム解析は、環境中で複雑に関わり合っている微生物のゲノム配列（遺伝子断片）を大雑把に決定し，その後アノテーションを付ける解析である。
メタゲノムプロジェクト
すでにいくつかのプロジェクトによって、多くの種遺伝子カタログが明らかになっている。これらの情報を利用して解析することも可能である。
メタゲノムアセンブラ
通常の単離された微生物とは異なり，メタゲノム解析では多種間で共通した配列（chimeric node）が存在するため，通常のアセンブラと異なる方法でアセンブルを行う必要性がある。同じサンプルでいくつか結果を比較したところ、SPAdesかIDBA-UDが比較的に良い成績であった（答えが分からないので、どちらが良いのかはわからない）。ただ、SPAdesは、IDBA-UDと比較して、contigを分けようとする傾向が強い印象を受ける。
このページの先頭へ
6-8) 遺伝子予測
メタアセンブル後に遺伝子予測を行います。バクテリアの遺伝子予測で紹介したprokkaやGeneMarkは，メタゲノム用にも使用できます。実行時間に大きく差があり，数によってはprokkaの場合だと数日かかる実行時間が，MetaGeneMarkを使用すれば，数分で完了することが可能です。
このページの先頭へ
6-9) アノテーション
得られたドラフトゲノム配列の中にどのような遺伝子があるのか，データベースと比較してアノテーションを付けていきます。
このページの先頭へ
生物種の同定
得られた配列の中にどのような生物がいるのか確認するため，rRNAの配列を検出するソフトウェアです。予測された配列をNCBIでBLASTすることで生物種を同定します。
- Metaxa
Small Subunitをコードする遺伝子の配列を検出するソフトウェア
このページの先頭へ
6-10) メタトランスクリプトーム解析
Trinityのようにアセンブルしてから解析するというよりは、Paired-endリードを繋げて、既知のデータベースと比較する解析方法です。
このページの先頭へ
6-11) 生態学的なデータ解析
このページの先頭へ
ジェノタイピングシーケンシング解析～RAD-seq解析、MIG-seq解析、GRAS-Di解析～
7-0) 各解析の違いについて
7-1) 解析パイプライン
RADシーケンス(RAD-Seq)は，Restriction Site Associated DNA Sequenceの略で，制限酵素認識サイトの周辺領域から塩基多型を検出する手法です。この方法は，ゲノム配列を制限酵素で切断し，切断した部位から
100-300 bp 程度の配列を読み，その配列から SNPs を探索する。
一方で、MIG-seqは、multiplexed ISSR genotyping by sequencingの略で、Simple Sequence
Repeat (SSR)の間の配列を増幅し、塩基多型を検出する手法です。RAD-seq法と異なり純度の高いDNAを必要と必要とせず、またPCRベースなのでDNA量も少なくて済む。ただ、原理上RAD-seq法より得られる多型の数は少ない。[
PDF]
これら解析により得られたSNPマーカーは，遺伝子座間の距離の測定、高密度遺伝子連鎖地図の構築やQTLマッピングなどに用いることができる。
このページの先頭へ
7-2) データの補完
RAD-seq解析は一度に多くのSNPsがゲノムワイドに同定/ゲノタイピングすることが可能だが，一方で欠測が多いという問題点がある。それを補うために，imputation(データの補完)を行う必要があります。
このページの先頭へ
7-3) 連鎖地図の作成
連鎖地図を作成するには，以下のソフトウェアを使用します。onemapは無償で使用することができます。一方でjoinmapは，一部無償で使用できますがさまざまな制約があるようです。それを解除するためには、お金が必要になります。
このページの先頭へ
7-4) 集団遺伝学的な解析
このページの先頭へ
7-5) 系統解析
このページの先頭へ
分子マーカーの開発
8-1) 分子マーカー
このページの先頭へ
8-2) SSRマーカーの開発
Simple Sequence Repeat (SSR)マーカーを予測するソフトウェアです。
このページの先頭へ
8-3) dCPASマーカーの開発
dCAPSマーカーを予測するソフトウェアです。
このページの先頭へ
8-4) SNPマーカーの開発
このページの先頭へ
8-5) フラグメント解析
予測されたマーカーが正確に機能するかサンガーシーケンサーで確認を行います。
このページの先頭へ
メチル化解析
whole-genome bisulfite sequencingにより、ゲノムレベルでメチル化シトシンの有無を検出することが可能です。バイサルファイト処理は、メチル化されていないシトシン(C)は、ウラシル(U)に変換します。一方でメチル化したシトシン(mC)は、Uに変換させません。この処理によりメチル化の有無を区別することができます。あとは、NGSで配列を決定するのみです。
9-1)メチル化解析
このページの先頭へ
その他
10-1）インストールや初期設定に必要なもの
プログラミングを実行する際に必要なものをインストールしてください。何か警告が出た時にインストールするのが
いいかと思います。
このページの先頭へ
10-2）ファイル形式
次世代シークエンスの結果を解析する上でいくつかの形式があります。いくつかはバイナリーファイル（人間が理解できない機械語のようなもの）で，今回は人間が理解できる形式ファイルのsamとbedファイルについて紹介します
このページの先頭へ
10-3）ファイル・データ変換ツール
このページの先頭へ
10-4）DDBJ登録
- DRA登録方法
NGSから得られた生リードを登録します。
- MSS登録方法
生リードをアセンブルして得られたcontig/scaffoldを登録します。
このページの先頭へ
10-5）Circosの使い方
このページの先頭へ
10-6）生物の分類
このページの先頭へ
10-7）大量プライマー設計
このページの先頭へ
