Differences between published differential gene expression results and own analysis on RNA-seq data

Question

[vernonlim98]

Does anyone else experience their results of differential gene expression analysis being vastly different from what has been published? I am still very new to R programming and bioinformatics, and am just trying to find differentially expressed genes between platinum resistant and sensitive samples in the TCGA ovarian cancer dataset.

When I try to run either limma or DESeq2, I cannot seem to replicate the results that have been published by multiple papers, even when I use the same datasets and try to follow what has been done in their methodology, and their code doesn't seem to be publicly published. Objectively, I would trust published results over my own analysis, but when I try to do individual t-tests on the genes of interest, the results tend to lie closer to my own analysis compared to what has been published.

Anyone else facing the same issue, or have any possible insights? Any help will be greatly appreciated

[ATpoint]

Please link the paper, lets see their methods. With human data it can well be that they included something like surrogate variables to dampen unwanted variation, or tested against fold changes, and if so then results can be very different compared to a standard pairwise analysis. That all assumes that your code and input was correct, so one would need to see your minimal code too. Otherwise it's just guessing.

Answer 1

[yura.grabovska]

When you say vastly different, what do you mean exactly? Does your top 100 up/down not match theirs at all?

[vernonlim98]

Yes it does not match at all, and those genes that are common in both analysis results have very different logFC and p values

Answer 2

[Ales]

I second the comment about comparing the top 100 genes - if you see vastly different results, I would lean closer to examining protocols themselves and verifying there are no hidden errors.

If the top 100 show some similarity, the differences might be due to a number of reasons:

1. Software. take for example pseudoalignment methods (Salmon/Kallisto) vs spliced alignment (Hisat2/STAR). The former will bias the results more towards what is in your known/reference transcriptome and will sometimes force allocate extra reads into known genes (eg. multimapping reads from homologous gene not included in your reference)
2. Subset of the previous point - parameters matter a lot, and alas, are sometimes poorly documented. As you commented - if there are no scripts available to supplement the analysis, your guess for which versions and arguments were used might not be entirely accurate. For example, in Salmon - did they augment the index with decoy sequences or not? This single difference could amount to considerable difference in the results alone.
3. Preprocessing of the data - did you perform the same filtering of the data before performing the analysis as the authors?
4. Reference gene annotations can differ very (and I mean very) significantly not just from source to source (RefSeq vs Gencode) but from version to version (especially true for Ensembl and Gencode). If you are working with anything other than human/mouse/drosophila/thaliana - the annotations are probably even more divergent and using a version/source different from the original study will likely yield vastly different results. Tangentially, is the same version of the genome assembly being used? Many human studies (especially clinical) still rely on the outdated GRCh37 reference genome assembly.

There are many other reasons others might offer below, but these are the first few I would look into (especially genome annotations, since it often goes overlooked under the assumption that its a solved problem). There's a reason so many papers get published comparing different methods for the same computational task - often times there are considerable differences at all parts of the analysis and many go overlooked.

[dsull]

Honestly, those things shouldn't result in a considerable discrepancy.

The largest difference will come from the post-processing & analysis (normalization type, limma vs. deseq2, batch correction methodology, etc.).

Just from anecdotal experience.

For the OP: Start by making plots (plot your logFC and p-values against what the paper reports, make Venn Diagram overlaps of DE genes, etc.).

[vernonlim98]

Thanks both for your input. It got me to look at the published methodologies more closely, and you guys were right. One paper that I was referencing used FPKM data from TCGA and standardised with quantiles algorithm use preprocessCore before doing limma, and another paper used TPM data directly for limma. I'm no expert in R and bioinformatics by any means (have been self-learning just a few months) but I know that limma should be using raw counts and limma-voom transforms into CPM as normalisation.

Maybe that's why my results are so different, because I'm just using raw counts for limma-voom?

[dsull]

If those logFCs are different, then that's a normalization issue.

The logFC is straight forward arithmetic (it's simply dividing mean of [normalized] counts in one condition by mean of [normalized] counts in another condition, and taking the log).

For those genes that are discrepant, try calculating the logFCs by hand, attempting different methods of normalization (raw counts, CPM, FPKM, TPM, deseq, etc.)

[vernonlim98]

So is there any correct way for normalisation? I know limma recommends CPM via limma-voom and DESeq2 just takes in raw counts and does their own normalisation. I'm just confused why published papers use already normalised counts for limma/deseq when I think it's fairly established raw counts should be used.

[dsull]

There are many ways to do normalization and one of the correct ways is simply to plug raw counts into deseq2 (the math that deseq2 uses to perform normalization is very straightforward).

I’m not aware of any paper that plugs normalized counts into deseq2.