シークエンス解析実習2
〜第7日目：RNA-seq解析〜 - 担当：鈴木研究室
RNA-seq解析 illumina vs MinION
※ group2 のようなファイル名が出てきますが、僕が2班だったというだけです。
謝辞
今回の実験では、東京大学 医科学研究所 ヒトゲノム解析センターにスーパーコンピュータをお借りしました。この場を借りて御礼申し上げます。
# スパコンにログイン
$ ssh XXX@slogin.hgc.jp
# スパコンからファイルをコピー
$ scp XXX@slogin.hgc.jp:/path/to/file .原理
illumina
以下の流れでシークエンスを行なっています。
- ターゲットとなるDNAを抽出後、断片化
- 両端に2種類のアダプターを付与
- ブリッジPCR（ブリッジ結合→伸長→変性→ブリッジ結合→…）による増幅
- 蛍光標識したdNTPの取り込みを蛍光顕微鏡により解析
以下の動画（引用：illumina） が参考になります。
MinION
MinIONはロングリードをダイレクトに読むシークエンサーであり、イルミナが100~200 bの配列を扱うのに対し、10 kb以上の配列を扱うことができます。しかし、精度が \(90\) %程度と低いのが難点です。
原理としては、ナノポアと呼ばれるナノスケールの穴の中や、ナノポア付近を生体分子が通過するときに発生する電流の変化を計測することで、その生体分子を特定するという仕組みです。以下の動画（引用：Oxford Nanopore Technologies） が参考になります。
MinIONの特徴として、以下の両方が行えることが挙げられます。
- cDNA sequencing（今回はこっち） RNAを逆転写してcDNAにしてから増幅してシークエンスする。少ないサンプルでシークエンスが可能。
- Direct RNA sequencing RNAをダイレクトにシークエンスする。多くのインプットが必要だが、増幅しないので、RNA修飾が解析可能。
なお、分子バーコーディングによる再構築（長いDNAにバーコードを入れてからショートリードにして読み、バーコードを手がかりに再構築する）というアイデアでロングリードを読むシークエンサーもあります。
解析の流れ
illumina
以下の流れに沿って解析を行います。
graph LR;
A((BCL))-->B;
B((fastq))-->C;
C[マッピング]-->D;
D((BAM))-->F[リードカウント];※ 6/28,7/2 の二日間でウェットの実験を行い、サンプルをアプライし、シークエンスを行いました。
この時、生データファイルはバイナリベースコール(BCL)形式で生成されるようで、それをbcl2fastqというソフトウェアを使ってFASTQ形式というNGSデータを保管する標準形式（生シーケンスデータとクオリティスコアを保存）に変換します。（これは、研究室の方々にやっていただきました。）
続いて、これを参照ゲノム配列（今回用いたのはUCSCのhg38）にSTARというライブラリを用いてマッピングします。
マッピングした結果は、SAM形式のファイルになります。このファイルはFASTQの中身も含み、数十ギガバイトの巨大なファイルとなるので、一般にBAM形式に圧縮して提供されます。
この時特殊な圧縮のされ方をしてるため、専用のツール（samtools）を用いて展開する必要がありますが、通常はIGVなどのゲノムブラウザにそのまま読み込ませて結果を確認します。
最後に、発現量を算出し、illuminaとMinIONの結果を比べます。これは featureCounts という、Subreadパッケージ内の、BAMファイルを入力として各feature領域（gene,CDS,exon）ごとのリード数をカウントするプログラムを利用します。
MinION
以下の流れに沿って解析を行います。
illuminaと異なり、MinIONで得られるのはイオンシグナル情報なので、それを塩基配列の情報に変換するのに Guppy というツールを用いてベースコール情報を付加する必要があるのと、マッピング時に長いリードや高いエラー率に対応したツールを用いる必要がある点のみが異なります。
graph LR;
A((fast5))-->G;
G[ベースコール]-->B;
B((fastq))-->C;
C[マッピング]-->D;
D((BAM))-->F[リードカウント];※ 6/26,6/27 の二日間でウェットの実験を行い、サンプルをアプライし、シークエンスを行いました。
マッピング
illuminaとMinIONでは扱うツールが異なりますが、ソフトウェア化されているので、ユーザー側からすれば大差がありません。
何をやっているのかを理解したいので、できるだけいただいたスクリプトの中身を丁寧に見たいと思います。
STAR(illumina)
STARを用いて、ヒトリファレンスゲノム配列へのマッピングを行います。
なお、この際東大医科研HGCスパコンは複数のユーザーで利用しているため、キューイングシステムを利用しています。
キューイングシステムは、実行したいプログラムを必要なメモリやCPUの数と共に実行ファイル（基本的にはシェルスクリプト）に書き込んでキューイングシステムにサブミットし、キューイングシステムが総合判断して複数ユーザーが投入したジョブを円滑に実行する、という仕組みです。
今回用意していただいたスクリプトは以下のようになっていました。
#!/bin/bash
#$ -S /bin/bash
#$ -l s_vmem=35G,mem_req=35G
#$ -l os7
source /home/lect01/2019_jisyu.source # ツールの場所を指定
STAR --genomeDir /home/lect01/reference/STAR_hg38 \ # リファレンスゲノムが置いてあるディレクトリ
--readFilesIn 2_Read1.fastq 2_Read2.fastq \ # 入力fastqファイル
--chimSegmentMin 30 \ # リードの前半と後半で異なる染色体にマッピングされるキメラ(fusion)リードを検出するために指定。長さの短い方がこの値よりも長い時、そのリードをキメラリードとみなす。
--chimOutType WithinBAM \ # キメラリードのマッピング結果を通常のマッピング結果と同じファイルに保存する。（別ファイルの場合はSeparateSAMold）
--outSAMtype BAM SortedByCoordinate \ # 遺伝子座でソートされた、Aligned.sortedByCoord.out.bam というファイルで出力する
--sjdbGTFfile /home/lect01/reference/STAR_hg38/genes.gtf \ # 遺伝子のアノテーション情報
--outFilterMultimapNmax 1 \ # 1つのリードがマッピングできる位置の最大値。Uniquely mapped readsのみを解析対象としているので、今回は1
--outFileNamePrefix group2_ # ファイルの接尾辞を追加
exit2019_jisyu.source ファイル
#-----------------------------------------------------------
# TOOLS
#-----------------------------------------------------------
TOOLS=/home/lect01/tools/bin
#-----------------------------------------------------------
# PATH
#-----------------------------------------------------------
export PATH=/PATH/to/where/STAR/is:${TOOLS}:${PATH}
export PYTHONHOME=/path/to/Python/home # ex.) /usr/local/package/python/version/
export PATH=${PYTHONHOME}/bin:${PATH}
export LD_LIBRARY_PATH=:${PYTHONHOME}/lib:${LD_LIBRARY_PATH}なお、group2_Aligned.sortedByCoord.out.bam というファイルができたら、インデックス付けをします。
$ samtools index group2_Aligned.sortedByCoord.out.bamこれで、group2_Aligned.sortedByCoord.out.bam.bai というファイルが作成されます。
Minimap2(MinION)
ロングリードシークエンスデータをマッピングするツールは様々なものが開発されていますが、今回はMinimap2を利用しました。
| 名前 | 特徴 |
|---|---|
| LAST | 我らがMartin Frith先生 |
| Minimap2 | Direct RNA-seq用のパラメータが存在する。 |
| Minialign | 速くて正確 |
| GraphMap |
今回用意していただいたスクリプトは以下のようになっていました。
#!/bin/bash
#$ -S /bin/bash
#$ -l s_vmem=20G,mem_req=20G
#$ -l os7
source /home/lect01/2019_jisyu.source # ツールの場所を指定
minimap2 -ax splice --secondary=no /home/lect01/reference/minimap2/genome.mmi 2.pass.fastq > pass.sam
# /home/lect01/reference/minimap2/genome.mmi: リファレンスゲノム
# 2.pass.fastq: 入力fastqファイル
exit同様に、pass_sort.bam というファイルができたら、インデックス付けをします。なお、Minimap2では染色体順にソートがされていないので、ソートも行います。
$ samtools view -bS pass.sam | samtools sort - pass_sort
$ samtools index pass_sort.bamこれで、pass_sort.bam.bai というファイルが作成されます。
データの可視化
データの可視化には、Integrative Genomics Viewer(IGV)というツールを用います。IGVにBAMファイルを選択するだけで結果が可視化されます。
今回は、「ヒト肺腺癌の培養細胞のうち各班に割り当てられたのが何だったのか？」を当てる実習であるので、「各細胞で既知の変異が、自分の細胞に起こっているか？」を調べることで推測を行いました。
なお、既知の変異が染色体上のどこで起こっているのかはCatalogue Of Somatic Mutations In Cancer(COSMIC)を参照しました。
発現量の算出と比較
異なるプラットフォームの遺伝子発現量における相関を解析します。なお、使用する遺伝子についてだけでも様々な手法があります。
- 全遺伝子を使う
- タンパク質コード遺伝子だけを使う
- どちらかで5以上発現している遺伝子を使う
正規化
リード数がデータによって異なるので、総リード数を用いて各遺伝子にマップされたリード数を正規化します。
- ppm: parts per million reads
なお、illuminaショートリードの場合、遺伝子領域が長いほどより多くのリードがマップされることが知られているので、遺伝子の長さでも補正する必要があります。
- rpkm: reads per kilobase per million reads
illumina(rpkm)
$ featureCounts -p -B -t exon -g gene_id \
-a /home/lect01/reference/genes.gtf \
-o illumina_counts.txt group2_Aligned.sortedByCoord.out.bam ========== _____ _ _ ____ _____ ______ _____
===== / ____| | | | _ \| __ \| ____| /\ | __ \
===== | (___ | | | | |_) | |__) | |__ / \ | | | |
==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
v1.6.2
//========================== featureCounts setting ===========================\\
|| ||
|| Input files : 1 BAM file ||
|| P group2_Aligned.sortedByCoord.out.bam ||
|| ||
|| Output file : illumina_counts.txt ||
|| Summary : illumina_counts.txt.summary ||
|| Annotation : genes.gtf (GTF) ||
|| Dir for temp files : ./ ||
|| ||
|| Threads : 1 ||
|| Level : meta-feature level ||
|| Paired-end : yes ||
|| Multimapping reads : not counted ||
|| Multi-overlapping reads : not counted ||
|| Min overlapping bases : 1 ||
|| ||
|| Chimeric reads : counted ||
|| Both ends mapped : required ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//
//================================= Running ==================================\\
|| ||
|| Load annotation file genes.gtf ... ||
|| Features : 535393 ||
|| Meta-features : 26485 ||
|| Chromosomes/contigs : 255 ||
|| ||
|| Process BAM file group2_Aligned.sortedByCoord.out.bam... ||
|| Paired-end reads are included. ||
|| Assign fragments (read pairs) to features... ||
|| ||
|| WARNING: reads from the same pair were found not adjacent to each ||
|| other in the input (due to read sorting by location or ||
|| reporting of multi-mapping read pairs). ||
|| ||
|| Pairing up the read pairs. ||
|| ||
|| Total fragments : 9265506 ||
|| Successfully assigned fragments : 6667566 (72.0%) ||
|| Running time : 0.55 minutes ||
|| ||
|| Read assignment finished. ||
|| ||
|| Summary of counting results can be found in file "illumina_counts.txt.sum ||
|| mary" ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//
MinION(ppm)
$ featureCounts -L -t exon -g gene_id \
-a /home/lect01/reference/genes.gtf \
-o minion_counts.txt pass_sort.bam========== _____ _ _ ____ _____ ______ _____
===== / ____| | | | _ \| __ \| ____| /\ | __ \
===== | (___ | | | | |_) | |__) | |__ / \ | | | |
==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
v1.6.2
//========================== featureCounts setting ===========================\\
|| ||
|| Input files : 1 BAM file ||
|| S pass_sort.bam ||
|| ||
|| Output file : minion_counts.txt ||
|| Summary : minion_counts.txt.summary ||
|| Annotation : genes.gtf (GTF) ||
|| Dir for temp files : ./ ||
|| ||
|| Threads : 1 ||
|| Level : meta-feature level ||
|| Paired-end : no ||
|| Multimapping reads : not counted ||
|| Multi-overlapping reads : not counted ||
|| Min overlapping bases : 1 ||
|| Long read mode : yes ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//
//================================= Running ==================================\\
|| ||
|| Load annotation file genes.gtf ... ||
|| Features : 535393 ||
Status pass_sort.bam
|| Meta-features : 26485 ||
|| Chromosomes/contigs : 255 ||
|| ||
|| Process BAM file pass_sort.bam... ||
|| Single-end reads are included. ||
|| Assign reads to features... ||
|| Total reads : 515045 ||
|| Successfully assigned reads : 369980 (71.8%) ||
|| Running time : 0.28 minutes ||
|| ||
|| Read assignment finished. ||
|| ||
|| Summary of counting results can be found in file "minion_counts.txt.summa ||
|| ry" ||
|| ||
\\===================== http://subread.sourceforge.net/ ======================//
Pearson相関係数算出(R)
ILLUMINA <- read.table("illumina/illumina_counts.txt", header=T)
total_illumina <- 9265506
ppm_i <- ILLUMINA$group2_Aligned.sortedByCoord.out.bam / total_illumina * 1000000
rpkm_i <- ppm_i / ILLUMINA$Length * 1000
MINION <- read.table("minion/minion_counts.txt", header=T)
total_minion <- 515045
ppm_m <- MINION$pass_sort.bam / total_minion * 1000000
log_rpkm_i <- log2(rpkm_i + 1)
log_ppm_m <- log2(ppm_m + 1)
cor.test (log_rpkm_i, log_ppm_m, method="pearson")
# Pearson's product-moment correlation
# data: log_rpkm_i and log_ppm_m
# t = 309.92, df = 26483, p-value < 2.2e-16
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
# 0.8827357 0.8879421
# sample estimates:
# cor
# 0.8853666
pdf("comparison.pdf")
plot(log_rpkm_i, log_ppm_m, xlab="illumina, log2 (rpkm+1)", ylab="MinION, log2 (ppm+1)")
dev.off()
# null device
# 1