DEoverlap_FDR_forPaper.Rmd

---
title: "Differential Gene Expression overlap"
---
  
```{r setup, echo=F,message=F,warning=F}
library(cmapR)
library(colorspace)
library(scales)
source("~/Dropbox/GDB/line2user.R")
knitr::opts_chunk$set(
  echo=F,
  message=F,
  warning=F
)
```


```{r load_data}
load("~/Dropbox/GDB_archive/CMapCorr_files/lvl4_inputs_allgenes.RData")
lvl4_data <- lvl4_data_all
rm(lvl4_data_all)

temp_ligcountsperct <- sapply(unique(lvl4_data@cdesc$cell_id),function(CT)
  length(unique(lvl4_data@cdesc[lvl4_data@cdesc$cell_id == CT,"pert_iname"])))
ct9 <- names(temp_ligcountsperct)[temp_ligcountsperct > 100]
names(ct9) <- sapply(ct9,function(X) names(ct14)[ct14 == X])

temp_ctcountsperlig <- sapply(unique(lvl4_data@cdesc$pert_iname),function(LIG)
  length(unique(lvl4_data@cdesc[lvl4_data@cdesc$pert_iname == LIG,"cell_id"])))

lig295 <- names(temp_ctcountsperlig)[temp_ctcountsperlig >= 9]

rm(list=grep("^temp",ls(),value=T))
KEEP <- c("KEEP",ls())
```

# Figure 1

```{r fig1, fig.height=10,fig.width=9}
temp <- load("~/Dropbox/GDB_archive/CMapCorr_files/lig295_DE_lig_FDR.RData")
Z_all <- meanZ_lig
FDR_all <- FDR_lig
pDE_all <- pDE_lig
data_name <- "Ligands"
aCut <- 0.05

ALPHA <- c(0.1,0.05,0.01,0.001)
names(ALPHA) <- colnames(pDE_all)

# trim to genes with <Zp% probability of change by chance in at least one ligand
if (is.na(aCut)) {
  temp_cut <- ALPHA[1]
} else {
  temp_cut <- aCut
}
FDR <- FDR_all[apply(FDR_all,1,function(X) any(X <= temp_cut)),,drop=F]
if (ncol(FDR) >= 50) {
  FDR <- FDR[,apply(FDR,2,function(X) any(X <= temp_cut)),drop=F]
}
if (nrow(FDR) < 2 | ncol(FDR) < 2) {
  stop(paste0("Less than two samples reached FDR threshold of ",aCut * 100,"%"))
}

# order ligands
temp_hROW <- hclust(dist(t(FDR)),method="ward.D2") 
#### Euclidean distance in -log10 space? ####

# order genes
temp_hGENE <- hclust(dist(FDR),method="ward.D2")

# order genes by DE prior
# all(rownames(FDR) %in% DEprior$Gene_EntrezID)
# FDR <- FDR[DEprior$Gene_EntrezID[DEprior$Gene_EntrezID %in% rownames(FDR)],]

FDR <- FDR[temp_hGENE$order,temp_hROW$order]
countDE <- t(sapply(ALPHA,function(Z) apply(FDR,2,function(X) sum(X <= Z))))
FDR <- -log10(FDR) * sign(Z_all[rownames(FDR),colnames(FDR)])
pDE <- pDE_all[colnames(FDR),]
scoreDE <- pDE + 1e-4
scoreDE[scoreDE > 1] <- 1
scoreDE <- t(-log10(scoreDE))


layout(rbind(1:2),widths=c(6,1))
if (nrow(FDR) >= 100) {
  par(mar=c(4.5,6,1,0.5),mgp=2:0)
  label_cex <- 0.9
} else if (data_name == "Ligands") {
  par(mar=c(4.5,6,1,0.5),mgp=2:0)
  label_cex <- 1
} else if (data_name == "Cell lines") {
  par(mar=c(4.5,10,1.5,0.5),mgp=2:0)
  label_cex <- 0.9
} else {
  stop("data_name must be 'Ligands' or 'Cell lines'")
}
image(z=FDR / max(abs(FDR)),
      x=1:nrow(FDR),y=1:ncol(FDR),
      col=diverging_hcl(1000,palette="Blue-Red3"),
      breaks=seq(-1,1,length.out=1001),
      xaxt="n",yaxt="n",xlab=NA,ylab=NA)
if (nrow(FDR) < 100) {
  mtext(lvl4_data@rdesc[rownames(FDR),"pr_gene_symbol"],
        side=1,las=2,at=1:nrow(FDR),line=0.1,cex=0.7)
  mtext("Genes",side=1,line=3.5,font=2,cex=1.5)
} else {
  mtext("Genes",side=1,line=1,font=2,cex=1.5)
}
if (ncol(FDR) < 100) {
  mtext(colnames(FDR),side=2,las=2,
        at=1:ncol(FDR),adj=1,line=0.1,cex=label_cex)
}
if (nrow(FDR) >= 100) {
  mtext(data_name,side=2,line=2,font=2,cex=1.5,las=0)
} else if (data_name == "Ligands") {
  mtext(data_name,side=2,line=4.5,font=2,cex=1.5,las=0)
} else if (data_name == "Cell lines") {
  mtext(data_name,side=3,line=0,font=2,cex=1.5,las=0,
        at=par("usr")[1],adj=1)
} 

segments(x0=seq(line2user(5,2),line2user(1,2),length.out=1001),
         x1=seq(line2user(5,2),line2user(1,2),length.out=1001),
         y0=rep(line2user(1,1),1001),y1=rep(line2user(1.7,1),1001),
         xpd=NA,col=diverging_hcl(1001,palette="Blue-Red3"))
text(c(line2user(5,2),line2user(3,2),line2user(1,2)),
     rep(line2user(2.2,1),3),
     labels=c(-1 * round(max(abs(FDR))),
              0,
              round(max(abs(FDR)))),
     xpd=NA)
text(line2user(3,2),line2user(3.5,1),xpd=NA,
     labels=paste("Signed","-log10 FDR",sep="\n"))

if (!is.na(aCut)) {
  rect(xleft=which(abs(FDR) >= -log10(aCut),arr.ind=T)[,1] - 0.5,
       xright=which(abs(FDR) >= -log10(aCut),arr.ind=T)[,1] + 0.5,
       ybottom=which(abs(FDR) >= -log10(aCut),arr.ind=T)[,2] - 0.5,
       ytop=which(abs(FDR) >= -log10(aCut),arr.ind=T)[,2] + 0.5,
       border="mediumseagreen",lwd=0.5)
  rect(xleft=seq(line2user(3,2),line2user(1,2),
                 length.out=500)[round((-log10(aCut) / max(abs(FDR))) * 1001 / 2)],
       xright=line2user(1,2),ybottom=line2user(1.7,1),ytop=line2user(1,1),
       xpd=NA,col=NA,border="mediumseagreen")
  rect(xright=seq(line2user(5,2),line2user(3,2),
                  length.out=500)[500 - round((-log10(aCut) / max(abs(FDR))) * 1001 / 2)],
       xleft=line2user(5,2),ybottom=line2user(1.7,1),ytop=line2user(1,1),
       xpd=NA,col=NA,border="mediumseagreen")
  text(line2user(3,2),line2user(0.5,1),adj=0.5,xpd=NA,
       labels=paste0("FDR \u2264 ",aCut * 100,"%"),col="mediumseagreen",cex=0.9)  
}

if (nrow(FDR) >= 100) {
  par(mar=c(4.5,0.5,1,1),mgp=2:0)
} else if (data_name == "Ligands") {
  par(mar=c(4.5,0.5,1,1),mgp=2:0)
} else if (data_name == "Cell lines") {
  par(mar=c(4.5,0.5,1.5,1),mgp=2:0)
}
image(z=scoreDE / max(scoreDE),
      x=1:nrow(scoreDE),y=1:ncol(scoreDE),
      col=sequential_hcl(99,palette="PinkYl",rev=T),
      breaks=seq(0,1,length.out=100),
      xaxt="n",yaxt="n",xlab=NA,ylab=NA)
abline(v=seq(1.5,by=1,length.out=3),col="grey50")
mtext(c("10.0","5.0","1.0","0.1"),
      side=1,at=1:4,las=2,line=0.1,cex=0.9)
mtext("FDR %",side=1,at=0.2,las=2,line=0.1)
if (ncol(FDR) < 100) {
  text(as.vector(sapply(1:ncol(countDE),function(X) 1:nrow(countDE))),
       as.vector(sapply(1:ncol(countDE),function(X) rep(X,nrow(countDE)))),
       labels=as.vector(countDE),cex=0.8,col="grey20")
}
segments(x0=seq(1.5,4,length.out=1001),
         x1=seq(1.5,4,length.out=1001),
         y0=rep(line2user(2.1,1),1001),y1=rep(line2user(2.6,1),1001),
         xpd=NA,col=sequential_hcl(1001,palette="PinkYl",rev=T))
text(x=c(1.2,4.3),y=rep(line2user(2.4,1),2),xpd=NA,
     labels=c(0,max(scoreDE)))
mtext("-log10(P)",side=1,line=2.7,at=2.75)
mtext("# DE at FDR thresholds",side=4)


rm(list=ls()[!ls() %in% KEEP])
```


# Figure 2

```{r Fig2a, fig.height=4,fig.width=9}
temp <- load("~/Dropbox/GDB_archive/CMapCorr_files/lig295_DE_ligct_FDR.RData")
temp_FDRcut <- 0.05

pDEvector <- pDE_ligct[,as.character(temp_FDRcut * 100)]
names(pDEvector) <- sub("large_","",names(pDEvector))

scoreDE <- pDEvector + 1e-4
scoreDE[scoreDE > 1] <- 1
scoreDE <- -log10(scoreDE)

scoreMAT <- tapply(scoreDE,as.factor(sub("^.+_","",names(scoreDE))),c)
for (X in names(scoreMAT)) {
  names(scoreMAT[[X]]) <- sub(paste0("_",X),"",names(scoreMAT[[X]]))
  scoreMAT[[X]] <- scoreMAT[[X]][sort(names(scoreMAT[[X]]))]
}
rm(X)
scoreMAT <- t(do.call(rbind,scoreMAT))

if (ncol(scoreMAT) > 100) {
  scoreMAT <- scoreMAT[,apply(scoreMAT,2,function(X) any(X >= -log10(0.05 + 1e-4)))]
}

temp_hLIG <- hclust(dist(t(scoreMAT)),method="ward.D2")
temp_hCT <- hclust(dist(scoreMAT),method="ward.D2")
scoreMAT <- scoreMAT[temp_hCT$order,temp_hLIG$order]

par(mar=c(4.5,1,2,5.5),mgp=2:0)
image(t(scoreMAT),
      col=sequential_hcl(100,palette="PinkYl",rev=T),
      x=1:ncol(scoreMAT),y=1:nrow(scoreMAT),
      xaxt="n",yaxt="n",xlab=NA,ylab=NA)
mtext(paste0("Ligands (",ncol(scoreMAT)," / ",length(scoreDE) / nrow(scoreMAT),")"),
      side=1,line=0.5,font=1,cex=1.2)

text(x=rep(line2user(0.3,4),nrow(scoreMAT)),y=1:nrow(scoreMAT),
     labels=rownames(scoreMAT),xpd=NA,adj=0,srt=315,cex=0.8)
mtext("Cell lines",side=4,las=2,line=0.5,at=line2user(0,3),cex=1.2)

rect(xleft=which(t(scoreMAT) >= -log10(0.05 + 1e-4),arr.ind=T)[,1] - 0.5,
     xright=which(t(scoreMAT) >= -log10(0.05 + 1e-4),arr.ind=T)[,1] + 0.5,
     ybottom=which(t(scoreMAT) >= -log10(0.05 + 1e-4),arr.ind=T)[,2] - 0.5,
     ytop=which(t(scoreMAT) >= -log10(0.05 + 1e-4),arr.ind=T)[,2] + 0.5,
     border="dodgerblue")

temp_left <- par("usr")[1] + (par("usr")[2] - par("usr")[1]) * 0.2
temp_right <- par("usr")[1] + (par("usr")[2] - par("usr")[1]) * 0.8

segments(x0=seq(temp_left,temp_right,length.out=1000),
         x1=seq(temp_left,temp_right,length.out=1000),
         y0=rep(line2user(2,1),1000),y1=rep(line2user(2.6,1),1000),
         xpd=NA,col=sequential_hcl(1000,palette="PinkYl",rev=T))
mtext(c("100%","< 0.01%"),side=1,line=1.8,las=0,adj=c(1.1,-0.1),
      at=c(temp_left,temp_right))
mtext(paste0("Probability of at least # DE at ",
             as.character(temp_FDRcut * 100),
             "% FDR occuring by chance"),
      las=0,side=1,line=2.7,adj=0.5,
      at=temp_left + (temp_right - temp_left) / 2)

rect(xleft=temp_left + (-log10(0.05 + 1e-4) / max(scoreMAT)) * (temp_right - temp_left),
     xright=temp_right,ybottom=line2user(2,1),ytop=line2user(2.6,1),
     xpd=NA,col=NA,border="dodgerblue")
mtext("p \u2264 0.05",side=1,las=0,col="dodgerblue",cex=0.9,line=1,adj=1,at=temp_right)

mtext("A",font=2,cex=1.5,line=0.5,side=3,at=line2user(0.5,2))

rm(list=ls()[!ls() %in% KEEP])
```

Version 2, ligands with > 1 signif DE incl, # DE added.

```{r Fig2a_v2, fig.height=4,fig.width=9}
temp <- load("~/Dropbox/GDB_archive/CMapCorr_files/lig295_DE_ligct_FDR.RData")
temp_FDRcut <- 0.05

countDEvector <- apply(FDR_ligct,2,function(X) sum(X <= temp_FDRcut))
pDEvector <- pDE_ligct[,as.character(temp_FDRcut * 100)]
names(countDEvector) <- names(pDEvector) <- sub("large_","",names(pDEvector))

scoreDE <- pDEvector + 1e-4
scoreDE[scoreDE > 1] <- 1
scoreDE <- -log10(scoreDE)

scoreMAT <- tapply(scoreDE,as.factor(sub("^.+_","",names(scoreDE))),c)
countMAT <- tapply(countDEvector,as.factor(sub("^.+_","",names(countDEvector))),c)
for (X in names(scoreMAT)) {
  names(scoreMAT[[X]]) <- sub(paste0("_",X),"",names(scoreMAT[[X]]))
  names(countMAT[[X]]) <- sub(paste0("_",X),"",names(countMAT[[X]]))
  scoreMAT[[X]] <- scoreMAT[[X]][sort(names(scoreMAT[[X]]))]
  countMAT[[X]] <- countMAT[[X]][sort(names(countMAT[[X]]))]
}
rm(X)
scoreMAT <- t(do.call(rbind,scoreMAT))
countMAT <- t(do.call(rbind,countMAT))

if (ncol(scoreMAT) > 100) {
  temp_filt <- apply(scoreMAT,2,function(X) sum(X >= -log10(0.05 + 1e-4)) > 1)
  scoreMAT <- scoreMAT[,temp_filt]
  countMAT <- countMAT[,temp_filt]
}

temp_hLIG <- hclust(dist(t(scoreMAT)),method="ward.D2")
temp_hCT <- hclust(dist(scoreMAT),method="ward.D2")
scoreMAT <- scoreMAT[temp_hCT$order,temp_hLIG$order]
countMAT <- countMAT[temp_hCT$order,temp_hLIG$order]


par(mar=c(7,1,2,5.5),mgp=2:0)
image(t(scoreMAT),
      col=sequential_hcl(100,palette="PinkYl",rev=T),
      x=1:ncol(scoreMAT),y=1:nrow(scoreMAT),
      xaxt="n",yaxt="n",xlab=NA,ylab=NA)

text(x=rep(1:ncol(countMAT),nrow(countMAT)),
     y=as.integer(sapply(1:nrow(countMAT),rep,ncol(countMAT))),
     labels=as.integer(t(countMAT)),cex=0.8)

text(x=1:ncol(scoreMAT),y=rep(line2user(0.3,1),ncol(scoreMAT)),
     labels=colnames(scoreMAT),xpd=NA,adj=0,srt=315,cex=0.8)
mtext(paste0("Ligands (",ncol(scoreMAT)," / ",length(scoreDE) / nrow(scoreMAT),
             " with significant # DE in >1 cell line)"),
      side=1,line=2.5,font=1,cex=1.2)

text(x=rep(line2user(0.3,4),nrow(scoreMAT)),y=1:nrow(scoreMAT),
     labels=rownames(scoreMAT),xpd=NA,adj=0,srt=315,cex=0.8)
mtext("Cell lines",side=4,las=2,line=0.5,at=line2user(0.5,3),cex=1.2)

rect(xleft=which(t(scoreMAT) >= -log10(0.05 + 1e-4),arr.ind=T)[,1] - 0.5,
     xright=which(t(scoreMAT) >= -log10(0.05 + 1e-4),arr.ind=T)[,1] + 0.5,
     ybottom=which(t(scoreMAT) >= -log10(0.05 + 1e-4),arr.ind=T)[,2] - 0.5,
     ytop=which(t(scoreMAT) >= -log10(0.05 + 1e-4),arr.ind=T)[,2] + 0.5,
     border="dodgerblue")

temp_left <- par("usr")[1] + (par("usr")[2] - par("usr")[1]) * 0.2
temp_right <- par("usr")[1] + (par("usr")[2] - par("usr")[1]) * 0.8

segments(x0=seq(temp_left,temp_right,length.out=1000),
         x1=seq(temp_left,temp_right,length.out=1000),
         y0=rep(line2user(4.5,1),1000),y1=rep(line2user(5.1,1),1000),
         xpd=NA,col=sequential_hcl(1000,palette="PinkYl",rev=T))
mtext(c("100%","< 0.01%"),side=1,line=4.3,las=0,adj=c(1.1,-0.1),
      at=c(temp_left,temp_right))
mtext(paste0("Probability of at least # DE at ",
             as.character(temp_FDRcut * 100),
             "% FDR occuring by chance"),
      las=0,side=1,line=5.2,adj=0.5,
      at=temp_left + (temp_right - temp_left) / 2)

rect(xleft=temp_left + (-log10(0.05 + 1e-4) / max(scoreMAT)) * (temp_right - temp_left),
     xright=temp_right,ybottom=line2user(4.5,1),ytop=line2user(5.1,1),
     xpd=NA,col=NA,border="dodgerblue")
mtext("p \u2264 0.05",side=1,las=0,col="dodgerblue",cex=0.9,line=3.5,adj=1,at=temp_right)

mtext("A",font=2,cex=1.5,line=0.5,side=3,at=line2user(0.5,2))

rm(list=ls()[!ls() %in% KEEP])
```


```{r pDE_comparison_lig295, fig.width=7,fig.height=4}
D <- "lig295_DE"

load(paste0("~/Dropbox/GDB_archive/CMapCorr_files/",D,"_lig_FDR.RData"))
scoreDElig <- pDE_lig + 1e-4
scoreDElig[scoreDElig > 1] <- 1
scoreDElig <- -log10(scoreDElig)
load(paste0("~/Dropbox/GDB_archive/CMapCorr_files/",D,"_ligct_FDR.RData"))
scoreDEligct <- pDE_ligct + 1e-4
scoreDEligct[scoreDEligct > 1] <- 1
scoreDEligct <- -log10(scoreDEligct)
load(paste0("~/Dropbox/GDB_archive/CMapCorr_files/",D,"_rep_FDR.RData"))
scoreDErep <- pDE_rep + 1e-4
scoreDErep[scoreDErep > 1] <- 1
scoreDErep <- -log10(scoreDErep)

par(mar=c(3,7,0.5,2),mgp=2:0)
plot(NA,NA,xlim=c(0,4),ylim=c(1,15),
     xlab="Probability of at least # DE at FDR threshold occurring by chance",
     xaxt="n",yaxt="n",ylab=NA)
axis(1,at=-log10(c(1,0.1,0.05,0.01,0.001,0.0001)),
     labels=c("100%","10%","5%","1%","0.1%","0.01%"))
boxplot(scoreDElig,at=seq(1,by=4,length.out=4),
        add=T,horizontal=T,
        xaxt="n",yaxt="n",pch=".",cex=2,
        border=qualitative_hcl(3,palette="dark3")[1],
        col=qualitative_hcl(3,alpha=0.5,palette="dark3")[1])
boxplot(scoreDEligct,at=seq(2,by=4,length.out=4),
        add=T,horizontal=T,
        xaxt="n",yaxt="n",pch=".",cex=2,
        border=qualitative_hcl(3,palette="dark3")[2],
        col=qualitative_hcl(3,alpha=0.5,palette="dark3")[2])
temp_stat <- sapply(colnames(scoreDElig),function(X) 
  wilcox.test(scoreDElig[,X],scoreDEligct[,X])$p.value)
if (any(temp_stat <= 2.2e-16)) {
  mtext("*",side=2,line=3.3,font=2,
        at=seq(2,by=4,length.out=4)[temp_stat <= 2.2e-16])
}
boxplot(scoreDErep,at=seq(3,by=4,length.out=4),
        add=T,horizontal=T,
        xaxt="n",yaxt="n",pch=".",cex=2,
        border=qualitative_hcl(3,palette="dark3")[3],
        col=qualitative_hcl(3,alpha=0.5,palette="dark3")[3])
temp_stat <- sapply(colnames(scoreDElig),function(X) 
  wilcox.test(scoreDElig[,X],scoreDErep[,X])$p.value)
if (any(temp_stat <= 2.2e-16)) {
  mtext("*",side=2,line=3.3,font=2,
        at=seq(3,by=4,length.out=4)[temp_stat <= 2.2e-16])
}
# points(seq(1,by=5,length.out=4),colMeans(scoreDEct),
#        pch="-",cex=3,col=alpha("red",0.5))
# points(seq(2,by=5,length.out=4),colMeans(scoreDElig),
#        pch="-",cex=3,col=alpha("red",0.5))
# points(seq(3,by=5,length.out=4),colMeans(scoreDEligct),
#        pch="-",cex=3,col=alpha("red",0.5))
# points(seq(4,by=5,length.out=4),colMeans(scoreDErep),
#        pch="-",cex=3,col=alpha("red",0.5))
mtext(rep(c("Ligand","Lig / Line","Replicate",NA),4),
      col=c(qualitative_hcl(3,palette="dark3"),NA),
      side=2,line=0.1,las=2,at=1:16)
mtext(paste0(c(10,5,1,0.1),"% FDR"),
      side=4,line=0,at=seq(2,by=4,length.out=4),cex=0.9)

mtext("B",font=2,cex=1.5,line=-1,side=3,at=line2user(6,2))
```

The * denotes p < 2.2e-16 from Wilcoxon rank-sum test vs. "Ligand".  

# Supp Fig 2
Example plot to clarify what "averaging per ligand" means:  

```{r demo_meanZ, fig.width=6,fig.height=4}
LIG <- "TNF"
temp_gene <- rev(rownames(FDR_lig)[order(FDR_lig[,LIG])[seq(1,sum(FDR_lig[,LIG] < 1))]])
temp_Z <- sapply(temp_gene,function(X) 
  lvl4_data@mat[X,lvl4_data@cdesc$pert_iname == LIG],
  simplify=F)

par(mar=c(3,4,1,4),mgp=2:0)
plot(NA,NA,xlim=range(unlist(temp_Z)),ylim=c(0.5,length(temp_Z) + 0.5),
     yaxs="i",yaxt="n",ylab=NA,
     xlab=paste("Z-scores from all samples treated with",LIG))
abline(v=0,lty=2,col="red")
boxplot(temp_Z,horizontal=T,pch=".",cex=2,add=T,yaxt="n")
points(meanZ_lig[temp_gene,LIG],seq_along(temp_gene),
       pch=20,col=alpha("dodgerblue",0.8))
legend("topleft",bty="n",pch=20,col=alpha("dodgerblue",0.8),
       legend="Mean Z-score")
mtext(lvl4_data@rdesc[temp_gene,"pr_gene_symbol"],
      side=2,at=seq_along(temp_gene),line=0.1,las=2)
mtext(signif(FDR_lig[temp_gene,LIG],2),
      side=4,at=seq_along(temp_gene),line=0.1,las=2)
mtext("FDR",side=4,at=length(temp_gene) + 0.7,line=0.5,las=2,font=2)

rm(LIG)
rm(list=grep("^temp",ls(),value=T))
```

# Comparing prob to num
```{r fig.width=8,fig.height=6}
load("~/Dropbox/GDB_archive/CMapCorr_files/lig295_DE_lig_FDR.RData")
load("~/Dropbox/GDB_archive/CMapCorr_files/lig295_DE_ligct_FDR.RData")

temp_FDRcut <- 0.05

nDE_lig <- apply(FDR_lig,2,function(X) sum(X <= temp_FDRcut))
nDE_ligct <- apply(FDR_ligct,2,function(X) sum(X <= temp_FDRcut))

pDE_ligct <- pDE_ligct[,as.character(temp_FDRcut * 100)]
pDE_ligct <- pDE_ligct + 1e-4
pDE_ligct[pDE_ligct > 1] <- 1
pDE_ligct <- -log10(pDE_ligct)
pDE_lig <- pDE_lig[,as.character(temp_FDRcut * 100)]
pDE_lig <- pDE_lig + 1e-4
pDE_lig[pDE_lig > 1] <- 1
pDE_lig <- -log10(pDE_lig)

pDE_ligBYct <- tapply(pDE_ligct,as.vector(sub("^.+_","",names(pDE_ligct))),c)
pDE_ligBYct <- pDE_ligBYct[names(pDE_lig)]

layout(matrix(c(2,1,0,3),2),widths=c(3,1),heights=c(1,3))
par(mar=c(3,3,0,0),mgp=2:0)
plot(nDE_ligct,pDE_ligct,
     pch=".",cex=2,col=alpha("red",0.5),yaxt="n",
     xlim=range(c(nDE_lig,nDE_ligct)),ylim=range(c(pDE_lig,pDE_ligct)),
     xlab="Number of DE genes at FDR < 0.05 for each ligand",
     ylab="Probability of seeing at least that number of DE genes by chance")
axis(2,at=-log10(c(1,0.1,0.05,0.01,0.001,0.0001)),
     labels=c("100%","10%","5%","1%","0.1%","0.01%"))
# points(nDE_lig,pDE_lig,pch=".",cex=3,col=alpha("dodgerblue",0.5))
points(nDE_lig[sapply(pDE_ligBYct,max) > pDE_lig],
       pDE_lig[sapply(pDE_ligBYct,max) > pDE_lig],
       pch=".",cex=3,col=alpha("dodgerblue",0.5))
points(nDE_lig[sapply(pDE_ligBYct,max) <= pDE_lig],
       pDE_lig[sapply(pDE_ligBYct,max) <= pDE_lig],
       pch=20,col=alpha("dodgerblue",0.5))

par(mar=c(0,3,1,0))
boxplot(list(LIG=nDE_lig,LIGCT=nDE_ligct),horizontal=T,xaxt="n",
        col=alpha(c("dodgerblue","red"),0.5),border=c("dodgerblue","red"),pch=20)
par(mar=c(3,0,0,1))
boxplot(list(LIG=pDE_lig,LIGCT=pDE_ligct),horizontal=F,yaxt="n",
        col=alpha(c("dodgerblue","red"),0.5),border=c("dodgerblue","red"),pch=20)
```

# Ligands with consistent DE
```{r byCTdiff_num,fig.height=4,fig.width=8.5}
nDE_ligBYct <- tapply(nDE_ligct,as.vector(sub("^.+_","",names(nDE_ligct))),c)
nDE_ligBYct <- nDE_ligBYct[names(nDE_lig)]

byCTdiff <- sapply(names(nDE_lig),function(LIG) nDE_ligBYct[[LIG]] - nDE_lig[[LIG]],simplify=F)
byCTdiff <- byCTdiff[order(sapply(byCTdiff,median))]

par(mar=c(1.5,4,1,1),mgp=2:0)
boxplot(byCTdiff,pch=".",cex=2,xaxt="n")
abline(h=0)
mtext("Ligands",side=1,line=0.1)
mtext("Change in # DE at 5% FDR with cell line context",side=2,line=2)

```

# Supp Fig 3

```{r byCTdiff,fig.height=4,fig.width=8.5}
byCTdiff <- sapply(names(pDE_lig),function(LIG) pDE_ligBYct[[LIG]] - pDE_lig[[LIG]],simplify=F)
byCTdiff <- byCTdiff[order(sapply(byCTdiff,median))]

par(mar=c(1.5,4.5,1,1),mgp=2:0)
boxplot(byCTdiff,pch=".",cex=2,xaxt="n")
abline(h=0)
mtext("Ligands (ordered by median change in p-value)",side=1,line=0.1)
mtext(paste("Change in -log10 probability of # DE at 5% FDR ",
            "with cell line context vs. all ligand-treated samples",
            sep="\n"),
      side=2,line=2.5)
```


```{r include=FALSE}
plot(sapply(pDE_ligBYct,max),pDE_lig,pch=20)
abline(0,1)
table(sapply(pDE_ligBYct,max) > pDE_lig)

consistentLIG <- names(which(sapply(pDE_ligBYct,max) <= pDE_lig))

```

There were `r length(consistentLIG)` ligands out of 295 (`r round(100 * length(consistentLIG) / 295)`%) that **did not** show improvement in consistency of transcriptional response *in at least one cell line* when considering cell line context.  The following figures show the responsive genes in each case:  


```{r consistentLIG, fig.height=3,fig.width=8.5}
# LIG <- "ADM"
for (LIG in consistentLIG) {
  FDR <- FDR_ligct[,grep(paste0("_",LIG,"$"),colnames(FDR_ligct))]
  temp_cut <- 2
  temp_genes <- apply(FDR,1,function(X) sum(X <= 0.05) > temp_cut)
  while (sum(temp_genes) <= 1) {
    temp_cut <- temp_cut - 1
    if (temp_cut < 0) { break }
    temp_genes <- apply(FDR,1,function(X) sum(X <= 0.05) > temp_cut)
  }
  if (sum(temp_genes) <= 1) {
    print(paste("Ligand:",LIG))
    print(paste0("By ligand: ",nDE_lig[LIG]," DE at 5% FDR (p=",10^-pDE_lig[LIG],")"))
    print(paste0("By ligand per cell line: max ",
                 nDE_ligct[names(which.max(pDE_ligBYct[[LIG]]))],
                 " DE (p=",10^-pDE_ligct[names(which.max(pDE_ligBYct[[LIG]]))],")"))
    next
  }
  FDR <- FDR[temp_genes,]
  FDR <- -log10(FDR) * sign(meanZ_ligct[rownames(FDR),colnames(FDR)])
  temp_hROW <- hclust(dist(t(FDR)),method="ward.D2") 
  temp_hGENE <- hclust(dist(FDR),method="ward.D2")
  FDR <- FDR[temp_hGENE$order,temp_hROW$order]
  temp_max <- round(quantile(abs(FDR),.9))
  FDR[FDR > temp_max] <- temp_max; FDR[FDR < -temp_max] <- -temp_max
  
  
  par(mar=c(5,7,2,1),mgp=2:0)
  image(z=FDR / max(abs(FDR)),
        x=1:nrow(FDR),y=1:ncol(FDR),
        col=diverging_hcl(1000,palette="Blue-Red3"),
        breaks=seq(-1,1,length.out=1001),
        xaxt="n",yaxt="n",xlab=NA,ylab=NA)
  mtext(paste("Ligand:",LIG),line=0.4,font=2,cex=1.5,at=line2user(-6,2),adj=0)
  mtext("Genes",side=1,line=3.5,font=2,cex=1.2)
  mtext(lvl4_data@rdesc[rownames(FDR),"pr_gene_symbol"],
        side=1,las=2,at=1:nrow(FDR),line=0.1,cex=0.8)
  mtext("Cell lines",side=3,line=0.1,at=line2user(0.1,2),adj=1,font=2,cex=1.2)
  mtext(sub("^large_","",sub("_[A-Z0-9]+$","",colnames(FDR))),
        side=2,las=2,at=1:ncol(FDR),adj=1,line=0.1,cex=0.8)
  
  mtext(paste0("By ligand: ",nDE_lig[LIG]," DE at 5% FDR (p=",10^-pDE_lig[LIG],")"),
        side=3,line=1,at=line2user(0,4),adj=1,cex=0.8)
  mtext(paste0("By ligand per cell line: max ",
               nDE_ligct[names(which.max(pDE_ligBYct[[LIG]]))],
               " DE (p=",10^-pDE_ligct[names(which.max(pDE_ligBYct[[LIG]]))],")"),
        side=3,line=0.1,at=line2user(0,4),adj=1,cex=0.8)
  
  segments(x0=seq(line2user(5,2),line2user(1,2),length.out=1001),
           x1=seq(line2user(5,2),line2user(1,2),length.out=1001),
           y0=rep(line2user(1,1),1001),y1=rep(line2user(1.7,1),1001),
           xpd=NA,col=diverging_hcl(1001,palette="Blue-Red3"))
  text(c(line2user(5,2),line2user(3,2),line2user(1,2)),
       rep(line2user(2.2,1),3),
       labels=c(paste0("\u2264",-1 * round(max(abs(FDR)))),
                0,
                paste0("\u2265",round(max(abs(FDR))))),
       xpd=NA)
  text(line2user(3,2),line2user(3.5,1),xpd=NA,
       labels=paste("Signed","-log10 FDR",sep="\n"))
  
}
```