SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
What is the meaning density in CummeRbund plots? Carlos Borroto Bioinformatics 16 12-29-2016 09:05 PM
volcano plots cummeRbund godzilla07 Bioinformatics 8 02-04-2015 12:55 PM
Strange Dispersion and Density Plots in cummeRbund DrS Bioinformatics 5 03-25-2013 05:24 AM
CummeRbund plots godzilla07 Bioinformatics 1 08-10-2012 04:00 AM

Reply
 
Thread Tools
Old 06-14-2013, 02:04 PM   #1
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default Fruit Fly CummeRbund Plots Are Not The Same

Hello all,

After several days of moving through the TopHat Cufflinks workflow with the example fruit fly data (http://www.ncbi.nlm.nih.gov/pubmed/22383036), I have finally made it to CummeRbund.

When I execute the following commands in R, I get plots that are not the same as the ones found in the publication. Is that a problem?
library(cummeRbund)
cuff_data <- readCufflinks('diff_out')
csDensity(genes(cuff_data))
csScatter(genes(cuff_data), 'C1', 'C2')
csVolcano(genes(cuff_data), 'C1', 'C2')

Attached is the volcano plot. I just took a screen shot and saved it as ppt since I didn't know how to export or save the R graph.

Thanks and God bless,
Jason
Attached Files
File Type: ppt 2013-6-14_DiffGeneValcano.ppt (239.0 KB, 51 views)
jmwhitha is offline   Reply With Quote
Old 06-14-2013, 02:08 PM   #2
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default csDensity is really different

I also wanted to share my csDensity plot because it is really different.
Attached Files
File Type: ppt 2013-6-14_DiffGenecsDensity.ppt (287.5 KB, 36 views)
jmwhitha is offline   Reply With Quote
Old 06-14-2013, 02:20 PM   #3
jparsons
Member
 
Location: SF Bay Area

Join Date: Feb 2012
Posts: 62
Default

This does look like a problem. Assuming that you are using the right original data (GSE32038), you should get the same volcano plot. Most concerning to me is that you're missing any significant genes.

I imagine that you are using updated versions of the software relative to the paper, but in principle this shouldn't change the results so dramatically. Several default options have changed recently, though. Without looking over everything, it's tough to say for sure.

An interesting and potentially educational for you approach would be to run the same data through Galaxy's tophat/cufflinks/cummerbund pipeline to see if you get the same results (or if perhaps you or the tutorial had done something incorrectly).
jparsons is offline   Reply With Quote
Old 06-17-2013, 10:03 AM   #4
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default

I'd like to check the assumption of using the correct data with you.

I went to http://0-www.ncbi.nlm.nih.gov.elis.t...i?acc=GSE32038, and copied the links for the RAW.tar and simulated_fastq_files.tar.gz.

My commands were:
(in GSE32038 directory)
wget ftp://ftp.ncbi.nlm.nih.gov/geo/serie...E32038_RAW.tar -O GSE32038_RAW.tar
tar -xvf GSE32038_RAW.tar
gunzip GSM79448*
(in the my_rnaseq_exp directory)
wget ftp://ftp.ncbi.nlm.nih.gov/geo/serie...q_files.tar.gz -O GSE32038_simulated_fastq_files.tar.gz #Could not find any other fastq files in the GEO Assession
tar -xvf GSE32038_simulated_fastq_files.tar.gz
gunzip GSM79448*

The TopHat Cufflinks journal article says to download the "raw sequencing reads, alignment reads, assembled transfrags and differential analysis", but I was not exactly sure which of files were the ones I needed. Did I download the correct ones?

Thank you!
Jason
jmwhitha is offline   Reply With Quote
Old 06-17-2013, 11:36 AM   #5
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,077
Default

Based on the name of the file (GSE32038_simulated_fastq_files.tar.gz) you seem to have downloaded the right data file.

Figure 8 in the paper has different numeric ranges for axes (specially the Y-axis) that may be why your volcano plot is looking different. We can't eliminate the possibility that you have different results in the underlying data. Same is true for the other graph.

Also keep in mind that the example in the nature protocols paper used previous versions of the Tophat/Cufflinks suite. You are probably using the newer versions (v.2.x).

Last edited by GenoMax; 06-17-2013 at 11:38 AM.
GenoMax is offline   Reply With Quote
Old 06-17-2013, 03:58 PM   #6
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default

Okay, and for the fruit fly genome I went to ftp://ussd-ftp.illumina.com/.
wget ftp://igenome:G3nom3s4u@ussd-ftp.ill...DGP5.25.tar.gz -O Drosophila_melanogaster_Ensembl_BDGP5.25.tar.gz
tar zxvf Drosophila_melanogaster_Ensembl_BDGP5.25.tar.gz

Is that also correct?

I also noticed the axis were different in Figure 6. If FPKM (fragments per kilobase of exon per million fragments mapped) is a decimal, log(FPKM) is negative. This is what we see in my figure but not theirs. If FPKM is a decimal than that means there are fewer fragments per kb of exon than million fragments mapped.

In a previous thread (http://seqanswers.com/forums/showthread.php?t=31074), I posted several warnings that I got from ran Cuffmerge. These warnings were for missing fasta records such as heterochromatin and mitochondrial genome. Perhaps, without these index files, fragments mapped incorrectly all the way back at the TopHat step. This may account for negative values for log10(FPKM) in the csDensity graph.

If that is the case, for me to get the same results as in the publication, I will need a fasta index with everything including heterochromatin and mitochondrial genome. If the website (ftp://ussd-ftp.illumina.com/) is not the correct source, does anyone know what is the website with the full genome?
jmwhitha is offline   Reply With Quote
Old 07-08-2013, 05:08 AM   #7
rubbertjes
Member
 
Location: Netherlands

Join Date: Jun 2013
Posts: 13
Default

I actually have similar problems to yours.
I also analyzed the simulated data, in combination with the iGenomes data. My results also differ from the results shown in the paper. However, my results also differ from yours.
The reason you're not showing any significant genes in your volcano plots, is because you did not switch on the significant genes. For some reason you have to do this explicitly:

e.g. csVolcano(genes(cuff_data), 'WT', 'KO', alpha=0.05, showSignificant=T)

My error bars/standard deviation is also really high for each gene. Much bigger than one standard deviation, I don't know why this is.

Did you ever figure out if you did something wrong, or what the reason for the changes are?

Cheers

Last edited by rubbertjes; 07-08-2013 at 08:42 AM.
rubbertjes is offline   Reply With Quote
Old 07-08-2013, 08:20 AM   #8
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default

Hello rubbertjes,

I was never able to test my hypothesis since I could never find (or was given) a full fruit fly fasta index with the heterochromatin and mitochondrial genomes.

Thanks for tip about showSignificant=T. I will give that a try this week to see what happens.

God bless
jmwhitha is offline   Reply With Quote
Old 09-26-2014, 07:35 AM   #9
Santi
Junior Member
 
Location: Buenos Aires

Join Date: Sep 2014
Posts: 2
Default

Hello jmwhitha,

Did you finally guessed what were the reasons for the different plot outcomes? I am getting exactly the same "wrong" plots as you. If I run cummerbund using the reference result files from the tutorial, some of the plots are corrected (such as the volcano and the scatter plots), but others are not (eg the density plot and the barplots, which only show the error bars). I get some warnings/errors from cummerbund for these last commands (these, even when using the result data provided as reference, thus might be something due to R packages version I guess):

> csDensity(genes(cuff_data))
Warning messages:
1: Removed 4075 rows containing non-finite values (stat_density).
2: Removed 4082 rows containing non-finite values (stat_density).

> mygene <- getGene(cuff_data,'regucalcin')
> expressionBarplot(mygene)
Error : Mapping a variable to y and also using stat="bin".
With stat="bin", it will attempt to set the y value to the count of cases in each group.
This can result in unexpected behavior and will not be allowed in a future version of ggplot2.
If you want y to represent counts of cases, use stat="bin" and don't map a variable to y.
If you want y to represent values in the data, use stat="identity".
See ?geom_bar for examples. (Defunct; last used in version 0.9.2)


I am using Tophat2, cufflinks2 and R3.1.1 for these analysis, while I think in the tutorial they used older versions for these tools...however, I guess resutls should not differ so much, right?

Santi
Santi is offline   Reply With Quote
Old 09-26-2014, 08:50 AM   #10
jmwhitha
Senior Member
 
Location: NC State, Raleigh, NC

Join Date: Mar 2013
Posts: 107
Default

Hello Santi,

I found a variety of issues with this pipeline so I went with DEGseq. Besides not having as many issues, the processing time was much faster. So, actually I'd recommend you try DEGseq or another pipeline rather than trying to get TopHat Cufflinks... to work.

God bless,
Jason
jmwhitha is offline   Reply With Quote
Old 09-26-2014, 09:10 AM   #11
Santi
Junior Member
 
Location: Buenos Aires

Join Date: Sep 2014
Posts: 2
Default

Hi Jason,

Thanks for the feedback! Cheers,

Santi
Santi is offline   Reply With Quote
Reply

Tags
csdensity, cufflinks, cummerbund, diff_out, volcano

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 07:20 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO