RNA-seq, in vivo study, 39 DE genes

[xulong]

Dear all,

I have an in vivo whole tissue study with RNA-seq. With the Tophat->Cufflinks workflow, I see in the dendrogram plot that the treatment replicates are grouped together, and un-treatment replicates are groups together too.
Then I did Tophat->htseq-count->DESeq workflow. This time the dendrogram plot is completely messed up, and of course there are no significantly DE genes at all.

There likely be that I did something wrong in the DESeq part. Or, because the fact that Cufflinks use FPKM while DESeq use raw count number makes the difference?

Has anyone experienced the same problem?

Thanks,

[sdriscoll]

I think whether 39 genes sounds reasonable or not should be up to you and your experiment. Number of DE genes is always variable in my experience.

Regarding the dendrogram a and pipeline differences its hard to judge without knowing every last detail of the steps you took and the commands you ran (as they were typed). If you can provide those details someone will probably be able to help you out. The difference in DE lists is not surprising since Cufflinks and htseq-count are completely different tools however the dendrogram difference sounds fishy.

[xulong]

Hi sdriscoll, appreciate your reply. You are absolutely right that it is hard to say without knowing all the commands I executed.

I am confident that the tophat output is fine, and the htseq-count procedure should be good too. Most likely is that I did not do something critical in the DESeq part, that messed up the statistics. Here is my code for DESeq.

----
# Construct the CountDataSet object for DESeq
countsTable <- read.delim("genes_merged.txt", sep="", header=F);
samplelist <- list("gene_ID","s1","s2","s3","s4","s5","s6");
colnames(countsTable) <- samplelist;
rownames(countsTable) <- countsTable$gene_ID;
countsTable <- countsTable[,-1];
cond <- c("treated","treated","untreated","untreated","treated","untreated");
cds <- newCountDataSet(countsTable,cond);

# Library size normalization
# This is anothe way
cds <- estimateSizeFactors(cds);
sizeFactors(cds);

# Dispersion estimation
cds <- estimateDispersions(cds);

# Call Differentaial Expression
res = nbinomTest(cds,"untreated","treated");

[sdriscoll]

This looks fine...the only odd thing is in your read.delim call. You've specified that the fields are separated by nothing (sep=""). Also the count test realm recommend filtering out low count features. It sounds like this is nicely explained in the new DESeq2 manual. Otherwise this looks fine.

[xulong]

Hey sdriscoll, thanks for the suggestions. I inspected the count table after the read.delim command. The table looks fine. I filtered out low reads count features too, with total of 10 reads for 5 samples as cutoff. However, I still got similar results. Distance plot messed up. Not like the nice cufflinks output.

Here is my whole code:
# DESeq pipeline

# Working directory
setwd("~/rnaseq_exp/DESeq");

# Construct the CountDataSet object for DESeq
countsTable <- read.delim("./htseq-count_union/genes_merged.txt", sep="", header=F);
samplelist <- list("gene_ID","s1","s2","s3","s4","s5");
colnames(countsTable) <- samplelist;
rownames(countsTable) <- countsTable$gene_ID;
# Take out gene id column
countsTable <- countsTable[,-1];

# Filter out genes with low cumulative read counts
countsTableFilter <- countsTable[ (rowSums(countsTable)>10), ];
cond <- factor(c("treated","treated","untreated","untreated","treated"));
cds <- newCountDataSet(countsTableFilter,cond);

# Library size normalization
cds <- estimateSizeFactors(cds);
sizeFactors(cds);

# Dispersion estimation
cds <- estimateDispersions(cds);

# Call Differentaial Expression
res = nbinomTest(cds,"untreated","treated");
sig = 0.1;

# Extract the significant genes
resSig <- res[ res$padj < sig, ];

# Order the top DE genes
resTop <- order(res$pval)[1:100];

# Distance plot
cdsblind = estimateDispersions(cds, method="blind");
vsd = getVarianceStabilizedData(cdsblind);
dists = dist( t(vsd) );

heatmap( as.matrix(dists), symm=TRUE, scale="none", margins=c(10,10));

[sdriscoll]

I don't think you're doing anything wrong. Maybe you can get Simon's attention on this...he may have suggestions to improve your DESeq results/pipeline. I know they have been working out new filtering methods which seem to improve power and maintain the FDR. I've certainly seen variance in the number of DE genes between tools. There was a time when cuffdiff reported 100's more DE genes than DESeq. Then they made it WAY more conservative and it was reporting far fewer. I haven't tested the lastest release which, as I understand it, is basically doing the same thing as DESeq. It IS concerning that not a single DE gene remained in the significant list.

Have you tried edgeR?

[xulong]

Appreciate your time sdriscoll. I am going to try edgeR tonight. As my understanding, the statistic methods that both Cufflinks and DESeq use are the same. It is pairwise t-test, followed by FDR for multiple comparison correction. What is different is DEseq takes the raw read counts, while Cufflinks use FPKM. Personally I think FPKM makes more sense because it gives normalization on both gene length and library size, while DESeq only normalize to the library size.
The nice dendrogram figure that Cufflinks gave me makes me feel confident with the Cufflinks results. But I want to make sure I did not do anything wrong in DESeq, and also want to understand more on why DESeq's results are so dramatically different.
Who is Simon? How can I reach him?

[sdriscoll]

Simon is one of the authors behind DESeq. Actually cufflinks and DESeq are even more similar than you think. Cufflinks performs its differential expression test on the counts just like DESeq and has also incorporated the exact same library normalization method that DESeq uses (at least in the current version). It reports FPKMs but it doesn't base any of the DE on them. It's been shown that you can't have reliable FPKMs unless you also have isoform level resolution of read count assignment. Their pipeline can claim to be able to do that and therefore they can report FPKMs. I like having FPKMs as well.

[xulong]

Hi Shawn, the fact that DESeq and cufflinks use the input data, same normalization, and same statistics, makes me more curious about why I got so different result with the two pipelines. I just PM Simon on this. I will put you posted from here on how it goes.

[Simon Anders]

The relevant difference between cuffdiff and DESeq is not so much the normalization or the statistic, but the variance estimation. Hence maybe post the result of plotDispEst to see what DESeq did.

Also consider trying our new DESeq2 package whose improved variance estimation method is especially beneficial for medium-sized experiments such as yours.

[xulong]

Hi Simon, thanks for the hint. I drew the dispersion results from both cufflinks and DESeq, they are totally different.
Rplot_disp.png is the one from cufflinks, and Rplot_plotDispEsts is from DESeq.
I, however, have no clue how to interpret this difference.

Another possibility is the htseq-count command that messed things up. Below is my code, but I can not see problems from this.
1 #!/bin/sh
2 # RNA-seq data analysis pipeline with Tophat and DESeq
3 # Author: XuLong Wang ([email protected])
4
5 # Step 1: Reads count over gene features
6 # Usage: sh myhtseq-count.sh > logfile
7
8 echo $0;
9 echo "Start: `date`";
10
11 for file_ID in 1 2 3 4 5 6
12 do
13 # Convert .bam files to .sam format
14 # samtools view -h -o ../Tophat-Cufflinks/tophat_out_s$file_ID/accepted_hits.sam \
15 # ../Tophat-Cufflinks/tophat_out_s$file_ID/accepted_hits.bam
16 # echo "BAM to SAM conversion for hs43l8_$file_ID finished!"
17 # Run htseq-count to get count table
18 # mode: union, intersection-strict, intersection-nonempty
19 htseq-count --mode=intersection-nonempty \
20 --stranded=no \
21 ../Tophat-Cufflinks/tophat_out_s$file_ID/accepted_hits.sam \
22 ../UCSC_mm10_genes.gtf \
23 > genes_s$file_ID.txt
24 echo "Reads counts for hs43l8_$file_ID finished!"
25 done
26
27 echo "Stop: `date`";

[Simon Anders]

I did not know that cuffdiff outputs the dispersion estimates, but to me, the value range rather looks like variances than dispersions. If so, calculate the dispersion a from the variance v with a=(v-µ)/µ², where µ is normalized mean count.

Your DESeq dispersion plot looks quite fine to me: For genes with reasonable expression strength, it is below 0.01, i.e., these genes differ typicaly by less than 10% (sqrt(0.01)=0.1) between replicates. So, if they differ by something like 30% or 40% between conditions, you should see significance.

Hence, maybe also post an MA plot (with "plotMA") and a sample clustering heatmap (see vignette).

[xulong]

Hey Simon, cuffdiff ouputs a file called genes.count_tracking, which contains columns named: count, count_variation, and count_dispersion_var. They make the dispersion plot based on count (x) and count_dispersion_var (y). Formula for the dispersion_var is obviously not what you provided. I attached the head of this file here.
Anyway, I calculated the dispersion from the variance column that cuffdiff provided with the formula you showed, and got the new dispersion plot: Rplot_deseq_disp_C.png for control and Rplot_deseq_disp_N.png for treatment.
They are different from the dispersion plot in deseq. Not surprised because cufflinks use different strategy for reads assembling. However, I see almost half of read counts falls below 10. Do not whether this is an issue. This is probably because cuffdiff does not filtered out the low read counts while in deseq I did.
I also attached the dendrogram plot from both cufflinks and deseq, where s1, s2, and s5 from deseq represents Nic_0, Nic_1, and Nic_2 in cufflinks. Do you think it is possible that deseq did not assigned the conditions properly? I also attached my R code for deseq pipeline here.
The MA plot from deseq is also attached here.
Some files will be attached in a separate reply due to the 5 file limit.

[xulong]

Here is my total code for deseq pipeline and head of the gene.count_tracking file.

[Simon Anders]

According to your heat map, samples s1 and s4 are very similar, and samples s2, s3, and s5 form another group of similar samples. So if you compared s1/s4 against s2/s3/s5, you should get some good hits.

However, you seem to be comparing s1/s2/s5 (treated) against s3/s4 (untreated):

cond <- factor(c("treated","treated","untreated","untreated","treated"));

Maybe you accidentally swapped samples s1 and s3 when supplying the data to DESeq?