Cufflinks

Cufflinks は RNA-seq を利用した遺伝子発現解析を行うためのプログラムをパッケージしている。例えば、マッピング結果からアセンブリーを行う cufflinks やマッピング結果を利用して発現変動遺伝子を検出する cuffdiff などのプログラムが含まれている。Cufflinks は TopHat と一緒に利用されることが多い。

Cufflinks

Cufflinks は TopHat のマッピング結果からアセンブリーを行うプログラムである。アセンブリーの結果は GTF 形式のファイルとして保存される。例えば、現在ディレクトリ中に SRR1976498, SRR1976499, SRR1976500, SRR1976501 のマッピング結果があるとする(TopHat2 のよるマッピング)。これらの TopHat のマッピング結果に対して Cufflinks を利用してアセンブリーを行う。

cufflinks -o ./SRR1976498 ./SRR1976498/accepted_hits.bam
cufflinks -o ./SRR1976499 ./SRR1976499/accepted_hits.bam
cufflinks -o ./SRR1976500 ./SRR1976500/accepted_hits.bam
cufflinks -o ./SRR1976501 ./SRR1976501/accepted_hits.bam

アノテーションファイルが存在するとき、以下のように、アノテーションファイルをアセンブリーのガイドとして使うことも可能である。

cufflinks -g Arabidopsis_thaliana.TAIR10.26.gff3 -o ./SRR1976498 ./SRR1976498/accepted_hits.bam
cufflinks -g Arabidopsis_thaliana.TAIR10.26.gff3 -o ./SRR1976499 ./SRR1976499/accepted_hits.bam
cufflinks -g Arabidopsis_thaliana.TAIR10.26.gff3 -o ./SRR1976500 ./SRR1976500/accepted_hits.bam
cufflinks -g Arabidopsis_thaliana.TAIR10.26.gff3 -o ./SRR1976501 ./SRR1976501/accepted_hits.bam

cufflinks コマンドを 4 つの BAM ファイルに対して実行したため、アセンブリーの結果も 4 つのファイルに保存される。そこで、cuffmerge コマンドを利用してこれら 4 ファイルをマージする。cuffmerge を実行するには、マージするファイルを指定する必要がある。この指定はテキストファイルに記述する必要がある。そこで、上述 4 つのアセンブリ結果ファイルの場所を merge_info.txt という名前のファイルに保存する。

echo -e "./SRR1976498/transcripts.gtf
./SRR1976499/transcripts.gtf
./SRR1976500/transcripts.gtf
./SRR1976501/transcripts.gtf" > merge_info.txt

merge_info.txt のファイルを準備出来たら、次に cuffmerge を実行する。-o オプションで出力ファイル名のプレフィックスを指定する。

cuffmerge -o all merge_info.txt

cuffmerge の実行が終えると、all ディレクトリが作られ、その中に merged.gtf ファイルが生成される。merged.gtf が上述 4 つのアセンブリー結果をマージしたアノテーションファイルである。遺伝子の発現量などを取得するときにこの merged.gtf ファイルを利用する。

Cuffdiff

Cuffdiff はマッピング結果から発現変動遺伝子を検出する際にりようするプログラムである。入力ファイルとしてマッピング結果の BAM ファイル、遺伝子のアノテーションファイル(GTF ファイル)を必要とする。GTF ファイルは、Ensembl などのデータベースからダウンロードしたファイルを利用してもよいし、Cufflinks を利用して作成したファイルを用いてもよい。

ここで SRR1976498, SRR1976499, SRR1976500, SRR1976501 のマッピング結果を利用して発現変動遺伝子を検出してみる。このサンプルデータは、SRR1976498 と SRR1976499 が野生型、SRR1976500 と SRR1976501 が変異型のデータである。野生型と変異型を比較する。同一群のデータはカンマ区切りで指定し、この際に空白(スペース)を入れないようにすること。

Ensembl などでダウンロードしたアノテーションファイルを利用する場合は以下のようにする。解析結果は de_result_1 ディレクトリに保存される。

cuffdiff -o ./de_result_1 Arabidopsis_thaliana.TAIR10.26.gff3 \
         ./SRR1976498/accepted_hits.bam,./SRR1976499/accepted_hits.bam \
         ./SRR1976500/accepted_hits.bam,./SRR1976501/accepted_hits.bam

例えば、上の例で自作した merged.gtf を利用する場合は、merged.gtf は all ディレクトリの下にあるので、以下のように実行する。解析結果は de_result_2 ディレクトリに保存される。

cuffdiff -o ./de_result_2 ./all/merged.gtf \
         ./SRR1976498/accepted_hits.bam,./SRR1976499/accepted_hits.bam \
         ./SRR1976500/accepted_hits.bam,./SRR1976501/accepted_hits.bam

解析結果のディレクトリ(de_result_1 または de_result_2)には複数の結果ファイルが保存される。

FPKM tracking files

マッピング結果から計算された、各サンプルにおける各遺伝子の PFKM(genes.fpkm_tracking)、各トランスクリプトの FPKM(isoforms.fpkm_tracking)などが保存されている。このページで取り上げた解析例では、サンプル数は野生型と変異型の 2 サンプルである。
出力ファイルはタブ区切りで、q1_FPKM と q2_FPKM の列にそれぞれのサンプルの FPKM が記載されている。

head genes.fpkm_tracking 
## tracking_id	class_code	nearest_ref_id	gene_id	gene_short_name	tss_id	locus	length	coverage	q1_FPKM	q1_conf_lo	q1_conf_hi	q1_status	q2_FPKM	q2_conf_lo	q2_conf_hi	q2_status
## 	-	-	-	-	-	1:3630-5899	-	-	123847	0	2.5738	OK	145017	0	0.444586	OK
## gene:AT1G01010	-	-	gene:AT1G01010	-	-	1:3630-5899	-	-	2.57921	1.34594	3.81249	OK	0.852512	0.315646	1.38938	OK
## gene:AT1G01020	-	-	gene:AT1G01020	-	-	1:5927-8737	-	-	12.1259	7.8105	16.4413	OK	14.2707	9.33157	19.2098	OK
## gene:AT1G01030	-	-	gene:AT1G01030	-	-	1:11648-13714	-	-	1.20879	0.52826	1.88932	OK	1.77622	0.879587	2.67285	OK
## gene:AT1G01040	-	-	gene:AT1G01040	-	-	1:23145-33153	-	-	11.1449	6.18508	16.1046	OK	11.4139	6.27133	16.5565	OK
## gene:AT1G01046	-	-	gene:AT1G01046	-	-	1:23145-33153	-	-	54.7357	0	203.58	OK	51.1632	0	198.065	OK
## gene:AT1G01050	-	-	gene:AT1G01050	-	-	1:23145-33153	-	-	86.5093	48.4253	124.593	OK	89.9796	50.6423	129.317	OK
## gene:AT1G01060	-	-	gene:AT1G01060	-	-	1:33378-37871	-	-	432.409	291.176	573.642	OK	488.693	329.853	647.534	OK
## gene:AT1G01070	-	-	gene:AT1G01070	-	-	1:38751-40944	-	-	4.92294	2.69579	7.15009	OK	6.0619	3.57497	8.54884	OK
Count tracking files

マッピング結果から推定された、各サンプルにおける各遺伝子のリードカウント(genes.count_tracking)、各トランスクリプトの リードカウント(isoforms.count_tracking)などが保存されている。paired-end リードの場合はフラグメントのカウントとなる(つまり 2 read が 1 count である)。このページで取り上げた解析例では、サンプル数は野生型と変異型の 2 サンプルである。
出力ファイルはタブ区切りで、q1_count と q2_count の列にそれぞれのサンプルの FPKM が記載されている。

head genes.count_tracking 
## tracking_id	q1_count	q1_count_variance	q1_count_uncertainty_var	q1_count_dispersion_var	q1_status	q2_count	q2_count_variance	q2_count_uncertainty_var	q2_count_dispersion_var	q2_status
## 	426378	43	0	42.1016	OK	503047	2	0	1.4031	OK
## gene:AT1G01010	138.574	1097.62	0	1067.59	OK	45.8032	208	0	207.484	OK
## gene:AT1G01020	386.771	4737.32	0	4735.41	OK	456.106	6224.03	0	6222.41	OK
## gene:AT1G01030	74.4418	439.104	0	437.215	OK	109.386	762.25	0	754.612	OK
## gene:AT1G01040	2327.79	268592	0	267616	OK	2435.43	299247	0	288549	OK
## gene:AT1G01046	62.6147	7248	0	7247.06	OK	58.5279	7060	0	7059.46	OK
## gene:AT1G01050	2417.94	283262	0	277305	OK	2514.93	302213	0	296482	OK
## gene:AT1G01060	37473.7	3.76327e+07	0	3.71215e+07	OK	42367.2	4.78193e+07	0	4.77579e+07	OK
## gene:AT1G01070	197.304	1992.05	0	1990.16	OK	243.311	2487.57	0	2447.96	OK
Read group tracking files

マッピング結果から推定された、各ライブラリー(または replicate)における各遺伝子のリードカウント(genes.read_group_tracking)、各トランスクリプトの リードカウント(isoforms.read_group_tracking)などが保存されている。paired-end リードの場合はフラグメントのカウントとなる。このページで取り上げた解析例では、2 サンプル 4 ライブラリーである。

head genes.read_group_tracking 
## tracking_id	condition	replicate	raw_frags	internal_scaled_frags	external_scaled_frags	FPKM	effective_length	status
## 	q1	0	0	0	0	0	-	OK
## 	q1	1	7.66182	9.94303	9.94303	1.41495	-	OK
## 	q2	0	0.705594	0.589297	0.589297	0.0838604	-	OK
## 	q2	1	0.00213391	0.00215615	0.00215615	0.000306833	-	OK
## gene:AT1G01010	q1	0	110	98.9196	98.9196	1.84114	-	OK
## gene:AT1G01010	q1	1	137.338	178.229	178.229	3.31728	-	OK
## gene:AT1G01010	q2	0	61.2944	51.1918	51.1918	0.952807	-	OK
## gene:AT1G01010	q2	1	39.9979	40.4147	40.4147	0.752218	-	OK
## gene:AT1G01020	q1	0	414.272	372.542	372.542	11.6798	-	OK
Differential expression tests

Cuffdiff による発現変動遺伝子の同定結果が保存されている。遺伝子レベルでの同定結果は gene_exp.diff、トランスクリプトレベルの同定結果は isoform_exp.diff に保存されている。
解析結果はタブ区切りのテキスト形式で保存され、2 列目に遺伝子またはトランスクリプトの名前、7 列目にサンプル 1 の FPKM、8 列目にサンプル 2 の FPKM、10 列目に統計量、11 列目に p-value、12 列目に FDR が記載されている。

head gene_exp.diff 
## test_id	gene_id	gene	locus	sample_1	sample_2	status	value_1	value_2	log2(fold_change)	test_stat	p_value	q_value	significant
## 	-	-	Mt:296819-298204	q1	q2	OK	0.707476	0.0420836	-4.07135	-0.568876	0.42	0.997966	no
## gene:AT1G01010	gene:AT1G01010	-	1:3630-5899	q1	q2	OK	2.57921	0.852512	-1.59714	-2.80016	0.0003	0.0096132	yes
## gene:AT1G01020	gene:AT1G01020	-	1:5927-8737	q1	q2	OK	12.1259	14.2707	0.234962	0.656142	0.35215	0.993575	no
## gene:AT1G01030	gene:AT1G01030	-	1:11648-13714	q1	q2	OK	1.20879	1.77622	0.555242	1.01795	0.1487	0.823331	no
## gene:AT1G01040	gene:AT1G01040	-	1:23145-33153	q1	q2	OK	11.1449	11.4139	0.0344158	0.0753381	0.91485	0.997966	no
## gene:AT1G01046	gene:AT1G01046	-	1:23145-33153	q1	q2	OK	54.7357	51.1632	-0.0973762	-0.0341355	0.89905	0.997966	no
## gene:AT1G01050	gene:AT1G01050	-	1:23145-33153	q1	q2	OK	86.5093	89.9796	0.0567425	0.126787	0.8567	0.997966	no
## gene:AT1G01060	gene:AT1G01060	-	1:33378-37871	q1	q2	OK	432.409	488.693	0.176534	0.531108	0.44575	0.997966	no
## gene:AT1G01070	gene:AT1G01070	-	1:38751-40944	q1	q2	OK	4.92294	6.0619	0.300251	0.68155	0.3347	0.983331	no
Differential splicing tests, Differential coding output, Differential promoter use

その他の出力ファイル。

発現変動遺伝子

各群の遺伝子発現量は FPKM またはカウントデータとしてファイルに保存されている。また、発現変動遺伝子の同定結果もファイルに保存されている。これらのファイルからデータを取り出して、プロットすることができる。

fpkm <- read.table("genes.fpkm_tracking", sep = "\t", header = T)
stat <- read.table("gene_exp.diff", sep = "\t", header = T)

x <- fpkm[, "q1_FPKM"]
y <- fpkm[, "q2_FPKM"]
cols <- ifelse(stat[, "q_value"] < 0.01, "orange", "grey")

plot(log2(x), log2(y), col = cols, pch = 19, cex = 0.5)
Cuffdiffによる発現変動遺伝子の同定結果

References

  1. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012, 7(3):562-78. PubMed Abstract