tophatでマッピングされたリードから、遺伝子発現量を計算させるステップがcufflinksによる計算ステップ
これによって生成されるファイルは
genes.fpkm_tracking isoforms.fpkm_tracking skipped.gtf transcripts.gtf
の4つ
genes.fpkm_trackingは遺伝子単位の発現量、isoforms.fpkm_trackingはアイソフォーム単位の発現量
でも、このあとの発現量比較にはtranscript.gtfを使う。
cuffmergeですべてのサンプルのtranscript.gtfを1つに合わせ、cuffdiffでサンプル間のFPKMを比較
cuffdiff_resultsに出力されるファイルは以下のもの
bias_params.info cds_exp.diff cds.count_tracking cds.diff cds.fpkm_tracking cds.read_group_tracking gene_exp.diff genes.count_tracking genes.fpkm_tracking genes.read_group_tracking isoform_exp.diff isoforms.count_tracking isoforms.fpkm_tracking isoforms.read_group_tracking promoters.diff read_groups.info run.info splicing.diff tss_group_exp.diff tss_groups.count_tracking tss_groups.fpkm_tracking tss_groups.read_group_tracking var_model.info
Rでデータ解析するときはcuffdiff_resultsディレクトリからデータベース(cuffData.db)を作成することになる。
cuff <- readCufflinks("hogehoge/tophat_results/cuffdiff_result")
GSEA法でエンリッチメント解析を行うときにはgene_exp.diffファイルから遺伝子を抽出することになる
2つずつのサンプル間で発現の差を遺伝子ごとに比較していっているファイルである。
比較したい遺伝子の組み合わせを5列目、6列目で見て、8,9列目がそれぞれのFPKM値、14列目に有意差判定結果が記載されている
遺伝子名がgene1, gene2の2サンプル間の比較評価から、有意に発現変動している遺伝子名を拾い出すには
awk '$5 =="gene1" && $6 == "gene2" && $14 == "yes" {print $3}' gene_exp.diff | grep -v '-' >gene_exp_diff12.txt
という感じで書き出す。