added extra plots and proteomics

ETHZ-INS · Jul 24, 2023 · 5832289 · 5832289
1 parent 09df7b8
commit 5832289
Show file tree

Hide file tree

Showing 5 changed files with 7,256 additions and 2 deletions.
diff --git a/01_DEA_exploration.Rmd b/01_DEA_exploration.Rmd
@@ -81,6 +81,7 @@ sechm(se, topDegs, assay="scaledLFC", gaps_at="condition", breaks=TRUE,
 ```
 
 ```{r, fig.width=8, fig.height=6}
+se <- readRDS("SE.processed.rds")
 rowData(se)$logFC.2hDEX <- rowData(se)$logFC.18hDEX <- NA 
 rowData(se)$propInhib.prior <- rowData(se)$propInhib.after <- NA 
 
@@ -97,11 +98,13 @@ d$propInhib <- -(1-2^abs(d$logFC-d$DEX.logFC))
 d$propInhib[d$propInhib>1] <- 1
 ag <- aggregate(d[,c("DEX.logFC","propInhib")], by=d[,"gene",drop=FALSE], FUN=median)
 row.names(ag) <- ag$gene
+ag <- ag[intersect(row.names(ag),row.names(se)),]
 rowData(se)[row.names(ag), "logFC.2hDEX"] <- ag$DEX.logFC
 rowData(se)[row.names(ag), "propInhib.before"] <- ag$propInhib
+d$condition <- factor(d$condition, c("16h MIF  > 2h DEX+MIF", "16h Cort113 > 2h DEX+Cort113", "16h KH-103 > 2h DEX+KH-103"))
 p1 <- ggplot(d, aes(order, logFC)) + geom_line(aes(order, DEX.logFC), lwd=1.5, colour="red") + 
   geom_hline(yintercept=0) + geom_point(alpha=0.2, size=0.5) + facet_wrap(~condition) + ylim(-3,3) +
-  xlab("logFC rank in 16h DMSO > 2h DEX+DMSO") + 
+  xlab("logFC rank in 16h DMSO > 2h DEX+DMSO") + geom_smooth() +
   scale_x_continuous(breaks=c(0,max(d$order)), labels=c("down", "up"))
 
 tmp <- as.list(rowData(se)[,grep("^DEA\\.2h",colnames(rowData(se)))])
@@ -117,13 +120,19 @@ d$propInhib <- -(1-2^abs(d$logFC-d$DEX.logFC))
 d$propInhib[d$propInhib>1] <- 1
 ag <- aggregate(d[,c("DEX.logFC","propInhib")], by=d[,"gene",drop=FALSE], FUN=median)
 row.names(ag) <- ag$gene
+ag <- ag[intersect(row.names(ag),row.names(se)),]
 rowData(se)[row.names(ag), "logFC.18hDEX"] <- ag$DEX.logFC
 rowData(se)[row.names(ag), "propInhib.after"] <- ag$propInhib
+d$condition <- factor(d$condition, c("2h DEX > 16h DEX+MIF", "2h DEX > 16h DEX+Cort113", "2h DEX > 16h DEX+KH-103"))
 p2 <- ggplot(d, aes(order, logFC)) + geom_line(aes(order, DEX.logFC), lwd=1.5, colour="red") + 
   geom_hline(yintercept=0) + geom_point(alpha=0.2, size=0.5) + facet_wrap(~condition) + ylim(-3,3) +
-  xlab("logFC rank in 2h DEX > 16h DEX+DMSO") + 
+  xlab("logFC rank in 2h DEX > 16h DEX+DMSO") + geom_smooth() +
   scale_x_continuous(breaks=c(0,max(d$order)), labels=c("down", "up"))
 
+pdf("inhibition_curves.pdf", width=8, height=6)
+plot_grid(p1, p2, nrow=2, labels="AUTO")
+dev.off()
+
 plot_grid(p1, p2, nrow=2, labels="AUTO")
 ```
 (The red line indicates the foldchanges upon DEX (2h in A, 18h in B), while the dots indicate the foldhanges with the respective inhibitor)

diff --git a/01_DEA_stringent.Rmd b/01_DEA_stringent.Rmd
@@ -0,0 +1,163 @@
+---
+title: "MG A549 experiments - DEA and exploration"
+author: "Pierre-Luc Germain"
+date: "2023/07/17"
+output:
+  html_document:
+    toc: true
+    toc_float: true
+    code_folding: hide
+---
+
+```{r}
+suppressPackageStartupMessages({
+  library(SummarizedExperiment)
+  library(SEtools)
+  library(sechm)
+  library(edgeR)
+  library(ggplot2)
+  library(grid)
+  library(cowplot)
+})
+theme_set(theme_classic())
+lfc.th <- 1
+```
+
+# Overview 
+
+```{r, fig.width=8, fig.height=6}
+se <- readRDS("salmon/geneLevel.SE.rds")
+dds <- calcNormFactors(DGEList(assay(se)))
+assays(se)$logcpm <- log1p(cpm(dds))
+dds <- dds[filterByExpr(dds,model.matrix(~se$condition),min.count=20),]
+se <- se[row.names(dds),]
+se <- svacor(se, ~condition, n.sv=2)
+se$condition <- relevel(factor(se$condition), "18h untreated")
+se$condition2 <- se$condition
+levels(se$condition2) <- gsub("^18h DMSO$|18h untreated","control",levels(se$condition2))
+se$isInhibited <- grepl("Cort113|KH-103|MIF", se$condition2)
+se$isDEX <- grepl("DEX", se$condition2)
+se$condType <- factor(paste0(as.integer(se$isDEX),as.integer(se$isInhibited)))
+levels(se$condType) <- c("control", "inhibited", "DEX", "DEX+inhibited")
+metadata(se)$default_view <- list( assay="scaledLFC", groupvar="condType", colvar="condition",
+                                   top_annotation=c("condType") )
+
+mm <- model.matrix(~SV1+SV2+condition2, data=as.data.frame(colData(se)))
+dds <- estimateDisp(dds,mm)
+conds <- grep("condition2",colnames(mm),value=TRUE)
+names(conds) <- gsub("condition2","",conds)
+fit <- glmFit(dds,mm)
+res <- lapply(conds, FUN=function(x){
+  y <- as.data.frame(topTags(glmLRT(fit,x), Inf))[,c(1,4,5)]
+  attr(y, "description") <- paste0(gsub("condition2","",x), "  vs.  DMSO & untreated")
+  y
+})
+for(f in names(res)){
+  rowData(se)[[paste0("DEA.",f)]] <-res[[f]][row.names(se),]
+}
+topDegs <- unique(unlist(lapply(res, FUN=function(x){
+  head(row.names(x)[which(x$FDR<0.01 & abs(x$logFC)>lfc.th)],15)
+})))
+levels(se$condition) <- gsub("> ",">\n", levels(se$condition))
+se$condition <- factor(se$condition, 
+                       unique(c("18h untreated","18h DMSO",levels(se$condition))))
+se <- log2FC(se, "logcpm", se$condition2=="control")
+```
+
+```{r, fig.width=8, fig.height=8}
+sechm(se, topDegs, assay="scaledLFC", gaps_at="condition", breaks=TRUE,
+      top_annotation="batch", row_title="Union of top DEGs", show_rownames=TRUE,
+      row_names_gp=gpar(fontsize=9), column_title_gp=gpar(fontsize=10),
+      column_title_rot=90)
+```
+
+```{r, fig.width=8, fig.height=6}
+rowData(se)$logFC.2hDEX <- rowData(se)$logFC.18hDEX <- NA 
+rowData(se)$propInhib.prior <- rowData(se)$propInhib.after <- NA 
+
+tmp <- as.list(rowData(se)[,grep("\\.16h",colnames(rowData(se)))])
+tmp <- tmp[grep("DMSO",names(tmp),invert=TRUE)]
+names(tmp) <- gsub("^DEA\\.","",names(tmp))
+dea <- rowData(se)[["DEA.16h DMSO > 2h DEX+DMSO"]]
+degs <- row.names(dea)[which(dea$FDR<0.05 & abs(dea$logFC)>lfc.th)]
+d <- dplyr::bind_rows(lapply(tmp, FUN=function(x) data.frame(gene=degs, logFC=x[degs,"logFC"])), .id="condition")
+o <- sort(rank(setNames(dea[degs,"logFC"], degs)))
+d$order <- o[d$gene]
+d$DEX.logFC <- dea[d$gene,"logFC"]
+d$propInhib <- -(1-2^abs(d$logFC-d$DEX.logFC))
+d$propInhib[d$propInhib>1] <- 1
+ag <- aggregate(d[,c("DEX.logFC","propInhib")], by=d[,"gene",drop=FALSE], FUN=median)
+row.names(ag) <- ag$gene
+ag <- ag[intersect(row.names(ag),row.names(se)),]
+rowData(se)[row.names(ag), "logFC.2hDEX"] <- ag$DEX.logFC
+rowData(se)[row.names(ag), "propInhib.before"] <- ag$propInhib
+d$condition <- factor(d$condition, c("16h MIF  > 2h DEX+MIF", "16h Cort113 > 2h DEX+Cort113", "16h KH-103 > 2h DEX+KH-103"))
+p1 <- ggplot(d, aes(order, logFC)) + geom_line(aes(order, DEX.logFC), lwd=1.5, colour="red") + 
+  geom_hline(yintercept=0) + geom_point(alpha=0.2, size=0.5) + facet_wrap(~condition) + ylim(-3,3) +
+  xlab("logFC rank in 16h DMSO > 2h DEX+DMSO") + geom_smooth() +
+  scale_x_continuous(breaks=c(0,max(d$order)), labels=c("down", "up"))
+
+tmp <- as.list(rowData(se)[,grep("^DEA\\.2h",colnames(rowData(se)))])
+tmp <- tmp[grep("DMSO",names(tmp),invert=TRUE)]
+names(tmp) <- gsub("^DEA\\.","",names(tmp))
+dea <- rowData(se)[["DEA.2h DEX > 16h DEX+DMSO"]]
+degs <- row.names(dea)[which(dea$FDR<0.05 & abs(dea$logFC)>lfc.th)]
+d <- dplyr::bind_rows(lapply(tmp, FUN=function(x) data.frame(gene=degs, logFC=x[degs,"logFC"])), .id="condition")
+o <- sort(rank(setNames(dea[degs,"logFC"], degs)))
+d$order <- o[d$gene]
+d$DEX.logFC <- dea[d$gene,"logFC"]
+d$propInhib <- -(1-2^abs(d$logFC-d$DEX.logFC))
+d$propInhib[d$propInhib>1] <- 1
+ag <- aggregate(d[,c("DEX.logFC","propInhib")], by=d[,"gene",drop=FALSE], FUN=median)
+row.names(ag) <- ag$gene
+ag <- ag[intersect(row.names(ag),row.names(se)),]
+rowData(se)[row.names(ag), "logFC.18hDEX"] <- ag$DEX.logFC
+rowData(se)[row.names(ag), "propInhib.after"] <- ag$propInhib
+d$condition <- factor(d$condition, c("2h DEX > 16h DEX+MIF", "2h DEX > 16h DEX+Cort113", "2h DEX > 16h DEX+KH-103"))
+p2 <- ggplot(d, aes(order, logFC)) + geom_line(aes(order, DEX.logFC), lwd=1.5, colour="red") + 
+  geom_hline(yintercept=0) + geom_point(alpha=0.2, size=0.5) + facet_wrap(~condition) + ylim(-3,3) +
+  xlab("logFC rank in 2h DEX > 16h DEX+DMSO") + geom_smooth() +
+  scale_x_continuous(breaks=c(0,max(d$order)), labels=c("down", "up"))
+plot_grid(p1, p2, nrow=2, labels="AUTO")
+```
+
+
+```{r, fig.width=8, fig.height=6}
+pdf("inhibition_curves.pdf", width=8, height=6)
+plot_grid(p1, p2, nrow=2, labels="AUTO")
+dev.off()
+```
+(The red line indicates the foldchanges upon DEX (2h in A, 18h in B), while the dots indicate the foldhanges with the respective inhibitor)
+KH-103 is *more efficient when administered before the activation*, but *less efficient when administered after*, where mifrepristone seems best.
+
+
+
+# Side effects of the inhibitors
+
+## Genes triggered by any blocker
+
+```{r, fig.width=5, fig.height=3.5}
+w <- se$condition2 %in% c("control","18h Cort113","18h KH-103","18h MIF")
+dds2 <- dds[,w]
+mm <- model.matrix(~droplevels(se$condition2[w]))
+fit <- glmFit(dds2,mm)
+res <- lapply(colnames(mm)[-1], FUN=function(x){
+  as.data.frame(topTags(glmLRT(fit, x),Inf))[,c(1,4,5)]
+})
+degs <- getDEGs(res,lfc.th=lfc.th, fdr.th=0.05,merge=TRUE)
+
+h <- sechm(se[,w], degs, assay="scaledLFC", gaps_at="condition", 
+      breaks=0.999, top_annotation="batch", row_title="Triggered by any blocker", 
+      column_title_gp=gpar(fontsize=10), column_title_rot=90, row_names_gp=gpar(fontsize=8))
+draw(h,merge=TRUE)
+```
+
+
+
+
+# Session info
+
+```{r}
+sessionInfo()
+```
+
diff --git a/03_MAplots.Rmd b/03_MAplots.Rmd
@@ -0,0 +1,166 @@
+---
+title: "MG A549 experiments - MA plots"
+author: "Pierre-Luc Germain"
+date: "2023/07/17"
+output:
+  html_document:
+    toc: true
+    toc_float: true
+    code_folding: hide
+---
+
+```{r}
+suppressPackageStartupMessages({
+  library(SummarizedExperiment)
+  library(sechm)
+  library(ggplot2)
+  library(ggrastr)
+  library(cowplot)
+  library(ComplexHeatmap)
+  library(ggvenn)
+  library(UpSetR)
+})
+theme_set(theme_classic())
+se <- readRDS("SE.processed.rds")
+```
+
+```{r}
+deas <- getDEA(se, homogenize = FALSE)
+d2h <- deas$`16h DMSO > 2h DEX+DMSO`
+d18h <- deas$`2h DEX > 16h DEX+DMSO`
+```
+
+```{r}
+masig <- function(dea, hl=NULL, ylim=NULL, FDR.th=0.01, lfc.th=log2(1.2)){
+  dea$logCPM <- rowMeans(assays(se)$logcpm)
+  if(is.null(ylim)) ylim <- range(dea$logFC[which(dea$logCPM>=1)],na.rm=TRUE)
+  print(ylim)
+  p <- ggplot(dea, aes(logCPM, logFC)) + geom_point_rast(aes(colour=-log10(FDR))) +
+    scale_color_viridis_c(trans="sqrt", option="B", direction=-1) + coord_cartesian(ylim=ylim)
+  if(is.null(hl) || (length(hl)==1 && is.integer(hl))){
+    if(is.null(hl)) hl <- 10
+    dea <- dea[which(dea$FDR<FDR.th & abs(dea$logFC)>=lfc.th),]
+    hl <- head(row.names(dea)[order(dea$FDR)],hl)
+  }
+  if(all(is.na(hl))) return(p)
+  dea <- dea[hl,]
+  dea$gene <- row.names(dea)
+  p + ggrepel::geom_text_repel(data=dea, aes(label=gene), min.segment.length=0, size=3,
+                               nudge_y=0.05*sign(dea$logFC)*abs(diff(ylim)))
+}
+```
+
+```{r, fig.width=9, fig.height=4}
+p1 <- masig(d2h) + ggtitle("16h DMSO > 2h DEX+DMSO vs 18h DMSO")
+p2 <- masig(d18h) + ggtitle("2h DEX > 16h DEX+DMSO vs 18h DMSO")
+leg <- cowplot::get_legend(p1)
+p <- cowplot::plot_grid(p1 + theme(legend.position="none"),
+                       p2 + theme(legend.position="none"),
+                       leg, rel_widths=c(4,4,1.5), nrow=1, labels=c("A","B",NA), scale=0.95)
+pdf("maPlots.pdf", width=9, height=4, pointsize = 10)
+p
+dev.off()
+```
+
+```{r, fig.width=9, fig.height=8}
+pl <- lapply(c("18h Cort113", "18h KH-103","18h MIF"), FUN=function(x){
+  masig(deas[[x]], ylim=c(-3.8,4.5)) + ggtitle(paste0(x, " vs 18h DMSO")) + theme(legend.position="bottom")
+})
+pp <- plot_grid(p, nrow = 2,
+          plot_grid(plotlist=pl, nrow=1, labels=c("C","D","E"), scale=0.95))
+pdf("maPlots2.pdf", width=10, height=8, pointsize = 10)
+pp
+dev.off()
+```
+
+# Venn Diagrams
+
+## Side effects
+
+```{r, fig.width=10, fig.height=6}
+degs <- getDEGs(deas[c(grep("18h",names(deas), value=TRUE), "2h DEX > 16h DEX+DMSO")], lfc.th=log2(1.2), merge=FALSE)
+names(degs)[length(degs)] <- "18h DEX"
+p1 <- ggvenn(degs, set_name_size = 4) + ggtitle("Side effects of the inhibitors")
+upset(fromList(degs))
+```
+
+## Inhibition
+
+```{r, fig.width=10, fig.height=6}
+degs <- getDEGs(deas[c(grep("^2h",names(deas), value=TRUE))], lfc.th=log2(1.2), merge=FALSE)
+names(degs) <- gsub("16h DEX+","",gsub("^2h DEX > ","DEX",names(degs)))
+degs <- c(degs[-2],degs[2])
+p2 <- ggvenn(degs, set_name_size = 4) + ggtitle("Inhibition after DEX")
+
+degs <- getDEGs(deas[c(grep("^16h",names(deas), value=TRUE))], lfc.th=log2(1.2), merge=FALSE)
+names(degs) <- gsub("16h ","",gsub(" > 2h DEX.+","+DEX",names(degs)))
+degs <- c(degs[-2],degs[2])
+p3 <- ggvenn(degs, set_name_size = 4) + ggtitle("Inhibition prior DEX")
+
+```
+
+```{r}
+pdf("venn.pdf", width=10, height=18)
+plot_grid(p1,p2,p3, nrow=3)
+dev.off()
+```
+```{r}
+
+```
+
+## Icaros
+
+```{r, fig.height=4, fig.width=10}
+# known targets in b cells
+ic <- c("CCND2", "CDKN1A", "CDK6", "MYC", "LIG4", "FOXO1", "CD79B", "BLNK", "SYKB", "VPREB1")
+se$condition
+sechm(se, ic, gaps_at="condition", column_title_rot=90, row_title="Ikaros targets")
+sechm(se, grep("^IKZF",row.names(se)), gaps_at="condition", column_title_rot=90)
+```
+
+
+```{r, fig.height=7, fig.width=7}
+ggplot(sechm::meltSE(se[,grep("^16h|^18h", se$condition)], intersect(ic,row.names(se)), rowDat.columns = NA), aes(condition, log2FC)) + geom_hline(yintercept=0, linetype="dashed") + geom_boxplot() + facet_wrap(~feature, scales = "free_x") + coord_flip()
+```
+
+
+```{r}
+ik <- readRDS("IKZF1.rds")
+ik <- ik[ik$score>8]
+library(AnnotationHub)
+ah <- AnnotationHub(localHub=TRUE)
+ensdb <- ah[["AH109336"]]
+tx <- transcripts(ensdb, columns=c("tx_biotype", "gene_name", "tx_support_level"))
+tx <- tx[which(tx$tx_biotype=="protein_coding" & tx$tx_support_level==1L)]
+proms <- promoters(tx, upstream=1000, downstream=200)
+proms$ikaros <- overlapsAny(proms, ik)
+gsets <- lapply(split(proms$gene_name, proms$ikaros), unique)
+names(gsets) <- c("none","Ikaros")
+gsets$none <- setdiff(gsets$none,gsets$Ikaros)
+gsets <- lapply(gsets, intersect, y=row.names(se))
+degs <- getDEGs(deas[c(grep("18h",names(deas), value=TRUE), "2h DEX > 16h DEX+DMSO")], lfc.th=log2(1.2), merge=FALSE)
+names(degs)[length(degs)] <- "18h DEX"
+degs$Ikaros <- gsets$Ikaros
+ggvenn(degs[-4], set_name_size = 4)
+
+overlap.prob <- function (set1, set2, universe, lower = F){
+    set1 <- unique(as.character(set1))
+    set2 <- unique(as.character(set2))
+    if (class(universe) == "character") {
+        set1 <- intersect(set1, universe)
+        set2 <- intersect(set2, universe)
+        universe <- length(unique(universe))
+    }
+    phyper(max(0, sum(set1 %in% set2)-1), length(set1), universe - length(set1), 
+        length(set2), lower.tail = lower)
+}
+sapply(degs, FUN=function(x) overlap.prob(x, gsets$Ikaros, unlist(gsets)))
+```
+
+
+# Session info
+
+```{r}
+sessionInfo()
+```
+