選択的スプライシングにより 1 つの遺伝子から複数種類の転写産物が生成される場合がある。そのため、マッピング結果から遺伝子発現量を見積もるときに、実験目的に応じて転写産物の isoform を無視して遺伝子ごとに計上する場合と、isoform を考慮して isoform 毎に計上する場合がある。前者の場合は、featureCounts が一般的に使われている。これに対して、後者の場合は Cufflinks や StringTie が一般的に使われている。このページでは、Cufflinks について述べるが、StringTie の方が比較的新しく開発されたプログラムであり、特別な理由がなければ Cufflinks よりも StringTie を使用することが推奨される。

Cufflinks を使用した発現量定量の流れとして、(1) BAM ファイルから新規 isoform と思われるもののアノテーションを作成し、(2) 既存のアノテーションと新しいアノテーションをマージしてから発現量の定量を行う、の 2 ステップからなる。ここで計算される遺伝子発現量は FPKM として出力される。この発現量に対して引き続き発現変動遺伝子を検出したい場合は Cuffdiff を使用する。

Cufflinks

Cufflinks を利用することで、TopHat などのマッピング結果結果を利用して、トランスクリプトのアセンブリを行い、転写産物 isoform 毎の発現量を推定できるようになる。トランスクリプトのアセンブリは、Cufflinks に BAM ファイルを与えることで、そのアセンブリの結果が GTF ファイルで出力される。この処理を個々のサンプルに対して行う必要がある。例えば、サンプルが 4 つあれば、次のように 4 回 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

なお、モデル生物のデータを解析するとき、すでにアノテーションが存在する場合がある。この場合、既存のアノテーションをガイドとして使うことこともできる。例えば、モデル生物シロイヌナズナの場合は、次のように既存のアノテーション gff3 を追加で与えることができる。

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 により、サンプル毎のアノテーションが作られる。次に、Cuffmerge を使用して、これらのアノテーションをマージしていく。まず、上で作った 4 つのアノテーションファイル（.gtf）へのパスをテキストファイルに記載し、次にこのテキストファイルを Cuffmerge に与えてマージを行う。

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

cuffmerge -o all merge_info.txt

プログラムが正しく実行されると、all ディレクトリが作られ、その中に merge.gtf が作られる。この merge.gtf を使って必要に応じて Cuffidff を使用して発現変動遺伝子を同定したりすることができる。

Cuffdiff

Cuffdiff は Cufflinks により定量した遺伝子発現量（FPKM）を利用して発現変動遺伝子を検出するプログラムである。入力ファイルとして、マッピング結果（BAM ファイル）と遺伝子のアノテーションファイル（GTF ファイル）である。GTF は Cuffdiff でマージしたファイルを用いる。

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

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

解析結果のディレクトリには複数の結果ファイルが出力される。

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

発現変動遺伝子

各群の遺伝子発現量は 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)