WSLを使ってゲノムのマッピングから比較ゲノム解析をしてみよう

初めまして、ウェット研究室から来たものです。
今回は、自分の勉強を兼ねて今までのドライ解析の勉強のメモをさせてもらいます。
正確、かつ、まともな情報が欲しい方は、この記事のキーワードで調べながら、先行する様々な論文や素晴らしいサイト様で確認してください。
また、比較ゲノム解析は資料が英語なことも多いので頑張ってください。

今回は多くのサイト様を参考にさせてもらっていますので、意図せずに盗用しているようなことがあれば、該当部分を修正します。

改善点などはもちろん受け付けておりますので、お時間がありましたら、ぜひともご指摘くださるとありがたいです。
特に、引用が不完全である、不正確であるなどの指摘は迅速に対応させていただきます。

#本記事の対象
・ゲノム解析に興味があるけれどドライ解析について分からない方
・それなりにいいパソコンとネットワーク回線を持っている方
・比較ゲノム解析に興味がある方
・がんばりや

#環境
Windows 10のノートパソコン（4CPU 8コア メモリ16GB）
容量の結構あるハードディスク
脆弱なネット回線
WSLのソフト（Windows for Linux）（Ubuntu 18.04 & Debian・マイクロソフトストアでダウンロードしてください）

#今回やること
・ゲノムの下処理
・DNAリードのファイル（sra FASTQ）の処理
・bwa-memでのマッピングなど
・Variant calling（欠損データのImputationは省略します）によるVCFファイルの作成
・VCFファイルの解析の一例

#前提条件として
公共データベースにこれがあれば解析は出来ます。
・目的生物のゲノム配列ファイル（拡張子は.fa、.fnaなど、NCBI等で手に入るはずです）
・近縁種などの目的生物のゲノムDNAリード（.sraファイル）

#まずは必要なファイルやプログラムをインストールしよう
まず、WSL（Ubuntu 18.04を普段自分は使っています）をマイクロソフトストアからダウンロードします。具体的な処理については他の記事などを参考にしてください。
名前とパスワードの設定が必要です（特にパスワードはよく使うためめ覚えやすいものにしてください）。

自分はGドライブにデータを入れているため

cd /mnt/g

というコマンドを入力してGドライブでの操作を行っています。
また、新しいプログラムを入手するためには、

sudo apt install XX #パスワード入力が必要で、sudoをしないと動きません。セキュリティに注意してください。

と入力してください（python-pip・git・condaなどはよく使います）。

また、Githubからデータを使わせてもらうときには、Githubの「Cloneと書かれているボタン」をクリックすると、末尾に「.git」と書かれているURLが出るため、そのURLを使って

cd /mnt/g #Gドライブなどダウンロードしたい場所に移動する
git clone XX.git

として、あとはGithubのホームページなどに書かれている指示に従えばプログラムを使えるようになります（PATH=環境変数に注意してください）。

また、pythonで書かれているプログラム入手は、「miniconda3」が便利です。これは先ほどのsudo aptではなく、公式ホームページから入手します。

まず、以下のホームページにアクセスしてMiniconda3 Linux 64-bitをダウンロードすると、末尾に.shという拡張子がついたファイルが手に入ります。
https://docs.conda.io/en/latest/miniconda.html

このファイルを自分のインストールしたい場所（自分ならGドライブ）にコピーしてから、インストールします。プログラムの指示にいろいろ従ってください。

cd /mnt/g #Gドライブなどダウンロードしたい場所に移動する
./Miniconda3-latest-Linux-x86_64.sh #移動した場所（Gドライブ）に配置した.shファイルを実行するという意味になります。

うまくWSLの環境変数（PATH）が入ったなら、WSLを再起動すると、「conda」コマンドでいろいろPythonで書かれたプログラムをインストールできるようになります。（Pythonでは似たような方法でpipコマンドも使いますが、pipとcondaを共存させるのは禁忌とされています。）

condaでのプログラムの入手は以下のようになります。

conda install -c bioconda xpclr #例としてXP-CLRをダウンロード

#ゲノムDNA配列の入手
NCBIの公共データベースで様々な生物のゲノムが公開されていますが、最新の論文や目的の生物専門のデータベースで「Assembly」「Genome」と書かれているものからデータを取得することをお勧めします。末尾は「.fa」か「.fna」になっているものが正解です。
　lessコマンドなどで確認すると、最初の行に「>配列名」、次の行に「ATGCなどのゲノム配列」になっているはずです（メモ帳で開くとメモリがパンクするのでやめましょう）。

#DNAリードの入手
　ゲノムへのマッピングのためには、FASTQファイルが必要です。まずはそれを入手しましょう。
　こちらもNCBIのデータベースで公開されています。
　https://www.ncbi.nlm.nih.gov/
　素晴らしい研究者の方々が、研究の成果の元となった、DNAシーケンシングのデータを公開されています。NCBIの検索窓でSRAを選んで、「目的の生物 AND resequencing」などで検索すれば、様々なデータが見つかると思います（無い生物の場合はあきらめるか自作してください、RNA-seqのデータと間違えないようにしてください）。

　公共データベースに公開されているデータは、数GBあるSRAファイルで保存されています。そのため、ダウンロードが必要です。専用ソフトとして、SRA TOOLKITが一般的に使われています。以下のサイトからソフトをダウンロードしてください。
https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software
ソフトの解凍、インストール、PATH（環境変数）の通し方などは他の方の記事を見てください。

cdコマンドでsraファイル（数GB）をダウンロードする場所に移動してから、prefetchでデータ取得、fastq-dumpで解凍することをお勧めします（fastq-dumpや高速なfasterq-dumpからそのままデータを入手することもできます。）。

sraファイルを入手したら、以下のようにFASTQファイルにします。
FASTQファイルはそれぞれ10GBあるので、ハードディスクの容量には注意してください。
DNAシーケンシングデータにはsingle-endとpair-end（2つセット）がありますが、自分はよくpair-endを使うためそれをメインに解説します。

cd /mnt/g #ディレクトリ(=フォルダ)の移動
sra=SRR0000000 #目的のSRAのAccession Number、DRRの場合もある
./sratoolkit.2.10.8-ubuntu64/bin/prefetch-orig.2.10.8  --progress 1 ${sra} #ファイルをダウンロードします
./sratoolkit.2.10.8-ubuntu64/bin/fastq-dump-orig.2.10.8 /mnt/g/${sra}/${sra}.sra \
 -M 25 -E --skip-technical --split-3 -W #DRAというデータベースでのFASTQ解凍を参考にしました。
#https://www.ddbj.nig.ac.jp/faq/ja/read-number-fastq.html
#sed -e 's/forward//g' -e 's/reverse//g' ${sra}_1.fastq > ${sra}_1_sed.fastq & sed -e 's/reverse//g' -e 's/forward//g' ${sra}_2.fastq > ${sra}_2_sed.fastq #リード名がFASTQ間で一致して後で紹介するbwa-memにエラーが出たため、これをしました。通常は必要ありません。ファイルのリード名の変更を意味します。

#リードの確認とアダプタートリミング
リードの質の確認は、FASTQCというソフトを使ってください。リードの質を確認できます。
また、DNAシーケンシングにはアダプター配列がついています。そのため、trim-galore!というソフトでアダプターを取り除いてください。trim-galore!は一般的なアダプター配列なら自動で取り除いてくれます。
どちらもcondaでダウンロードできます。

(参考)
https://bi.biopapyrus.jp/rnaseq/qc/fastqc.html
http://kazumaxneo.hatenablog.com/entry/2017/06/21/111007

#ゲノムデータへのマッピング

　ついにFASTQファイルとゲノムデータがそろいました。まずはゲノムデータに下処理をして、それからマッピングします。マッピングしたデータはsam・bamファイルで出力されます。また、bamファイルは使いやすいようにsortとindexします。

sudo apt install bwa #マッピングに使います
bwa index Genome.fa #インデックスファイルができます。

　次に、マッピングを行います。自分の環境だと一つのFASTQペアをゲノムにマッピングするのに数時間かかりました。

bwa mem -t 8 Genome.fa SRR0000000_1.fastq SRR0000000_2.fastq > SRR0000000.sam #tはパソコンのコア数です。
samtools view -bS -@ 8 SRR0000000.sam > SRR0000000.bam #こちらも@もコア数です。bamは圧縮されたファイルです。
samtools sort SRR0000000.bam  > SRR0000000.sort.bam 
samtools index SRR0000000.sort.bam     

　パソコンも皆さまもお疲れさまでした。これでゲノムにDNAシーケンシングのデータがマッピングできました。さて、「比較ゲノム解析」にはこれを複数やる必要があります。頑張りましょう。

#Variant calling（例：Platypus）
　さて、比較ゲノム解析を行うには、DNAの個体データ間（SNPsとIndels）の比較を行う必要があります。複数のbamファイルを使うため、用意してください。
　また、Variant callingにはGATKなど様々なプログラムが研究・配布されています。自分は高速かつIndels検出に強いPlatypusを使っています。そのため今回はPlatypusで説明します。
　まずは、自分の使いたいプログラムをダウンロードしましょう。

（なぜかUbuntuで使えなかったので、WSLでもdebianのpython2のcondaを使っています。）

conda install cython
git clone https://github.com/andyrimmer/Platypus.git
cd Platypus
make 

インストールしたらVarinat callingの開始です。

 ./Platypus/bin/Platypus.py callVariants --bamFiles=A.sort.bam B.sort.bam C.sort.bam \
 --refFile=Genome.fa \
 --assemble=1 --nCPU=8 \
 --verbosity=1 \
 --output=result.vcf  --regions=Chromosome_1 #一度に1つの染色体しかできません。

（参考）
http://kazumaxneo.hatenablog.com/entry/2017/03/27/100759

#Imputation
　作ったVCFデータも完全ではなく、欠損データがあります。これはBEAGLE 5.0というソフトでデータの補完を行います。
（難しいので説明は省きます。ごめんなさい）
https://faculty.washington.edu/browning/beagle/b5_0.html

#ついに比較ゲノム解析
さて、VCFファイルさえ手に入ればもう比較ゲノム解析は可能です。今回は集団間比較の$F_{ST}$とXP-CLRを例としてやってみましょう。

#Weir and Cockerham's FST
　これは、集団間の中でどれだけゲノムの変異の頻度が分かれているかを示すための指標で、「上位N%」の部分が集団化の差を示すことに使われています。
マッピングデータさえあれば、簡単にこの値を計算できます。この値が大きければ集団間の差は大きいということらしいです。
詳しくはマニュアルを見てください。
http://vcftools.sourceforge.net/man_latest.html

conda install -c bioconda vcftools
vcftools --vcf result.vcf --weir-fst-pop population_1.txt --weir-fst-pop population_2.txt

#XP-CLR
　これは、染色体ごとの、集団間で選択圧に差がある領域（Selective sweeps）を探すための検定を行う手法です。この数値が高い領域（20000bpなどの窓を使う）は集団間で選択圧が異なることを示すということのようです。Pythonでプログラムが動くこと、$F_{ST}$とは別の理論で基づいていることなどから、$F_{ST}$とセットで使ってみてはいかがでしょうか。結果はテキストファイルで出力されます。

Chen, H., Patterson, N., & Reich, D. (2010). Population differentiation as a test for selective sweeps. Genome research, 20(3), 393–402. https://doi.org/10.1101/gr.100545.109

#conda版は古いため、Githubから使いましょう
#conda install -c bioconda xpclr
xpclr --input result.vcf \
 --samplesA SRA000000 SRA000001 \
 --samplesB SRA000002 SRA000003 \
 --out /mnt/g/OUT --chr Chromosome_1

#その後のお話
　染色体ごとに$F_{ST}$やXP-CLRが計算出来たら、あとはそれをRのqqmanでマンハッタンプロットのようなものを書く、ビッグデータのように解析をする、これらの統計量が高いのはどの遺伝子領域かを調べるなど、いろいろできます。詳しくは興味のある分野の比較ゲノム解析論文を読んで真似してください。もちろん解析には努力が必要ですが…
　ウェット研究者ではLinuxに慣れていない方も多いですが、少しでもドライ解析が広まればより研究の幅が広がると思うので頑張ってください。応援しています。

#注意点
リファレンスゲノムは完全ではない可能性があります。
具体的には、「リファレンスゲノムの生物のゲノムのエキソン部分が、近縁種ではかけていてマッピングできないのに気づかない」などの状況があります。その場合、NCBIのBLASTをかけるなど、近縁種の配列と一致するかどうかを調べるなどをして、いろいろ調べてください。どうかご注意ください。

#余談
ダウンロードしたプログラムをエクスプローラーで目的地に動かして実行と、Permission Deniedといわれることがあります。
自分はchmodで権限を与えているのですが、それで大丈夫なのかよくわかっていません。