inst/doc/cTRAP.R

In cTRAP: Identification of candidate causal perturbations from differential gene expression data

## ----load package-------------------------------------------------------------
library(cTRAP)

## ----load ENCODE samples, eval=FALSE------------------------------------------
#  gene <- "EIF4G1"
#  cellLine <- "HepG2"
#  
#  ENCODEmetadata <- downloadENCODEknockdownMetadata(cellLine, gene)
#  table(ENCODEmetadata$`Experiment target`)
#  length(unique(ENCODEmetadata$`Experiment target`))
#  
#  ENCODEsamples <- loadENCODEsamples(ENCODEmetadata)[[1]]
#  counts <- prepareENCODEgeneExpression(ENCODEsamples)

## ----load data, include=FALSE-------------------------------------------------
data("ENCODEmetadata")
data("counts")

## ----differential gene expression, eval=FALSE---------------------------------
#  # Remove low coverage (at least 10 counts shared across two samples)
#  minReads   <- 10
#  minSamples <- 2
#  filter <- rowSums(counts[ , -c(1, 2)] >= minReads) >= minSamples
#  counts <- counts[filter, ]
#  
#  # Convert ENSEMBL identifier to gene symbol
#  counts$gene_id <- convertENSEMBLtoGeneSymbols(counts$gene_id)
#  
#  # Perform differential gene expression analysis
#  diffExpr <- performDifferentialExpression(counts)

## ----t-statistics, eval=FALSE-------------------------------------------------
#  # Get t-statistics of differential expression with respective gene names
#  # (expected input for downstream analyses)
#  diffExprStat <- diffExpr$t
#  names(diffExprStat) <- diffExpr$Gene_symbol

## ---- echo=FALSE--------------------------------------------------------------
data("diffExprStat")

## ----CMap metadata conditions-------------------------------------------------
# Load CMap metadata (automatically downloaded if not found)
cmapMetadata <- loadCMapData("cmapMetadata.txt", type="metadata")

# Summarise conditions for all CMap perturbations
getCMapConditions(cmapMetadata)

# Summarise conditions for CMap perturbations in HepG2 cell line
getCMapConditions(cmapMetadata, cellLine="HepG2")

# Summarise conditions for a specific CMap perturbation in HepG2 cell line
getCMapConditions(
    cmapMetadata, cellLine="HepG2",
    perturbationType="Consensus signature from shRNAs targeting the same gene")

## ----CMap knockdown perturbations, eval=FALSE---------------------------------
#  # Filter CMap gene knockdown HepG2 data to be loaded
#  cmapMetadataKD <- filterCMapMetadata(
#      cmapMetadata, cellLine="HepG2",
#      perturbationType="Consensus signature from shRNAs targeting the same gene")
#  
#  # Load filtered data (data will be downloaded if not found)
#  cmapPerturbationsKD <- prepareCMapPerturbations(
#      metadata=cmapMetadataKD, zscores="cmapZscores.gctx",
#      geneInfo="cmapGeneInfo.txt")

## ----CMap small molecule perturbations, eval=FALSE----------------------------
#  # Filter CMap gene small molecule HepG2 data to be loaded
#  cmapMetadataCompounds <- filterCMapMetadata(
#      cmapMetadata, cellLine="HepG2", perturbationType="Compound")
#  
#  # Load filtered data (data will be downloaded if not found)
#  cmapPerturbationsCompounds <- prepareCMapPerturbations(
#      metadata=cmapMetadataCompounds, zscores="cmapZscores.gctx",
#      geneInfo="cmapGeneInfo.txt", compoundInfo="cmapCompoundInfo.txt")

## ----load CMap perturbations, include=FALSE-----------------------------------
data("cmapPerturbationsKD")
data("cmapPerturbationsCompounds")
cmapPerturbationsCompounds <- cmapPerturbationsCompounds[
    , grep("HEPG2", colnames(cmapPerturbationsCompounds))]

## ----compare knockdown, message=FALSE-----------------------------------------
compareKD <- rankSimilarPerturbations(diffExprStat, cmapPerturbationsKD)

## ----compare small molecules, message=FALSE-----------------------------------
compareCompounds <- rankSimilarPerturbations(diffExprStat, 
                                             cmapPerturbationsCompounds)

## ----order knockdown perturbation comparison, fig.width=6, fig.asp=1----------
# Most positively associated perturbations (note that EIF4G1 knockdown is the
# 7th, 1st and 2nd most positively associated perturbation based on Spearman's
# correlation, Pearson's correlation and GSEA, respectively)
head(compareKD[order(spearman_rank)], n=10)
head(compareKD[order(pearson_rank)])
head(compareKD[order(GSEA_rank)])
head(compareKD[order(rankProduct_rank)])

# Most negatively associated perturbations
head(compareKD[order(-spearman_rank)])
head(compareKD[order(-pearson_rank)])
head(compareKD[order(-GSEA_rank)])
head(compareKD[order(-rankProduct_rank)])

# Plot list of compared pertubations
plot(compareKD, "spearman", n=c(10, 3))
plot(compareKD, "pearson")
plot(compareKD, "gsea")
plot(compareKD, "rankProduct")

## ----order small molecule perturbation comparison, fig.width=6, fig.asp=1-----
# Most positively and negatively associated perturbations
compareCompounds[order(rankProduct_rank)]
plot(compareCompounds, "rankProduct")

## ----perturbation info, fig.width=5, fig.asp=1--------------------------------
# Information on the EIF4G1 knockdown perturbation
EIF4G1knockdown <- grep("EIF4G1", compareKD[[1]], value=TRUE)
print(compareKD, EIF4G1knockdown)

# Information on the top 10 most similar compound perturbations (based on
# Spearman's correlation)
print(compareKD[order(rankProduct_rank)], 1:10)

# Get table with all information available (including ranks, metadata, compound
# information, etc.)
table <- as.table(compareKD)
head(table)

# Obtain available raw information from compared perturbations
names(attributes(compareKD)) # Data available in compared perturbations
attr(compareKD, "metadata")  # Perturbation metadata
attr(compareKD, "geneInfo")  # Gene information

## ---- echo=FALSE--------------------------------------------------------------
attr(compareKD, "zscoresFilename") <- cmapPerturbationsKD

## ----plot knockdown comparison, message=FALSE, fig.width=5, fig.asp=1---------
# Plot relationship with EIF4G1 knockdown from CMap
plot(compareKD, EIF4G1knockdown, "spearman")
plot(compareKD, EIF4G1knockdown, "pearson")
plot(compareKD, EIF4G1knockdown, "gsea")

# Plot relationship with most negatively associated perturbation
plot(compareKD, compareKD[order(-spearman_rank)][1, 1], "spearman")
plot(compareKD, compareKD[order(-pearson_rank)][1, 1], "pearson")
plot(compareKD, compareKD[order(-GSEA_rank)][1, 1], "gsea")

## ---- echo=FALSE--------------------------------------------------------------
attr(compareCompounds, "zscoresFilename") <- cmapPerturbationsCompounds

## ----plot small molecule comparison, message=FALSE, fig.width=5, fig.asp=1----
# Plot relationship with most positively associated perturbation
plot(compareCompounds, compareCompounds[order(spearman_rank)][1, 1], "spearman")
plot(compareCompounds, compareCompounds[order(pearson_rank)][1, 1], "pearson")
plot(compareCompounds, compareCompounds[order(GSEA_rank)][1, 1], "gsea")

# Plot relationship with most negatively associated perturbation
plot(compareCompounds, compareCompounds[order(-spearman_rank)][1,1], "spearman")
plot(compareCompounds, compareCompounds[order(-pearson_rank)][1, 1], "pearson")
plot(compareCompounds, compareCompounds[order(-GSEA_rank)][1, 1], "gsea")

## ---- fig.width=5, fig.asp=1--------------------------------------------------
listExpressionDrugSensitivityAssociation()
ctrp      <- listExpressionDrugSensitivityAssociation()[[2]]
assoc     <- loadExpressionDrugSensitivityAssociation(ctrp)
predicted <- predictTargetingDrugs(diffExprStat, assoc)
plot(predicted, method="rankProduct")

# Plot results for a given drug
plot(predicted, predicted[[1, 1]], method="spearman")
plot(predicted, predicted[[1, 1]], method="gsea")

## ---- fig.width=5, fig.asp=1--------------------------------------------------
# Label by compound name
plotTargetingDrugsVSsimilarPerturbations(
  predicted, compareCompounds, column="spearman_rank")
# Label by compound's gene target
plotTargetingDrugsVSsimilarPerturbations(
  predicted, compareCompounds, column="spearman_rank", labelBy="target")
# Label by compound's mechanism of action (MOA)
plotTargetingDrugsVSsimilarPerturbations(
  predicted, compareCompounds, column="spearman_rank", labelBy="moa")

## ----drug set enrichment analysis, fig.width=5, fig.asp=1---------------------
descriptors <- loadDrugDescriptors("CMap", "2D")
drugSets    <- prepareDrugSets(descriptors)
dsea        <- analyseDrugSetEnrichment(drugSets, predicted)
# Plot the 6 most significant drug set enrichment results
plotDrugSetEnrichment(drugSets, predicted, selectedSets=head(dsea$pathway, 6))