I'm planning to use EdgeR for differential expression analysis and I read that using Gibbs sampling improves EdgeR's ability to correction for dispersion. I ran Salmon with the --numGibbsSamples 100 argument, but I can't find the counts for the Gibbs samples. What I do see is a subdirectory 'bootstrap' in 'aux_info' that contains two files: bootstraps.gz and names.tsv.gz. The former appears to be binary - if I use zcat or gunzip -c, gibberish appears on screen, but the former, although it has a .gz extension, appears to be an uncompressed text file. (I can view with less without uncompressing, anyway.) Is the boostraps.gz file supposed to be binary or does the gibberish indicate that my files are corrupted?
Are these the files from Gibbs sampling that EdgeR uses to correct the quantification of the biological samples in the quant.sf files, or are there other files that I have overlooked?
Thanks and best wishes,
Tonya
•
written
by
Tonya S. •
2 answers
edgeR's catchSalmon function automatically reads the aux_info files on the fly. You do not need to parse these files yourself. You just give catchSalmon the name of the directory containing the Salmon output, and it does the rest.
However, if you are working with a poorly annotated organism, and you don't have any confidence that your list of annotated transcripts is comprehensive, and your interest is in genes rather than transcripts, then I would recommend that you simply add up the Salmon counts by gene and do a regular gene-level analysis.
The transcript dispersions from catchSalmon do not help you to do a gene level DE analysis. It is possible to use Gibbs samples from Salmon to compute gene-level ambiguity dispersions, but the gene level dispersions tend to be low so I have not felt it to be worthwhile.
Log in to answer this question.
...on transcript level, are you aware of that? Not on gene level, and 99% of DE analysis are done on gene level, so not resolved by transcripts per gene. Do you really want and need transcript level, as it is notably more difficult to interpret than gene level?
Ah, so improving on the transcript level will not also improve on gene level. I guess that would make sense since all of the reads would just go to the gene - it wouldn't matter at all how they are subdivided among transcripts.
My transcripts were reconstructed using stringtie, so I wouldn't place much stock in any individual transcript.
There is gene- and transcript level. Gene level means that your counts are for genes not transcripts. Methods like STAR + featureCounts give gene level counts right away. Then there are the trascript quantifiers, such as salmon, kallisto, stringtie, etc, that give transcript abundances. This can be summarizes (essentially summed with some extras) to gene level using e.g. tximport. For DE you can ask whether your transcript levels are different, or the gene level is different. The latter is by far the most common. And for the latter taking into account mapping ambiguity between transcripts (because that is where most ambiguity comes from, because transcripts of a gene share most exons) can help DE inference. For gene level, it is not common to use this info, but there are few tools like swish (and probably edgeR with a similar strategy as they do for transcript level) that can use this information. You probably want a normal gene level analysis.