Limma for DEG analysis

Question

[rahma.khoualdia]

Hello everyone ! this is my first time using the limma package for running DEG analysis and I'm a bit lost cause when I compare my results with FindMarkers() function I don't get the same genes and logFC is so low comparing to the other method , i'll share the code with you if u can detect any modification needed. thank u !

### LimmaR
library(limma)
library(Seurat)
library(dplyr)
expression_matrix <- GetAssayData(seu_dblt_filt0_filtered1, assay = "RNA", slot = "data")
expression_matrix <- as.matrix(expression_matrix)
gene_variances <- apply(expression_matrix, 1, var)
zero_variance_genes <- which(gene_variances == 0)
expression_matrix_filtered <- expression_matrix[-zero_variance_genes, ]
design <- model.matrix(~0 + seu_dblt_filt0_filtered1$condition)
colnames(design) <- gsub("seu_dblt_filt0_filtered1$condition", "", colnames(design))
fit <- lmFit(expression_matrix_filtered, design)
fit2 <- eBayes(fit)
result <- topTable(fit2, adjust = "fdr", number = Inf)

# Rename columns in the design matrix to be syntactically valid
colnames(design) <- gsub("seu_dblt_filt0_filtered1\\$condition", "", colnames(design))  #
colnames(design) <- gsub(" ", "_", colnames(design))  
colnames(design) <- make.names(colnames(design))  
# Check the new column names
colnames(design)
# Define the contrast for FAP_injured vs FAP_non_injured
contrast_matrix_fap <- makeContrasts(FAP_injured - FAP_non.injured, levels = design)

# Apply the contrast to the linear model
fit_fap <- contrasts.fit(fit2, contrast_matrix_fap)
fit_fap <- eBayes(fit_fap, trend=TRUE)
contrast_results_fap <- topTable(fit_fap, adjust = "fdr", number = Inf)

# View FAP contrast results
head(contrast_results_fap)

# Define the contrast for ASC_injured vs ASC_non_injured
contrast_matrix_asc <- makeContrasts(ASC_injured - ASC_non.injured, levels = design)

# Apply the contrast to the linear model
fit_asc <- contrasts.fit(fit2, contrast_matrix_asc)
fit_asc <- eBayes(fit_asc)
contrast_results_asc <- topTable(fit_asc, adjust = "fdr", number = Inf)

# View ASC contrast results
head(contrast_results_asc)

# Add the gene names as a new column
contrast_results_asc$gene <- rownames(contrast_results_asc)

# Reorder the columns to make sure they appear in the desired order: gene, logFC, AveExpr, t, P.Value, adj.P.Val, B
contrast_results_asc <- contrast_results_asc %>%
  select(gene, logFC, AveExpr, t, P.Value, adj.P.Val, B)

# View the updated results
head(contrast_results_asc)

# Save to CSV
write.csv(contrast_results_asc, file = "contrast_results_asc.csv", row.names = FALSE)

# Add the gene names as a new column
contrast_results_fap$gene <- rownames(contrast_results_fap)

# Reorder the columns to make sure they appear in the desired order: gene, logFC, AveExpr, t, P.Value, adj.P.Val, B
contrast_results_fap <- contrast_results_fap %>%
  select(gene, logFC, AveExpr, t, P.Value, adj.P.Val, B)

# View the updated results
head(contrast_results_fap)

# Save to CSV
write.csv(contrast_results_fap, file = "contrast_results_fap.csv", row.names = FALSE)

[ATpoint]

Posting a wall of code is best practice. The limma user guide contains the best practice for the package. What I will say is that you should be sure to use proper input data which is usually for limma-trend the log2-transformed normalized counts or logCPMs. Are you using normalized data on log2 scale?

As for the results, Seurat is by default very primitive, meaning that logFCs are "raw" meaning not moderated for noise or variance and the pvalues are a simple WIlcox test that might (will) give different results to limma.

[rahma.khoualdia]

Thank you so much for your response ! i'm gonna check the structure of my data again and see the limma documentation