Only one gene is upregulated?

Question

[mskr_]

I just start learning bioinformatics from RNA-seq/DEGs. I am trying to replicate an article. I very much appreciate any help. I used salmon in combination with tximport and DESeq2. PCA plot showed 77%PC1 8%PC2. Apparently, different conditions are mixed and not well separated. Is it healthy to continue this analysis?

library(tidyverse) 
library(tximport) 
library(DESeq2)

sample_table <- read.csv("SraTable.txt") %>% arrange(sample_table$Run) %>% select(Run)

sample_table$Run <- c("pre_ischemia1","ischemia1","pre_ischemia2","ischemia2","pre_ischemia3", "ischemia3","pre_ischemia4","ischemia4","pre_ischemia5","ischemia5")

SRR10700832_quant/quant.sf sample_files <- paste0(c("SRR10700832_quant","SRR10700833_quant","SRR10700835_quant","SRR10700836_quant","SRR10700838_quant", "SRR10700839_quant","SRR10700841_quant","SRR10700842_quant","SRR10700844_quant","SRR10700846_quant"), "/quant.sf")

names(sample_files) <- pull(sample_table, Run)
gene_map <- read.csv("gene_map.csv", col.names = c("enstid", "ensgid")) 
count_data <- tximport(files = sample_files, type="salmon", tx2gene= gene_map, ignoreTxVersion = TRUE)
sample_table <- as.data.frame(sample_table)
sample_table$condition <- factor(rep(c("pre_ischemia", "ischemia"), times=5), levels = c("pre_ischemia", "ischemia"))

deseq_dataset <- DESeqDataSetFromTximport(txi= count_data, colData = sample_table, design = ~condition) 
deseq_dataset <- DESeq(deseq_dataset)
vst <- varianceStabilizingTransformation(deseq_dataset)

plotPCA(vst, intgroup='condition')
results <- results(deseq_dataset)

summary(results)

res <- results(deseq_dataset) yielded following;

out of 38464 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 1, 0.0026%
LFC < 0 (down)     : 0, 0%
outliers [1]       : 263, 0.68%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

There is only 1 gene upregulated, zero down and zero low counts which are not likely. Which step could be the problem ? I would like to know sanity checkpoint in the entire flow.

Thank you for your time

P.S sorry it this hurts your eyes. I couldnt format the code in the post.

[cpad0112]

what is the result of plotPCA(vst, intgroup='condition') from your code?

[mskr_]

Hi, 77%PC1 8%PC2 I added the plot

[tothepoint]

p-value < 0.1

[mskr_]

Hi thank you for the reply, I changed it to 0.05, but there is no difference since there is one gene that fulfills that condition.

[swbarnes2]

How can anyone help you with code we cannot see?

[mskr_]

Hi thank you for your reply, since it didn't give any error I thought there is something fundamentally wrong. I updated the code

[ATpoint]

@mskr_ You can enable proper code formatting following the markdown rules, so wrapping code blocks in ``` or by selecting the code chunk and then using the 10101 button. We made the changes for you now :)

[mskr_]

Hi, I placed the long code in a place after clicking 10101 button but it didnt update. Just trying to figure out, sorry and thank you :)

Answer 1

[ATpoint]

Hello,

I just played with these data a bit as the salmon quantification directories were available at GEO, so this was no effort.

It seems the authors filtered on p-value based on the methods text...

The significantly differentially expressed genes (DEGs) between groups were defined at cut-off criteria of |log2 fold-change| ≥ 1 and p-value < 0.05.

...instead of padj. When you check some of the top20 DEGs that they provide in Supp. Table S3 you will see that these "significant" ones only pass the p-value cutoff but all have a padj >>> 0.05, so there indeed simply are no DEGs at standard cutoffs, at least not without correcting for potential batch effects. I did not look into the paper at any detail, but this at least seems to "replicate" the results they have. Not discussing now whether this strategy they chose was good or not.

[mskr_]

ATpoint thank you for pointing out this detail. I totally missed that! Very much appreciated!

Answer 2

[Carlo Yague]

Hi mskr_,

Your interpretation of the PCA plot is correct. The dispersion of samples in the PC1-PC2 space, not clustering according to condition, indicates that the residual variability is much higher than the factorial variability (condition effect) in your sample. In consequence, it is not surprising that you find only one DEG in the differential expression analysis.

There are multiple hypothesis that might explain the lack of condition effect:

• There could be truly few difference between the transcriptome of pre-ischemia and ischemia cells. Not much you can do about this.
• The sample labeling might have been mixed up. You need to double check that SRR10700832 is really pre-ischemia and so forth.
• There could be another source of variability (batch effect, etc...) that strongly affects the measure of the transcriptome and hides the true signal. If you identify such an explainable source of variability, you might want to include that in the model to control for it which will possibly improve your results.

In the article you try to replicate, were there many DEG ? Have you read in detail how they analyzed their data ?

[mskr_]

Thank you for the detail explanation. As ATpoint pointed out , they filter based on pvalue, rather than padj. Then it fits.

Just a curiosity what pc1 /pc2 values are off the limits? Or you would consider too much?

[Carlo Yague]

Units of PC are difficult to interprets since they reflect the combination of multiple variable (thousands variable in this case). So usually I don't really look at PC values, just at the percentage of variance explained by each component, and how the samples are distributed according to the PCs. If you want more insight into PCA interpretation in the context of RNA-seq, there are multiple resources out there, including this one : https://hbctraining.github.io/DGE_workshop/lessons/03_DGE_QC_analysis.html