Hi everyone,
I am trying to identify differential gene expression in 64 Human Heart samples (21 normal & 43 heart-failure). I used GSNAP for aligning my reads to the reference genome (filtered, sorted & indexed the bam files):
gsnap -d <genome_file> --gunzip --batch=5 --nthreads=2 --distant-splice-penalty=10000 --clip-overlap -s <splice_file> -N 1 --npaths=100 -Q -v <snp_file> --read-group-id=sample1 --read-group-name=sample1 --read-group-library=sample1 --read-group-platform='Illumina' --format=sam sample1_R1.fastq.gz sample1_R2.fastq.gz 2> sample1.out | samtools view -bS - > sample1.bam
Following this, I filtered, sorted & indexed my bam files for all 64 samples. Then I ran cufflinks on the resulting bam files:
cufflinks -p 2 -o sample1-cufflinks-output sample1_filter_sort.bam 2> sample1.out
Then, I used cuffmerge to get the merged.gtf file:
cuffmerge -g <GTF-file> -s hg19.fa -p 25 <assembly_list.txt>
Following this, I used cuffdiff to get differentially expressed genes using:
cuffdiff -o heart_cuffdiff -L Control,HF --multi-read-correct --upper-quartile-norm -b hg19.fa -p 20 ./cuffmerge/merged.gtf <comma-separated-Normals> <comma-separated-HeartFailure>
Everything ran perfectly and I am interested in the following files in the cuffdiff output directory:
$ls *.diff
cds.diff cds_exp.diff gene_exp.diff isoform_exp.diff promoters.diff splicing.diff tss_group_exp.diff
$ls *.fpkm_tracking
cds.fpkm_tracking genes.fpkm_tracking isoforms.fpkm_tracking tss_groups.fpkm_tracking
The *.diff files will tell me how many genes are found to be differentially expressed between normals & heartfailing samples. And the *.fpkm_tracking files will give me metrics such as coverage and length etc. But I found no coverage information in the fpkm_tracking files. In fact, the class_code, nearest_ref_id, length, and coverage fields all have '-' value. All the *.fpkm_tracking files look like:
$head genes.fpkm_tracking
tracking_id class_code nearest_ref_id gene_id gene_short_name tss_id locus length coverage Control_FPKM Control_conf_lo Control_conf_hi Control_status HF_FPKM HF_conf_lo HF_conf_hi HF_status
XLOC_000001 - - XLOC_000001 - TSS1 chr1:12658-29348 - - 0.747812 0 2.09809 OK 0.900845 0 2.3693 OK
How can I get the these missing field values especially the coverage information? Just an FYI, in the cuffdiff commandline, I am supplying the merged GTF output from cuffmerge and not the cuffcompare output.
