はじめに
RNA-seq（RNAシーケンシング）は、サンプル中の遺伝子の発現レベルを分析するために使用される方法です！
近年ではRNA-seqが安価で行えるようになり、非常に盛んに行われる手法になりました！
それに伴ってSequence Read Archive (SRA）などのデータベースに登録されているRNA-seqの実験データの数もかなり増えました！
SRAなどの公共のデータベースに登録されている実験データを使用すれば、実験を行わずとも目的の解析を行うことができる時代になっています！いい時代ですね！
ただ、RNA-seq解析の手法には本当にたくさんの種類があって、「どれをどんな風に使って解析したらいいのか分からない！ 」という方もたくさんいらっしゃると思います！
今回の記事ではRNA-seq解析の一連の流れと、その実践方法についても徹底的に解説していきます！
今回の一連の解析全てを実行可能なDockerイメージを作りました！こちらをぜひご使用ください！https://hub.docker.com/r/koji11235/rnaseq_env/
RNA-seq全体の流れ
https://bioinformatics-core-shared-training.github.io/RNAseq-R/slides/rnaSeq_May2017.pd
まずはRNA-seq解析全体の流れについてご紹介します！
解析自体の目的や、使用可能な計算リソースによって各フェイズで使用するソフトウェアなどが変更されることはありますが、以下の大まかな流れは基本的には変わリません！
- 品質管理（Quality Control (QC)）：
Input： 生シーケンスリード（FASTQ ）
Output: クリーンアップされたシーケンシングリード（FASTQ） - 参照ゲノムへのマッピング(Reference genome mapping)：
Input： クリーンアップされたシーケンシングリード（FASTQ）
Output: リードの参照ゲノムへのアラインメント結果（BAMまたはSAM） - 遺伝子発現の定量化(Quantification)：
Input： リードの参照ゲノムへのアラインメント結果（BAMまたはSAM）
Output: 遺伝子発現量のカウントファイル（タブ区切りテキスト形式） - 発現変動遺伝子解析(Differential expression analysis)：
Input: 遺伝子発現量のカウントファイル（タブ区切りテキスト形式）
Output:発現変動した遺伝子や転写産物（タブ区切りのテキスト形式） - 機能アノテーション(Functional annotation)：
Input: 発現変動した遺伝子や転写産物（タブ区切りのテキスト形式）
Output: GO, パスウェイエンリッチメントなどの結果（形式は様々）
RNA-seq解析の実践！
全体の流れを理解したところで、早速RNA-seq解析を実践していきましょう！
使用するデータ: airway
Rのパッケージにもなっているairwayデータセットを題材として使用します！概要を簡単にご説明します！
dexamethasone（Dex）という強力な合成グルココルチコイドの気道平滑筋（airway smooth muscle (ASM)）細胞への作用を調べることを目的として行われたRNA-Seqのデータです！
サンプルはDexを投与したサンプルと投与していないものがあります。（今回は単純化のため使用しませんが、β2-agonist Albuterolを投与したサンプルもデータベース（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778）には登録されています。）
今回使用するのはDexを投与した4つ（SRR1039509, SRR1039513, SRR1039517, SRR1039521）と投与していない4つ（SRR1039508, SRR1039512, SRR1039516, SRR1039520）の全部で8サンプルです！
今回の解析の目的は元論文と同じで、Dexを投与して発現量が変動する遺伝子とその機能について調べることにしましょう！
データセットの詳細は以下の論文をご参考ください！
Himes BE, Jiang X, Wagner P, Hu R, Wang Q, Klanderman B, et al. (2014) RNA-Seq Transcriptome Profiling Identifies CRISPLD2 as a Glucocorticoid Responsive Gene that Modulates Cytokine Function in Airway Smooth Muscle Cells. PLoS ONE 9(6): e99625. https://doi.org/10.1371/journal.pone.0099625
今回の解析パイプライン
今回の解析パイプラインに使用するソフトウェアは以下の通りです！
- 品質管理（Quality Control (QC)）：FastQC
- 参照ゲノムへのマッピング(Reference genome mapping)：STAR
- 遺伝子発現の定量化(Quantification)：RSEM
- 発現変動遺伝子解析(Differential expression analysis)：edgeR
- 機能アノテーション(Functional annotation)：clusterProfiler(パスウェイ解析, Gene Ontology解析用)
docker imageの紹介
今回実践する一連の解析を実行可能なDocker imageをご用意しました！
よろしければご使用ください！
https://hub.docker.com/r/koji11235/rnaseq_env/
今回使用した一連のソフトウェアを使用できるDocker Imageです！
docker pull koji11235/rnaseq_env
docker run –rm -v ${HostPath}:/home -p 8787:8787 -e PASSWORD=112358 rnaseq_env
rocker/verseをベースに作成しているので、シェルに入ってもらって、rstudio-server start コマンドを実行し、ブラウザでhttp://localhost:8787/に接続すればRstudioを使用することもできます！
ぜひご使用ください！
RNA-seq解析 ~Phase0~ データの準備
fasterq-dumpでSRAからfastqファイルのダウンロード
まずはデータの準備です！
シーケンスリードアーカイブデータベース（SRA）からfastqファイルのダウンロードをしましょう！
SRA Toolkitの中に含まれる、fasterq-dumpというコマンドを使用してfastqファイルのダウンロードをします！
実際のコードは以下です！
# Download fastq files
srrids=(SRR1039509 SRR1039513 SRR1039517 SRR1039521 SRR1039508 SRR1039512 SRR1039516 SRR1039520)
output_dir=data/seq/
for srrid in ${srrids[@]}; do
output_file_1=${output_dir}${srrid}_1.fastq
output_file_2=${output_dir}${srrid}_2.fastq
fasterq-dump -O ${output_dir} --split-files --threads 4 ${srrid}
gzip ${output_file_1}
gzip ${output_file_2}
done
--split-files オプションを使用して、ペアエンドのリードを別々のファイルとして保存します。
また、保存されるファイルは圧縮されていないfastqファイルですのでファイルサイズがとても大きいです。ですので、gzipコマンドを使用して圧縮します。
gzipコマンドを２回使用しているのは --split-files オプションを使用して、ペアエンドのリードを別々のファイルとして保存したためです。
もしシングルエンドのfastqファイルでしたら、${srrid}.fastqという名前のファイルが保存されますのでgzip ${output_dir}${srrid}.fastqとして下さい。
fasterq-dumpとfastq-dumpについて
fasterq-dumpはこれまでに非常によく使用されていたfastq-dumpの後継で、より高速に動作するようになっています。
fastq-dumpはまだサポートされてはいますが、後々非推奨になると公表されていますので、fasterq-dumpを使用するのをお勧めします！
fasterq-dumpの使用方法については以下をご参照ください。
https://github.com/ncbi/sra-tools/wiki/HowTo:-fasterq-dump
RNA-seq解析 ~Phase1~ 品質管理（Quality Control (QC)
まずは生シーケンスリードの品質チェックをします！
品質チェックを行うためのソフトウェアには非常にたくさんの種類がありますが、今回はFastQCを使用します！
FastQCはWindowsでもMacでもLinuxでも使用することができます！
また、使い方が非常に簡単なのでRNA-seqのQCでは非常によく使用されています！
FastQCはGUIでもCUIでもどちらでも操作可能です！
数ファイル程度ならGUIで処理した方が早いと思いますが、今回はペアエンドの8サンプル分で、16ファイルの確認を行う必要があるのでCUIで処理します！
以下のように非常に簡単な処理で品質チェックを行うことができます！
# 全てのfastq.gzファイルをアレイに格納する
fastq_files=($(ls data/seq/*.fastq.gz))
# 出力ディレクトリを指定して、QCを実行
fastqc -o ./data/reports_fastqc ${fastq_files[@]}今回のデータだと、1ファイル当たり2分程度で処理されて、全体で35分程度かかりました！
これを実行すると./data/reports_fastqcディレクトリにSRR1039508_1_fastqc.htmlなどのファイルが作成されます！
このhtmlファイルがFastQCの出力結果です！ブラウザなどで結果を閲覧することができます！
早速中身を見ていきましょう！
FastQCの結果の見方はこちらのサイト（https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/）に詳細が記載されています！
ここでは簡単にだけ見方をご紹介しておきます！
FastQCの結果に11の評価項目があり、それぞれに以下の3段階の評価がなされます！
: OK
:Warning
: Failure
全ての項目で緑のチェックマークのOKになるのが好ましいですが、目的によってはいくつかFailureが混じっていても問題ないと言及されています！
画像でお示しした、Per base sequence qualityは横軸がリードの何番目の塩基かを示し、縦軸がその塩基のクオリティ（その塩基のベースコール精度）を示しています！
塩基のクオリティはFASTQファイルのクオリティスコアから算出されたもので、Phred quality scoreシステムという表現方法がとられています！スコアの最大は40で最小は0です。
クオリティスコアが40のとき、その塩基が正しい確率が99.99%となります！
FASTQファイルについてはこちらの記事でご紹介しています！
はじめに RNA-seq解析では本当にたくさんのファイル形式を扱うことになります！ どのファイルフォーマットがどんなデータを保存していて、どんな形式で記載されていて、この列は何を意味しているのかを把握するのは本当に大変で、私も忘れる[…]
今回のサンプルではFailureは含まれておらず、WarningとなっているSequence Duplication LevelsなどもRNA-seq解析を目的とするサンプルであればWarningになるのはとてもよく見られるため気にする必要のない項目ですのでこのまま解析を進めていきます！
RNA-seq解析 ~Phase2~ 参照ゲノムへのマッピング(Reference genome mapping)
生シーケンスリードのクオリティが確認できたら、次は参照ゲノムへマッピングを行います！
参照ゲノムへマッピングを行うソフトウェアはいくつも存在していて、それぞれに良さがあるのですが、今回はSTARというソフトウェアを使用します！
STARの概要
まずはSTARの概要についてご説明します
Spliced Transcripts Alignment to a Reference (STAR)の参考はこちらです。
Dobin, A., Davis, C.A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, P., Chaisson, M., Gingeras, T.R. (2013) STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29, 15–21.
https://academic.oup.com/bioinformatics/article/29/1/15/272537
STARは80 billion reads以上の巨大なRNA-SeqデータであるENCODE Transcriptome RNA-seq datasetを参照ゲノムへマッピングするために開発されたソフトウェアです。
STARは非常に高速に参照ゲノムへのマッピングを行うことができ、またその正確性も高いです！
- 参照ゲノムのダウンロード
- Indexファイルの作成
- 参照ゲノムへのマッピング
参照ゲノムのダウンロード
# ftpへ接続
lftp ftp.ensembl.org/pub/
# 目的のフォルダに移動し、ファイルをダウンロード
cd /pub/release-108/fasta/homo_sapiens/dna
get Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
cd /pub/release-108/gtf/homo_sapien/
get Homo_sapiens.GRCh38.108.gtf.gz
# ftpから離脱
exit
# ダウンロードしたファイルは圧縮されているので解凍する
gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip Homo_sapiens.GRCh38.108.gtf.gz
Homo_sapiens.GRCh38.dna.primary_assembly.fa.gzというファイルとHomo_sapiens.GRCh38.108.gtf.gz というファイルをダウンロードしました！
これら二つのファイルを参照ゲノムとして使用します！/pub/release-108 の108はリリースされたバージョンを示しているので、最新版のバージョンを確認して使用してください！
さて！それではシーケンスリードを参照ゲノムへマッピングしていきましょう！
STARによるIndexファイルの作成
STARによる参照ゲノムへのマッピングを行うためにはIndexファイルというファイルを作成する必要があります！
IndexファイルはSTARが内部的に使用するためのものです！
以下のコマンドでSTARのIndexファイルの作成を行います！
このコマンドの実行には20GB以上というかなりのメモリが必要なので、ノートPCでは実行できないかもしれないです。その際はAWSなどの利用をご検討ください！
# Indexファイルを作成
STAR --runMode genomeGenerate \
--genomeDir data/ref/STAR_reference/ \
--genomeFastaFiles data/ref/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
--sjdbGTFfile data/ref/Homo_sapiens.GRCh38.108.gtf
簡単に引数について説明しておきます！
--runMode genomeGenerate これはIndexファイル（参照ゲノム）を作成するモードで実行することを支持するものです！
このコマンドでdata/ref/STAR_reference/ にIndexファイルが作成されました！
STARによるreference genomeへのマッピング
先ほど作成したIndexファイルを使用してリードを参照ゲノムへマッピングしていきます！
コマンドは以下です！このコマンドの実行には最低でも32GBのメモリが必要になります！またかなりの時間がかかり、私の環境では一晩かかりました！
# output dirを設定
output_dir=result/BAM/
srrids=(SRR1039509 SRR1039513 SRR1039517 SRR1039521 SRR1039508 SRR1039512 SRR1039516 SRR1039520)
for srrid in ${srrids[@]}; do
STAR --runMode alignReads
--quantMode TranscriptomeSAM \
--genomeDir data/ref/STAR_reference
--readFilesCommand gunzip -c \
--readFilesIn data/seq/${srrid}_1.fastq.gz data/seq/${srrid}_2.fastq.gz \
--outSAMtype BAM SortedByCoordinate \
--outFileNamePrefix ${output_dir}${srrid} \
--runThreadN 1
done 簡単に引数の説明をしておきます！
--runMode alignReads これはシーケンスリードをマッピングするモードで実行することを指示するものです！上のIndexファイルの作成では--runMode genomeGenerate を指定しましたね！
特に重要な引数は--quantMode TranscriptomeSAM です！これはリードをGenomeではなくて、Transcriptomeにマップすることを指定するためのものです！
なぜ重要なのかというと、この後RSEMというソフトウェアを使用して定量を行うのですが、RSEMはシーケンスリードがTranscriptomeにマップされていることを前提にしたアルゴリズムで動いているからです！
これが終了すれば、参照ゲノムへのマッピングは終了です！
SRR1039509Aligned.toTranscriptome.out.bamというようなファイルが出力されます！
ファイル名の通り、シーケンスリードをTranscriptomeにマップした結果で、BAM形式で保存されています！
RNA-seq解析 ~Phase3~ 遺伝子発現の定量化(Quantification)
続いて、RSEMというソフトウェアを用いてシーケンスリードが各遺伝子にどれだけマッピングされたのかを定量していきます！
RSEMもSTARと同様にIndexファイルを作成する必要があります！以下のコードで作成できます！
rsem-prepare-reference --num-threads 4
--gtf data/ref/Homo_sapiens.GRCh38.108.gtf HOomo_sapiens.GRCh38.dna.primary_assembly.fa \
data/ref/RSEM_reference Indexファイルが作成できたら、遺伝子発現の定量を行なっていきます！
# 入力ファイルが含まれているディレクトリを指定
input_dir=result/BAM/
srrids=(SRR1039509 SRR1039513 SRR1039517 SRR1039521 SRR1039508 SRR1039512 SRR1039516 SRR1039520)
for srrid in ${srrids[@]}; do
rsem-calculate-expression --num-threeds 4 --paired-end
--bam ${input_dir}${srrid}Aligned.toTranscriptome.out.bam \
data/ref/RSEM_reference ${srrid}
done
STARの出力ファイルである${srrid}Aligned.toTranscriptome.out.bam を入力にして、定量を行います！
このコマンドの結果として、${srrid}.gemes.resultsと${srrid}.isoforms.resultsが出力されます！
.resultsファイルはテキスト形式のファイルで、7MBくらいの容量でそれほど大きくないのでエクセルなどでも開くことができます！
実際に中身を見てみましょう！
.resultsファイルの中には遺伝子のEnsenble gene IDとそれにいくつのリードがマップされているのかが記載されていますね！
これを使用して、次回！発現変動遺伝子解析（DEG解析）を行なっていきます！
参考
https://bioinformatics-core-shared-training.github.io/RNAseq-R/