0
$\begingroup$

I am using RNA seq data to analyze genes via a volcano plot comparing differential gene expression of bacteria with and without antibiotic in R. After having created my plot, I am unsure why some of my values which read "0" TPM in one condition and a TPM value that's not "0" in the other condition were not determined to be differentially expressed. Some of the genes on my volcano plot have this difference in TPM and show up as significant on my plot while others with this "0" value difference are not considered to be significant according to my plot.

Sample of my data: (UO1_D712_##### represents the locus number, the top column represents the different replicates where "1" is untreated and "2" is treated with antibiotic, the letters next to the numbers represent a different technical replicate):

structure(list(`1` = c("U01_D712_00001", 

"U01_D712_00002", "U01_D712_00003", 
"U01_D712_00004", "U01_D712_00005", "U01_D712_00006", "U01_D712_00007", 
"U01_D712_00008", "U01_D712_00009", "U01_D712_00010", "U01_D712_00011", 
"U01_D712_00012", "U01_D712_00013", "U01_D712_00014", "U01_D712_00015", 
"U01_D712_00016", "U01_D712_00017", "U01_D712_00018", "U01_D712_00019", 
"U01_D712_00020", "U01_D712_00021", "U01_D712_00022", "U01_D712_00023", 
"U01_D712_00024", "U01_D712_00025", "U01_D712_00026", "U01_D712_00027", 
"U01_D712_00028", "U01_D712_00029", "U01_D712_00030"), `1a` = c("6.590456502", 
"6.741694688", "7.110342585", "6.299482433", "2.77173982", "2.470330508", 
"3.827125186", "6.267229842", "5.524708663", "2.657913228", "2.87209859", 
"7.479820548", "4.980185572", "4.210955199", "0", "4.492492822", 
"3.611091371", "3.714270433", "7.455036914", "7.045203025", "6.061860857", 
"2.925268313", "6.544077039", "5.747318013", "9.97083444", "9.22523089", 
"4.20205383", "6.097040679", "2.621192351", "0"), `1b` = c("6.544454427", 
"6.601489488", "7.134224619", "5.814043553", "3.280958379", "2.649180803", 
"3.860542083", "6.256648363", "5.380427766", "2.581027705", "3.016132165", 
"7.405329447", "4.701503289", "4.073814818", "0", "4.304196924", 
"3.515977329", "3.843535649", "7.342972625", "6.966606769", "6.122878624", 
"3.007522306", "6.495797641", "5.965431621", "9.828050269", "9.219563915", 
"4.065778989", "6.105331066", "3.061209408", "0"), `1c` = c("6.608196006", 
"6.743010138", "7.102600793", "6.146518601", "3.555184202", "2.364971542", 
"4.034053983", "6.158523627", "5.656051812", "2.658660735", "3.054717455", 
"7.392164473", "4.950953264", "4.277770477", "0", "4.507936666", 
"3.794842979", "3.610794578", "7.471646548", "7.104792624", "6.484767016", 
"3.071205184", "6.584425715", "6.20466015", "9.986342122", "9.282943758", 
"4.179958213", "6.219551653", "2.984738345", "0"), `1d` = c("6.547155382", 
"6.558328892", "6.992501615", "6.449558793", "3.059801464", "2.418800257", 
"3.96498952", "6.013208538", "5.279919645", "2.893295085", "1.750510471", 
"7.408735671", "4.9425624", "3.804549986", "0", "4.174835979", 
"3.806006888", "3.570390524", "7.641137006", "6.976672494", "6.363030106", 
"3.083061726", "6.300910093", "6.007490342", "9.926316442", "9.09671588", 
"4.320556917", "6.153860107", "2.877230446", "0"), `1e` = c("6.626724417", 
"6.577345176", "7.156821278", "6.296873411", "3.618089702", "2.394444986", 
"4.129376392", "6.011715246", "5.31197869", "2.00754706", "2.695493528", 
"7.538910448", "4.606060035", "3.909472643", "0", "4.346616047", 
"3.468681284", "3.338231445", "7.559599613", "7.1527452", "6.232923513", 
"3.108624209", "6.535435309", "6.12922864", "10.13108497", "9.310331313", 
"3.959568571", "6.182335537", "2.902736258", "0"), `1f` = c("6.419219179", 
"6.650459302", "7.319125725", "5.570869357", "2.962845933", "2.55903176", 
"4.087597573", "5.995610538", "5.386268651", "2.750800859", "2.416572678", 
"7.579148955", "3.952633067", "3.615227674", "0", "4.562838935", 
"3.76942104", "3.920096905", "7.935320749", "7.501470652", "6.13700099", 
"3.123910608", "6.035971952", "5.706235015", "9.254751395", "8.379630979", 
"4.51391973", "5.6890651", "2.43285316", "0"), `1g` = c("6.553221221", 
"6.633949852", "7.182305386", "5.769365973", "2.721354972", "1.668390466", 
"4.148367057", "6.240883382", "5.458877133", "1.733842637", "2.723803355", 
"7.522249899", "4.149567197", "3.780763096", "0", "4.496306813", 
"3.645643535", "3.851001768", "7.678552875", "7.283411279", "6.591585956", 
"2.879345378", "6.389427003", "5.911222165", "9.851084493", "9.084575304", 
"4.272587776", "5.974762147", "2.98852705", "0"), `2a` = c("9.769887737", 
"4.550226652", "1.869021464", "4.944848987", "7.9678549", "2.682865013", 
"8.495575559", "2.234521659", "3.667196316", "9.180445037", "3.210107621", 
"7.21523691", "5.714579923", "5.423986751", "9.118981459", "8.635701597", 
"2.742889473", "3.712618983", "8.006057144", "4.999541279", "10.54351774", 
"5.880978085", "7.145433526", "6.982416661", "9.339651188", "5.360835327", 
"4.699680905", "3.423826225", "6.408271885", "5.038170992"), 
    `2b` = c("10.26397519", "1.945664005", "2.086763158", "2.800904763", 
    "8.583418657", "1.536094563", "9.32057547", "3.10685839", 
    "4.224502319", "10.1252842", "1.811175407", "5.714439316", 
    "6.039142559", "3.833361174", "9.757360286", "9.565906731", 
    "1.523640473", "2.315033488", "6.312524363", "4.889986456", 
    "10.23020108", "4.848727685", "7.533256999", "7.138160378", 
    "10.30380331", "3.955469283", "3.167940742", "3.599655687", 
    "4.828945262", "3.701043054"), `2c` = c("9.478481216", "2.789289131", 
    "3.393949527", "3.754810933", "8.710777154", "1.806170784", 
    "9.150005253", "2.612275457", "4.961073313", "9.802701699", 
    "2.933183115", "6.532384958", "6.919449225", "4.432699799", 
    "9.715063475", "9.265691356", "1.412064593", "3.330131873", 
    "6.665979896", "5.158526421", "9.417365584", "4.899531204", 
    "8.173459354", "7.271400938", "9.813068613", "4.384622077", 
    "3.700645365", "4.457874829", "5.649440022", "3.531010379"
    ), `2d` = c("8.199795497", "2.565711524", "2.202287889", 
    "3.856354444", "7.380224849", "2.192476466", "8.14446837", 
    "1.144258862", "3.31447122", "8.713146629", "2.697890381", 
    "6.304428859", "5.745291803", "4.898396114", "9.173362747", 
    "8.339933849", "3.159678152", "4.094587234", "7.608649692", 
    "5.280206424", "10.34630403", "5.098585806", "7.262400625", 
    "6.150190905", "9.316698845", "5.073027993", "4.695003229", 
    "2.485847024", "5.545300465", "4.350571411"), `2e` = c("8.74935033", 
    "3.739489484", "0.217205045", "4.413657999", "8.745525588", 
    "3.657060086", "9.279834921", "2.898951179", "4.282874018", 
    "9.610485827", "3.561102455", "7.228334332", "6.388491443", 
    "4.389908652", "8.781086564", "9.178866581", "2.374603596", 
    "3.961037408", "7.864369809", "3.654728044", "10.15284858", 
    "5.894439123", "7.68020282", "7.12243523", "9.998637438", 
    "5.092174395", "4.111530392", "3.776835632", "5.624523213", 
    "4.095011377"), `2f` = c("7.648310926", "4.345557215", "1.986576876", 
    "4.99426288", "7.087937177", "2.810917253", "7.77706637", 
    "2.62822773", "3.581811188", "8.470225989", "3.335757437", 
    "7.416094847", "5.208841128", "5.536128034", "8.255571138", 
    "7.993319997", "1.9209089", "3.573861828", "7.318814519", 
    "5.233804806", "11.05855833", "6.247720809", "6.673407583", 
    "6.029960625", "8.806867591", "5.459208493", "4.001428729", 
    "3.609936979", "5.876522973", "4.652839671"), `2g` = c("7.555235468", 
    "3.892899549", "1.726443458", "5.304546796", "7.039588042", 
    "3.027235295", "7.703852207", "1.753183519", "2.909288568", 
    "8.385169315", "3.902707541", "7.523315081", "4.978364017", 
    "5.49103181", "8.096218606", "7.944822989", "2.352609608", 
    "4.155433517", "7.227355741", "5.532668321", "11.24946953", 
    "6.159185473", "6.443203375", "5.931761874", "8.7421732", 
    "5.502000205", "4.652883503", "3.458323017", "6.566487449", 
    "4.89790353")), row.names = c(NA, 30L), class = "data.frame")

Design matrix:

structure(list(SampleID = c("1a ", "1b", "1c", 

"1d", "1e", "1f", 
"1g", "2a", "2b", "2c", "2d", "2e", "2f", "2g"), Group = c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), Replicate = c("a", 
"b", "c", "d", "e", "f", "g", "a", "b", "c", "d", "e", "f", "g"
)), row.names = c(NA, 14L), class = "data.frame")

Here is how I wrote my code to sample the RNA seq data (note ~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv refers to the the sample of data above and ~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv refers to the above design matrix):

R10LB_0vs0.5 <- read.csv2("~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv", sep=",", check.names = F)
R10LB_0vs0.5 <- janitor::remove_empty(R10LB_0vs0.5, which = "cols")
R10LB_0vs0.5_design <- read.csv2("~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv", sep=",")

rownames(R10LB_0vs0.5) <- R10LB_0vs0.5$`1`
R10LB_0vs0.5 <- R10LB_0vs0.5[,-1]

#colnames(R10LB_0ugvs0.5ug) <- R10LB_0ugvs0.5ug[1,]
#R10LB_0ugvs0.5ug <- R10LB_0ugvs0.5ug[-1,]

R10LB_0vs0.5 <- data.matrix(R10LB_0vs0.5)

R10LB_0vs0.5 <- DGEList(counts = R10LB_0vs0.5, group = R10LB_0vs0.5_design$Group)
str(R10LB_0vs0.5)

Group <- as.vector(as.character(R10LB_0vs0.5_design$Group))
Replicate <- as.vector((R10LB_0vs0.5_design$Replicate))
R10LB_0vs0.5_designmatrix <- model.matrix(~0+Group+Replicate)

keep <- filterByExpr(R10LB_0vs0.5)
R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]

R10LB_0vs0.5_Disp <- estimateDisp(R10LB_0vs0.5_filter, R10LB_0vs0.5_designmatrix)

CONTRASTS <- makeContrasts( Group1vs2 = Group1-Group2,
                            levels = R10LB_0vs0.5_designmatrix)

R10LB_0vs0.5fit <- glmFit(R10LB_0vs0.5_Disp, contrast = CONTRASTS[,1])
R10LB_0vs0.5lrt <- glmLRT(R10LB_0vs0.5fit)
R10LB_0vs0.5TT <- topTags(R10LB_0vs0.5lrt, n=nrow(R10LB_0vs0.5_Disp))
write.csv(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.csv")
saveRDS(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.RDS")

Volcano plot:

ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = PValue)) + geom_point()

p4 <- ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = -log10(PValue))) + geom_point() + theme_minimal()

#Add vertical lines for logFC thresholds, and one horizontal line for the p-value threshold
p5 <- p4 + geom_vline(xintercept = c(-0.6, 0.6), col = "red") + geom_hline(yintercept = -log10(0.05), col = "red")


# add a column of NAs
R10LB_0vs0.5_comparison$diffexpressed <- "NO"
# if logFC > 0.6 and PValue < 0.05, set as "UP" 
R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC > 0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "UP"
# if logFC < -0.6 and PValue < 0.05, set as "DOWN"
R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC < -0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "DOWN"

p4 <- ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed)) + geom_point() + theme_minimal()

p5 <- p4 + geom_vline(xintercept=c(-0.6, 0.6), col="red") + geom_hline(yintercept=-log10(0.05), col="red")

p6 <- p5 + scale_color_manual(values=c("blue", "black", "red"))
p6 <- p5 + scale_colour_manual(values = mycolors)

R10LB_0vs0.5_comparison$delabel <- NA
R10LB_0vs0.5_comparison$delabel[R10LB_0vs0.5_comparison$diffexpressed != "NO"] <- R10LB_0vs0.5_comparison$gene_symbol[R10LB_0vs0.5_comparison$diffexpressed != "NO"]

ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) + geom_point() + theme_minimal() +geom_text()

library(ggrepel)
ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) +geom_point() + theme_minimal() +geom_text_repel() + scale_color_manual(values=c("blue", "black", "red")) + geom_vline(xintercept=c(-0.6, 0.6), col="red") +geom_hline(yintercept=-log10(0.05), col="red")

This is the volcano plot I currently have with all of my data: volcano plot

Can you explain how some of these functions filtered my data (i.e. filterByExpr, estimateDisp, makeContrasts, glmFit, glmLRT)? Why are some of my data points that change from 0 TPM in one condition to some value in another condition showing up on my plot while others are not? Or should I accept that my data analysis is valid or "correct"?

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ I heard that low expression genes are often filtered out in a first past due to a high chance of false positives, so they could be dropped from filterByExpr. High variance within a condition might undermine the statistical significance of the mean difference across conditions. Is there a way to look at the variance of my data, and is there a threshold for variance between my replicates that causes them to be filtered out by filterByExpr or other functions? $\endgroup$
    – ry0-H
    Mar 12, 2021 at 20:57
  • $\begingroup$ It is a long time that I don't use R. However, I had a sinilar problem once. Using not cartesian coordinates, ggplot cut several points out of the plot. If you plan to continue using R, search about this problem. Otherwise, if you can accede to it go with MATLAB that has the best tool for volcano plot, by my opinion. It never fails. $\endgroup$
    – Spartan 117
    Jan 27, 2022 at 22:23

2 Answers 2

Reset to default
1
$\begingroup$

edgeR functions do not do any filtering, i.e., they do not remove any genes from the results or from the plots. The only filtering was that which you did explicitly yourself by:

keep <- filterByExpr(R10LB_0vs0.5)
R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]

For your data, filterByExpr will keep any gene that is expressed at a minimal level in at least 7 samples. You could relax that requirement, to allow genes to be expressed in only 5 or 6 samples say, by reducing the large.n argument to filterByExpr.

Improve this answer
$\endgroup$
0
$\begingroup$

[disclaimer: I understand DESeq2, not EdgeR]

Most likely, low expression values (as in raw counts, prior to any normalisation such as TPM). Expression that has high dispersion across all samples is less likely to be biologically significant.

It's easier to work out what's going on with an MA plot, with average gene expression on the X axis, and log2 fold change on the Y axis.

I'm used to DESeq2 results, so I'd guess the construction would be something like this (assuming that the 'baseMean' field exists in EdgeR results):

R10LB_0vs0.5lrt %>%
    ggplot() +
    aes(y = logFC, col = padj, x=log(baseMean)) +
    geom_point()

It's possible that EdgeR doesn't have anything like baseMean in its default results. As a proxy for this, you could use the count totals from the input matrix (e.g. R10LB_0vs0.5 %>% rowSums() %>% log2()).

Improve this answer
$\endgroup$
1
  • $\begingroup$ Do you mean that the TPM values are too low? Because I believe the ones that have been filtered out have a similar value to other "0" tpm genes that were significant. Also if I use the average of two conditions, will the genes I pull out be less significant because the MA plot will not account for variance between replicates? Thanks for the answer $\endgroup$
    – ry0-H
    Mar 12, 2021 at 7:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.