Hello BioC community,

Let me set the stage for my question first. We generated a list of 200 genes that correlate with certain clinical parameters using human RNAseq data. I'd like to analyze the expression patterns of these genes in published GEO datasets. One of the GEO datasets I'm interested in analyzing involves RNAseq of sorted human PBMCs (GSE107011). There are 30 different types of immune cells in this dataset with multiple biological replicates (127 samples total). The data I downloaded using GEOquery package are in TPM format as reported by Kallisto (somehow I can't find the raw data). I can generate heatmaps like this to visualize the gene expression patterns (samples in rows, genes in columns):

I would like to analyze the significantly different genes between data subsets and generate a summary table. I used limma package before for routine DE analyses, but I'm trying to develop more of a targeted analysis approach to investigate the ~200 genes (I don't care much about the other genes in the dataset for now). I performed the analysis, but the results I'm getting do not make sense to me.

These are the steps I performed, please let me know if my thinking is erroneous at any point.

Prepare expression data and sample metadata

# Read data
dat <- read.delim(gzfile("GSE107011/GSE107011_Processed_data_TPM.txt.gz"), check.names = F)

dat[1:5,1:5]
##                     DZQV_CD8_naive DZQV_CD8_CM DZQV_CD8_EM DZQV_CD8_TE
## 1 ENSG00000223972.5       0.214321     0.00000  0.00421414  0.00953342
## 2 ENSG00000227232.5       4.759630     3.49280  2.64877000  3.72232000
## 3 ENSG00000278267.1       0.000000     2.61518  2.45976000  0.00000000
## 4 ENSG00000243485.5       0.000000     0.00000  0.00000000  0.00000000
## 5 ENSG00000284332.1       0.000000     0.00000  0.00000000  0.00000000

# Do some processing including:
### Averaging counts from multiple transcripts matching the same gene
### Log2 transforming the linear TPM values
### Make a ggplot-friendly composite dataframe that contains info from sample metadata and gene expression (will be useful later on for graphing)

expdat[1:5,1:5]
##           sample                     cell_type gender  DNAJC11   CDK11A
## 1 DZQV_CD8_naive             Naive_CD8_T_cells      M 5.609300 5.799777
## 2    DZQV_CD8_CM     Central_memory_CD8_T_cell      M 5.547006 6.204563
## 3    DZQV_CD8_EM   Effector_memory_CD8_T_cells      M 5.778402 6.187289
## 4    DZQV_CD8_TE Terminal_effector_CD8_T_cells      M 5.186984 6.465826
## 5      DZQV_MAIT                    MAIT_cells      M 5.251571 6.152436

Subset the dataset to contain the genes of interest

I'm thinking that this will reduce the number of multiple comparisons (by discarding uninteresting genes), giving the analysis more power to detect differential expression.

# genes object is a character vector containing the names of genes-of-interest
genes[1:5]
## "PPP3CB_AS1"  "IL10RB_AS1"  "PPP1R26_AS1" "RRN3P2"      "AC099343.3" 

# Subset the data frame to contain info from genes and cell_type (group variable from sample metadata)
expsubset <- expdat[, colnames(expdat) %in% c("cell_type", genes)]

dim(expsubset)
## 127 200

# Groups in comparison
groups <- expsubset$cell_type

groups[1:5]
## [1] Naive_CD8_T_cells           Central_memory_CD8_T_cell   Effector_memory_CD8_T_cells
## 30 Levels: Central_memory_CD8_T_cell Classical_monocytes ... Vd2_gd_T_cells

# Design for DE analyses
design <- model.matrix(~0 + groups)
colnames(design) <- gsub("^groups", "", colnames(design))

# Extract expression data from data subset (discard cell_type column) and transform so samples are in columns
de_df <- as.data.frame(t(expsubset[, colnames(expsubset) != "cell_type"]))

rownames(de_df)[1:5]
## [1] "LINC00337"  "PIK3CD_AS2" "MATN1_AS1"  "LINC01778"  "AL353622.1"
colnames(de_df)[1:5]
## [1] "V1" "V2" "V3" "V4" "V5"

# Reassign column names. There will be duplicate column names because of biological replicates
colnames(de_df) <- expsubset$cell_type

colnames(de_df)[1:4]
## [1] "Naive_CD8_T_cells"             "Central_memory_CD8_T_cell"     "Effector_memory_CD8_T_cells"   "Terminal_effector_CD8_T_cells"

Limma analysis

# Create contrasts matrix
cont.matrix <- makeContrasts(Th1_vs_Nai4 = Th1_cells- Naive_CD4_T_cells,
                              Nai4_vs_Th1 = Naive_CD4_T_cells- Th1_cells, levels=colnames(design))  

# Output of cont matrix
cont.matrix
##                                Contrasts
##   Levels                          Th1_vs_Nai4 Nai4_vs_Th1
##   Central_memory_CD8_T_cell               0           0
## ...
##  Naive_CD4_T_cells                      -1           1
##  Naive_CD8_T_cells                       0           0
##  Natural_killer_cells                    0           0
##  Th1.Th17_cells                          0           0
##  Th1_cells                               1          -1
##  Th17_cells                              0           0
##  Th2_cells                               0           0
##  Vd2_gd_T_cells                          0           0             


# Linear fit
lmfit <- lmFit(de_df, design)
lmfit <- eBayes(lmfit, trend=T)

Tabulate results

topTable(lmfit, coef = which(colnames(cont.matrix2)=="Th1_vs_Nai4"), 
               p.value = 0.01, number = 4)
##               logFC  AveExpr        t      P.Value    adj.P.Val        B
## EBLN3P     6.960910 6.450265 48.99890 1.308402e-72 2.420544e-70 153.2290
## PSMB8_AS1  7.144795 7.284830 44.25446 2.626655e-68 2.429656e-66 143.9306
## PCED1B_AS1 8.126577 6.757791 40.95170 4.656414e-65 2.871455e-63 136.8323
## ITGB2_AS1  6.614067 5.535769 39.23635 2.807513e-63 1.298475e-61 132.9173

topTable(lmfit, coef = which(colnames(cont.matrix2)=="Nai4_vs_Th1"), 
               p.value = 0.01, number = 4)
##              logFC  AveExpr        t      P.Value    adj.P.Val        B
## PSMB8_AS1 7.103732 7.284830 44.00011 4.588982e-68 8.489617e-66 143.4033
## EBLN3P    5.927469 6.450265 41.72436 7.720048e-66 7.141045e-64 138.5431
## ITGB2_AS1 5.778358 5.535769 34.27871 1.008373e-57 6.218298e-56 120.6026
## VIM_AS1   6.128814 5.022750 32.85078 5.375385e-56 2.486116e-54 116.7481

The problem

The genes listed to be differentially expressed between two groups are pretty much the same regardless of which groups are compared or which direction I set up the comparisons (i.e. grp1-grp2 vs grp2-grp1). This doesn't make sense. For instance, take EBLN3P, which seems to be always "upregulated". When I plot the expression in a scatter plot, there is clearly no appreciable change between Th1 and Naive CD4 groups:

library(ggpubr)

ggdotchart(expdat, "cell_type", "EBLN3P", 
           color="gender", 
           sort="desc", add="segment", legend="top", rotate=T)

When I calculate the mean expression value, these "upregulated" genes are among the highest in the data frame:

meanvec <- sapply(expdat[, !colnames(expdat) %in% c("cell_type", "gender", "sample") & colnames(expdat) %in% lncs$gene_name], mean)

head(sort(meanvec, decreasing = T), 10)

## PSMB8_AS1 PCED1B_AS1     EBLN3P       HCP5 AL390728.6  LINC_PINT  ITGB2_AS1      
##   7.284830   6.757791   6.450265   6.195957   5.891494   5.762214   5.535769   

What am I doing wrong here?

I've got a few points / comments.

The first one I think will address the problem you are asking about directly. The biggest error here is that you have defined a contrast matrix, but you aren't using it. I'll take bits and pieces from your code, hopefully it makes sense.

cont.matrix <- makeContrasts(Th1_vs_Nai4 = Th1_cells- Naive_CD4_T_cells,
                             Nai4_vs_Th1 = Naive_CD4_T_cells- Th1_cells,
                             levels=colnames(design))
lmfit <- lmFit(de_df, design)
lmfit <- contrasts.fit(lmfit, cont.matrix)   # **you forgot this part**
lmfit <- eBayes(lmfit, trend = TRUE)

res <- topTable(lmfit, "Th1_vs_Nai4", n = Inf)

Testing the other contrast you have should just give you the same results as the first, but with the sign switched for the logFC.

Some other things you can consider:

1. I don't get what the "Averaging counts from multiple transcripts matching the same gene" step is in your pre-processing mojo. Looks like these reads are already rolled up into genes: I don't think you should be doing whatever it is else that you're doing here with this TPM data.
2. I'd suggest taking the raw data, and converting that right into an EList. Your $targets data.frame looks like it would just be the first three columns of your expdata data.frame you have already.
3. Shouldn't it rather be colnames(de_def) <- expdat$sample?

4. You might do better to include more of your data in the analysis, and only focus on the genes that you want in the final step. After I create the EList, I'd then only remove lowly expressed genes. I'd then do all the limma mojo on that and only subset down to genes you care about after your call to topTable. You can then readjust those p-values like so:

   res.all <- topTable(lmfit, "Th1_vs_Nai4", n = Inf)
   res <- res[rownames(res) %in% genes_I_care_about,]
   res$adj.P.Value <- p.adjust(res$P.Value, "BH")
   

   Having more data (genes) to work with is typically better for your eBayes step and should give you more (legitimate) power.

I think that's all I got for now ...