シークエンス解析はデータが大量にあってとても興味深いですが、専門性が高いソフトが使われているので敷居が高いです。
あまり知らないフトをたくさんPCにインストールするのは、管理が大変そうなので躊躇われます。
しかし、Biocondaを使えば全てAnacondaの仮想環境上で行うことが出来るようです。
そこで今回はBiocondaを用いてRNA-seqの基本的な部分を仮想環境上でやってみます。

＊Biocondaの論文(Nature Methods)は以下です。
Grüning, Björn, Ryan Dale, Andreas Sjödin, Brad A. Chapman, Jillian Rowe, Christopher H. Tomkins-Tinch, Renan Valieris, the Bioconda Team, and Johannes Köster. 2018. “Bioconda: Sustainable and Comprehensive Software Distribution for the Life Sciences”. Nature Methods, 2018 doi:10.1038/s41592-018-0046-7

解析の概要

RNA-seq解析の概要としては、以下のようにデータファイルが生成されていきます。

SRAファイル・・・データベースからダウンロードした時のファイル。保管用のファイル形式。
　　　↓
FASTQファイル・・・シーケンサでのシークエンスの結果(各リードの塩基配列)が入っている。
　　　↓
BAMファイル・・・リードをリファレンスにマッピングした結果が入っている。(生成しないソフトもある)
　　　↓
カウントデータ(CSV形式など)・・・転写産物ごとのリードカウント（発現量を表す）が表形式で入っている。

カウントデータにまでしてしまえば情報量は大幅に減りますが、後は通常のデータ解析のように解析できると思います。

(上の画像はhttps://bi.biopapyrus.jp/rnaseq/ より引用させて頂きました。)

0.環境構築

環境は
Ubuntu 16.04
Anaconda
で行いました。
(＊シークエンス解析ではWindowsはあまり使われないです。)

まずは仮想環境を構築しておきます。

conda create -n <仮想環境名> python=3.6 anaconda
source activate <仮想環境名>

Biocondaからのソフトのインストールは-cでチャンネルをすれば行えます。

conda install -c bioconda <インストールしたいソフト名>

次から必要なソフトのインストールをして、解析を行っていきます。

データの用意

データは、DRAやAOEなどからダウンロードすることができます。
今回はSRR1550981.sraをダウンロードしました。データ元の論文はこちらです。

解析にはリファレンスとなるデータも必要です。
リファレンスデータはGENCODEやENSEMBLからダウンロードできます。
今回のデータはヒトなので、リファレンスとしてもHomo_sapiensのものをダウンロードします。

1.FASTQファイルの生成

sra-toolsにはSRAファイルを扱う様々なツールが入っているので便利です。
pigzはgzへの圧縮を並列化して高速化するソフトです。
ソフトをインストールします。(fasterq-dumpを使いたいためsra-toolsのバージョンを指定しています。)

conda install -c bioconda sra-tools==2.10.0
conda install pigz

ダウンロードしたSRAファイルからFASTQファイルを生成します。
fasterq-dumpでは、ペアエンドを自動で認識して2つのFASTQファイルが生成されます。

fasterq-dump SRR1550981.sra

これで２つのFASTQファイルが生成されました。
headコマンドで1つ目のファイルの最初の10行を見てみます。

$ head ./SRR1550981.sra_1.fastq
@SRR1550981.sra.1 CHAUSSABELL-ILL:183:D16YPACXX:5:1101:2044:1949 length=50
CCCGCGTCTGCGGCTTCATGANGGTGCANNGNGCGNTGNGNANGTTGCAC
+SRR1550981.sra.1 CHAUSSABELL-ILL:183:D16YPACXX:5:1101:2044:1949 length=50
CCCFFFFFHHHHHIJJIJIIJ#0??FHI######################
@SRR1550981.sra.2 CHAUSSABELL-ILL:183:D16YPACXX:5:1101:2122:1981 length=50
CCAAACTGTTCCCGGGCAGCACTGAACATGGGCTCCTGGATGTCCGTGTA
+SRR1550981.sra.2 CHAUSSABELL-ILL:183:D16YPACXX:5:1101:2122:1981 length=50
CCCFFFFFHHHHHJJIHJJJJJJJJJJJJJJJJJJJCHIIJIIIJJGIII
@SRR1550981.sra.3 CHAUSSABELL-ILL:183:D16YPACXX:5:1101:2657:1962 length=50
CTGTATCAGGCACTCAAGTTGTGCAGGGACAGATCCAGACACTTGCCACC

リードの塩基配列(ATGC)が入っていることが分かります。

FASTQは非常に容量が大きいファイルです。今回のデータだと1つの大きさで5.9GBもあります。
そのため、pigzで圧縮しておきます。

pigz SRR1550981.sra_?.fastq

これで2つともが1.1GBまで圧縮できました。
圧縮したファイル形式(.fastq.gz)であってもシークエンス解析ソフトは入力として対応しています。

2.カウントデータの生成

方法は2通りあるので、2つの場合に分けて書きます。

2.A マッピングをする方法

この方法では、BAMファイルを生成して、その後に発現量をカウントします。
今回はマッピングにはHISAT2を使います。（他にはSTARを使う方法もありますが、PCのメモリが多くないとできません。）
ソフトをインストールします。

conda install -c bioconda trim-galore hisat2 samtools stringtie deeptools igv

2.A.1 マッピングデータの生成

まずはtrim-galoreでFASTQデータのトリミングを行います。

trim_galore --paired --fastqc SRR1550981.sra_1.fastq.gz  SRR1550981.sra_2.fastq.gz

トリミングされたFASTQデータと、FastQC(リードの品質管理)の結果も出力されました。

HISAT2でマッピングします。
まずは、リファレンスデータからインデックスを作成します。リファレンスデータは、gzで圧縮されたものは対応しておらず解凍する必要があるので注意です。（インデックスを作らずにHISAT2のサイトからダウンロードすることも可能です。）

hisat2-build <リファレンスのfastaファイル> gene_idx

これで名前に「gene_idx」がついたインデックスファイルが8つができました。
マッピングを行います。


hisat2 -1 SRR1550981.sra_1_val_1.fq.gz -2 SRR1550981.sra_2_val_2.fq.gz -x gene_idx -S mapping.sam

マッピングできました。
マッピング結果はsamファイルに保存されています。
samtoolsによって、ソートしてBAMファイルに変換します。その後にインデックス(.baiファイル)も作成します。

samtools sort mapping.sam -O bam -o sort.bam
samtools index sort.bam

BAMファイルはIGVというGUIツールで可視化することが出来ます。
まずはdeeptoolsのbamCoverageで処理をして、そこで出力されたファイルをigvでロードすれば可視化することができます。

bamCoverage -b sort.bam -bs 1 -o output.bigwig
igv

IGVで可視化が出来ました。
＊samファイルやbamファイルなどはサイズが大きくて処理に時間がかかります。その場合は-pのオプションで並列処理をすると速くなります。

2.A.2 カウントデータの生成

stringtieよってBAMファイルからカウントデータを生成できます。
この時にリファレンスのアノテーションファイル(gtf形式)が必要です。
カウントデータの生成を生成します。

stringtie sort.bam -G <gif形式のリファレンスのアノテーションデータ> -o output.gtf  -A output.csv

これでできました。
カウントデータはoutput.csvの中に入っています。(もうひとつの出力のoutput.gtfにも入っています。)
カウントデータの中身を見てみます。

$ head ./output.csv
Gene ID	Gene Name	Reference	Strand	Start	End	Coverage	FPKM	TPM
ENSG00000223972.5	DDX11L1	chr1	+	11869	14409	2.747866	1.533925	2.869029
ENSG00000227232.5	WASH7P	chr1	-	14404	29570	24.237690	11.716399	21.914165
ENSG00000278267.1	MIR6859-1	chr1	-	17369	17436	16.619198	8.033651	15.026013
ENSG00000243485.5	MIR1302-2HG	chr1	+	29554	31109	0.046033	0.022252	0.041620
ENSG00000284332.1	MIR1302-2	chr1	+	30366	30503	0.000000	0.000000	0.000000
STRG.1	-	chr1	+	11869	14409	3.151109	1.533925	2.869029
STRG.2	-	chr1	-	14632	24891	10.647043	15.779980	29.514622
STRG.3	-	chr1	-	17369	17436	16.619198	8.033651	15.026013
ENSG00000237613.2	FAM138A	chr1	-	34554	36081	0.072601	0.035095	0.065641

2.B マッピングしない方法

前項での方法は手順が多くて時間もかかりました。
解析においてカウントデータのみが欲しい場合は、マッピングしない方法の方が簡便です。
この方法ではマッピングをしないため、BAMファイルが生成されないです。
その代わりに高速で行うことが出来ます。
今回はKallistoを用います。(他にはSalmonでも同じことができます。)
ソフトをインストールします。

conda install -c bioconda kallisto

まずはKallistoのインデックスファイルを作成します。

kallisto index <リファレンスのfastaファイル> -i gene_idx

FASTQファイルから一気にカウントデータを生成します。(-oで出力先の指定が必須)

kallisto quant SRR1550981.sra_1.fastq.gz  SRR1550981.sra_2.fastq.gz -i gene_idx  -o ./

これだけでカウントデータが出来ました。しかもとても速いです。
カウントデータの中身を見てみます。

$ head ./abundance.tsv
target_id	length	eff_length	est_counts	tpm
ENST00000632684.1	12	6	0	0
ENST00000434970.2	9	5.5	0	0
ENST00000448914.1	13	7	0	0
ENST00000415118.1	8	4.5	0	0
ENST00000631435.1	12	6	0	0
ENST00000390583.1	31	14	0	0
ENST00000390577.1	37	16.4615	1	4.28747
ENST00000451044.1	17	8.5	0	0
ENST00000390578.1	31	14	0	0

3.カウントデータの解析

今回は割愛します。
テーブルデータなのでPythonやMATLABも使えるとは思いますが、Rを使うのが一般的のようです。
実はRのインストールもcondaで仮想環境上で行うことが出来ます。

conda install r
conda install rstudio

RにはBioconductorというシークエンス解析用のライブラリが充実しているのでそれが使えます。

また他には、Web上で解析を行うiDEPを使うという方法もあります。

最後に

最近のシークエンス解析は以下の2点でとても便利になったと思います。
１．全てを仮想環境上で行えること
２．解析にかかる時間が非常に短くなったこと
（昔TopHat2やcufflinksを使ったことがありましたが、とても時間がかかって大変でした。）

また有意義な解析をおこなうには、NGSの原理についても勉強する必要があると思います。
NGSは大量のシークエンススループットを実現できますが、精度は改善されたわけではなく、現在でもサンガー法を使うことがあるそうです。(参考：https://www.megabank.tohoku.ac.jp/genome/archives/132)

これからもシークエンス解析の勉強していきたいと思います。

＊参考文献
次世代シークエンサーDRY解析教本
RNA-Seqデータ解析　WETラボのための鉄板レシピ