runDEAnalysis: Differential expression analysis

In RCPA: Consensus Pathway Analysis

Differential expression analysis

Description

This function performs differential expression analysis using either limma, DESeq2 or edgeR.

Usage

runDEAnalysis(
  summarizedExperiment,
  method = c("limma", "DESeq2", "edgeR"),
  design,
  contrast,
  annotation = NULL
)

Arguments

summarizedExperiment	SummarizedExperiment object
method	Method to use for differential expression analysis. Can be "limma", "DESeq2" or "edgeR".
design	A design model output by model.matrix.
contrast	A contrast matrix. See limma::makeContrasts.
annotation	A data frame mapping between probe IDs and entrez gene IDs. If not provided, the function will try to get the mapping from the platform annotation in the SummarizedExperiment object. If the annotation is not available, the function will return the probe IDs. Regardless of the type of annotation, it must contains two columns: FROM and TO, where FROM is the probe ID and TO is the entrez gene ID.

Value

A SummarizedExperiment object with DE analysis results appended to the rowData slot with the following columns:

• ID: gene ID. If annotation is provided, this will be the entrez gene ID. Otherwise, it will be the probe ID.

• logFC: log2 fold change

• p.value: p-value from the DE analysis using the specified method

• pFDR: p-value adjusted for multiple testing using Benjamini-Hochberg method

• statistic: statistic from the DE analysis using the specified method. For limma, this is the t-statistic. For DESeq2, this is the Wald statistic. For edgeR, this is the log fold change.

• avgExpr: For limma, it is the average expression. For DESeq2, it is the log base mean. For edgeR, it is the log CPM.

• logFCSE: standard error of the log fold change.

• sampleSize: sample size used for DE analysis.

The assay slot will contain the input expression/count matrix, and the rownames will be mapped to the gene IDs if annotation is found in the input SummarizedExperiment object or in the annotation parameter. Other slots will be the same as in the input SummarizedExperiment object.

Examples

library(RCPA)
library(SummarizedExperiment)

# GSE5281
affyDataset <- loadData("affyDataset")
affyDesign <- model.matrix(~0 + condition + region + condition:region, 
                             data = colData(affyDataset))
colnames(affyDesign) <- make.names(colnames(affyDesign))
affyContrast <- limma::makeContrasts(conditionalzheimer-conditionnormal, 
                                     levels=affyDesign)

if (require("hgu133plus2.db", quietly = TRUE)){
affyDEExperiment <- RCPA::runDEAnalysis(affyDataset, method = "limma", 
                                        design = affyDesign,
                                        contrast = affyContrast, 
                                        annotation = "GPL570")
}

# check the DE analysis results

print(head(rowData(affyDEExperiment)))


# GSE61196
agilDataset <- loadData("agilDataset")
agilDesign <- model.matrix(~0 + condition, 
                            data = colData(agilDataset))
agilContrast <- limma::makeContrasts(conditionalzheimer-conditionnormal, 
                                     levels=agilDesign)

# Create Probe mapping
GPL4133Anno <- GEOquery::dataTable(GEOquery::getGEO("GPL4133"))@table
GPL4133GeneMapping <- data.frame(FROM = GPL4133Anno$SPOT_ID, 
                                 TO = as.character(GPL4133Anno$GENE),
                                 stringsAsFactors = FALSE)
GPL4133GeneMapping <- GPL4133GeneMapping[!is.na(GPL4133GeneMapping$TO), ]

agilDEExperiment <- RCPA::runDEAnalysis(agilDataset, method = "limma", 
                                        design = agilDesign,
                                        contrast = agilContrast, 
                                        annotation = GPL4133GeneMapping)
print(head(rowData(agilDEExperiment)))

# GSE153873
RNASeqDataset <- loadData("RNASeqDataset")
RNASeqDesign <- model.matrix(~0 + condition, data = colData(RNASeqDataset))
RNASeqContrast <- limma::makeContrasts(conditionalzheimer-conditionnormal, 
                                      levels=RNASeqDesign)


if (require("org.Hs.eg.db", quietly = TRUE)){
    ENSEMBLMapping <- AnnotationDbi::select(org.Hs.eg.db, 
                                            keys = rownames(RNASeqDataset),
                                            columns = c("SYMBOL", "ENTREZID"), 
                                            keytype = "SYMBOL")
    colnames(ENSEMBLMapping) <- c("FROM", "TO")

    RNASeqDEExperiment <- RCPA::runDEAnalysis(RNASeqDataset,
                           method = "DESeq2",
                           design = RNASeqDesign,
                           contrast = RNASeqContrast,
                           annotation = ENSEMBLMapping)
    print(head(rowData(RNASeqDEExperiment)))
}

RCPA documentation built on Nov. 21, 2023, 5:08 p.m.