Hi bioconductor community!

I have performed a DE analysis using DESeq2, with conditions being carriers vs. non-carriers of a risk allele identifed from a GWAS study. The sample size is 131 with 81 and 49 in each group. I get 18 down-regulated and 726 upregulated. The effect sizes have been shrunk by lfc. Of the up-regulated genes, the maximum logFC is 22.9 and minimum is 0.44 with mean 2.26 logFC. Of the down-regulated genes we have a range from -9.77 to -0.57 with mean -2.0 logFC.

In my opinion there is no overall trend in the MA plot, but it is extremely assymetric in terms of so many upregulated genes vs. downregulated, which I think is worrying. I have performed this experiment on 20-some genetic variants using the same expression dataset, with some of the variants showing the same behaviour. However, a lot of variants have completely symmetric MA Plots. The median ratio between up/ and down-regulated genes is 0.6 meaning that there are ca. twice as many downregulated genes on average. There are some extreme outliers with 30-40 times as many upregulated, and some with 50-70 times as many downregulated genes. But most (75 %) lie within 0.14 to 3 times as many upregulated genes vs. downregulated. The number of samples in each group is not widely different among the variants. It is quite difficult to say if this is biologically meaningful, given the limited knowledge about the genetic variants in general.

I just wanted to get your opinion on the matter. Does this MA plot indicate that something is wrong with the normalization or could it be due to low sample size, or something else?

Here is my code to generate the MA plot. The levels of the genetic variants are 1: carrier or 0: non-carrier. The column is called "SNP".

    dds <- DESeqDataSetFromTximport(txi = mRNA_list,
                                    colData = coldata,
                                    design = ~ SNP + age + gender)
    dds <- DESeq(dds)
    res <- results(dds)
    res <- results(dds, contrast=c("SNP","1","0"))

    resLFC <- lfcShrink(dds, contrast=c("SNP","1","0"), type="normal")
    plotMA(resLFC, ylim=c(-2,2))

so many upregulated genes vs. downregulated, which I think is worrying

Take a look at the plotCounts for a couple of the up-regulated genes, and a couple of the down-regulated genes.

You can identify these with:

with(res, plot(baseMean, log2FoldChange, log="x"))
with(res, identify(baseMean, log2FoldChange, log="x"))

My guess those is that there is something going on with the carriers, beyond just the risk allele. Were they treated in different batches and is that modeled here?