In RGLab/MASTdata: This is the data used for the MAST paper

Stimulated vs. Non-stimulated MAIT Cells

# library(devtools)
# install_git("http://github.com/RGLab/MASTDataPackage")
# install_git('http://github.com/RGLab/MAST')

# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <-
function(..., plotlist = NULL, file, cols = 1, layout = NULL) {
require(grid)

# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)

numPlots = length(plots)

# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots / cols)),
ncol = cols, nrow = ceiling(numPlots / cols))
}

if (numPlots == 1) {
print(plots[[1]])

} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

print(plots[[i]], vp = viewport(
layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col
))
}
}
}
suppressPackageStartupMessages({
library(MIMOSA)
library(plyr)
library(reshape)
library("ggplot2")
library("MASTDataPackage")
library("MAST")
library("GSEABase")
library("limma")
library(grid)
library("reshape2")
library(RColorBrewer)
library("data.table")
library("knitr")
library(psych)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
library(flowClust)
library(COMPASS)
library(scales)
library(grid)
library(stringr)
library(gridExtra)
library(class)
options(mc.cores = detectCores() - 1) #or set smaller if you find yourself running out of memory
})
LARGE_MEMORY_CORES <- 7 #for large memory operations
knitr::opts_chunk$set(list(
echo = FALSE,eval = TRUE,message = FALSE,error = FALSE,warning = FALSE,cache = TRUE,fig.width=10,fig.height=8))
plotheme <-
theme(
plot.background = element_blank(),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.line =
element_line(colour = "black")
)

opts_chunk$set(tidy=TRUE, cache=TRUE, messages=TRUE)

data(MASTDataPackage)
sca <- sca_mait
thres <- "adapt"
percent_expressed <- 0.2
FCTHRESHOLD <- log2(1.5)

sca <- subset(sca,nGeneOn > 4000)
eid <-
select(TxDb.Hsapiens.UCSC.hg19.knownGene,keys = fData(sca_mait)$entrez,keytype ="GENEID",columns = c("GENEID","TXNAME"))
ueid <- unique(na.omit(eid)$GENEID)
sca <- sca[,fData(sca)$entrez %in% ueid]

set.seed(123)
module <- "BTM"
min_gene_in_module <- 5

# thresholding
tt <- thresholdSCRNACountMatrix(2^exprs(sca) - 1, nbins = 20, min_per_bin = 30)

if (thres == "adapt") {
    exprs(sca) <- tt$counts_threshold
} else if (thres == "fixed") {
    mat <- exprs(sca)
    mat[mat < 1/log(2)] <- 0
    exprs(sca) <- mat
}

expressed_genes <- colMeans(exprs(sca) > 0) > percent_expressed
sca <- sca[, expressed_genes]
data <- exprs(sca)
cd <- cData(sca)
fd <- fData(sca)

## Data in long format for visualization
dt <- data.table(gene_id = fd$symbolid, t(data))
dt_long <- data.table(melt(dt))
dt_long[, `:=`(condition, ifelse(grepl("-Stim", variable), "Stim", "Unstim"))]

Hurdle model fitting

options(mc.cores=LARGE_MEMORY_CORES)
# ZLM (ridge regression for continuous, get standardized deviance residuals)
cond<-factor(cData(sca)$condition)
cond<-relevel(cond,"Unstim")
cData(sca)$condition<-cond
zlm <- zlm.SingleCellAssay(~condition + cngeneson, sca, method = "bayesglm", 
    ebayes = TRUE, ebayesControl = list(method = "MLE", model = "H1"))
zlm.nongeneson <- zlm.SingleCellAssay(~condition , sca, method = "bayesglm", 
    ebayes = TRUE, ebayesControl = list(method = "MLE", model = "H1"))
zlm.notrt <- zlm.SingleCellAssay(~cngeneson , sca, method = "bayesglm", 
    ebayes = TRUE, ebayesControl = list(method = "MLE", model = "H1"))
zlm.null <- zlm.SingleCellAssay(~1 , sca, method = "bayesglm", 
    ebayes = TRUE, ebayesControl = list(method = "MLE", model = "H1"))
# LRT
lrt <- lrTest(zlm, "condition")
lrt.nongeneson<-lrTest(zlm.nongeneson,"condition")

# Still need to merge on logFC
res_gene <- data.table(melt(lrt))
res_gene[,primerid:=gsub("\\.\\d+","",primerid)]
res_gene <- merge(res_gene, fd, by="primerid")
res_gene_hurdle <- res_gene[metric=="Pr(>Chisq)" & test.type=="hurdle"]
res_gene_hurdle[,adj:=p.adjust(value,"fdr")]
lfc<-getLogFC(zlm)[contrast=="conditionStim"]
setkey(lfc,primerid)
setkey(res_gene_hurdle,primerid)
res_gene_hurdle<-merge(lfc,res_gene_hurdle)

# merge on logFC
res_gene_nong <- data.table(melt(lrt.nongeneson))
res_gene_nong[,primerid:=gsub("\\.\\d+","",primerid)]
res_gene_nong <- merge(res_gene_nong, fd, by="primerid")
res_gene_hurdle_nong <- res_gene_nong[metric=="Pr(>Chisq)" & test.type=="hurdle"]
res_gene_hurdle_nong[,adj:=p.adjust(value,"fdr")]
lfc_nong<-getLogFC(zlm.nongeneson)[contrast=="conditionStim"]
setkey(lfc_nong,primerid)
setkey(res_gene_hurdle_nong,primerid)
res_gene_hurdle_nong<-merge(lfc_nong,res_gene_hurdle_nong)

nrow(res_gene_hurdle_nong[adj<0.01])
nrow(res_gene_hurdle[adj<0.01])

Comparison of ranks with and without CDR effect.

A=copy(res_gene_hurdle_nong)
B=copy(res_gene_hurdle)
A=na.omit(A[order(-abs(logFC),adj)])
B=na.omit(B[order(-abs(logFC),adj)])
C=merge(A[,.(.I,symbolid)],B[,.(.I,symbolid)],by="symbolid")
setkey(C,I.x)
ggplot(data.table(ngenes=seq(100,nrow(C),by=100),overlap=sapply(seq(100,nrow(C),by=100),function(x)mean(C[1:x,I.x<x&I.y<x]))))+geom_line(aes(x=ngenes,y=overlap))+theme_linedraw()+scale_x_continuous("Number of top genes Selected")+scale_y_continuous("Proportion overlap between CDR and non CDR models")+geom_vline(x = c(291,1491),lty=2)
#C[1:1491]
cat(paste(C[1:1491][!C[1:1491,I.x<1491&I.y<1491],symbolid],"\n"),file="../inst/extdata/withoutCDR.txt")
cat(paste(C[1:1491][C[1:1491,I.x<1491&I.y<1491],symbolid],"\n"),file="../inst/extdata/withCDR.txt")
cat(paste(C[,symbolid],"\n"),file="../inst/extdata/universe.txt")

#gsea for limma exclusive genes
gsea_exclusive_other_methods<-readRDS("../inst/extdata/significant_genes_q001_fc13.rds") #From the anlaysis by Masanao Yajima
setDT(gsea_exclusive_other_methods)

#limma only
cat(paste0(gsea_exclusive_other_methods[limma_ng&!zlm_ng,symbolid],"\n"),file="../inst/extdata/Limma_only.txt")
cat(paste0(gsea_exclusive_other_methods[edger_ng&!zlm_ng,symbolid],"\n"),file="../inst/extdata/edger_only.txt")
cat(paste0(gsea_exclusive_other_methods[deseq_ng & !zlm_ng,symbolid],"\n"),file="../inst/extdata/deseq_only.txt")

Analysis of deviance for ngeneson (Supplementary Figure 2)

library(scales)

panelb=rbind(data.table(melt((zlm.notrt@deviance-zlm@deviance)/zlm.null@deviance),model="trt"),data.table(melt((zlm.nongeneson@deviance-zlm@deviance)/zlm.null@deviance),model="cdr"))
setnames(panelb,c("gene","component","value","model"))
panelb=dcast(panelb,...~model)

#Shalek data for this SFig
filtered<-sca_alex
percent_expressed<-0.1
thres<-"adapt"
tt <- thresholdSCRNACountMatrix(2^exprs(filtered) - 1, nbins = 20, min_per_bin = 30)
if (thres == "adapt") {
    exprs(filtered) <- tt$counts_threshold
} else if (thres == "fixed") {
    mat <- exprs(filtered)
    mat[mat < 1/log(2)] <- 0
    exprs(filtered) <- mat
}
expressed_genes <- colMeans(exprs(filtered) > 0) > percent_expressed
filtered <- filtered[, expressed_genes]
cData(filtered)$Stim <- factor(cData(filtered)$Stim,levels=c("Unstimulated","LPS","PAM","PIC"))
cData(filtered)$Time <- factor(cData(filtered)$Time)
filtered_nobaseline <- subset(filtered,!Stim%in%c("Unstimulated")) #drop Time 0
cData(filtered_nobaseline)$Stim <- factor(cData(filtered_nobaseline)$Stim)
cData(filtered_nobaseline)$Time <- factor(cData(filtered_nobaseline)$Time)
ng <- zlm.SingleCellAssay(~cngeneson+Stim/Time,sca=filtered_nobaseline,method="ridge",ebayes=TRUE)
nong <- zlm.SingleCellAssay(~Stim/Time,sca=filtered_nobaseline,method="ridge",ebayes=TRUE)
notrt <- zlm.SingleCellAssay(~cngeneson,sca=filtered_nobaseline,method="ridge",ebayes=TRUE)
nullmodel = zlm.SingleCellAssay(~1,sca=filtered_nobaseline,method="ridge",ebayes=TRUE)


panelc=rbind(data.table(melt((notrt@deviance-ng@deviance)/nullmodel@deviance),model="trt"),data.table(melt((nong@deviance-ng@deviance)/nullmodel@deviance),model="cdr"))
setnames(panelc,c("gene","component","value","model"))
panelc=dcast(panelc,...~model)
#ggplot(panelc)+aes(x=trt,y=cdr,col=component)+geom_point()

toplot = rbind(cbind(panelc,dataset="mdc"),cbind(panelb,dataset="MAIT"))

p3 = ggplot(toplot)+geom_hex()+aes(x=cdr,y=trt)+facet_grid(component~dataset)+theme_linedraw()+scale_color_discrete("Component",labels=c("Continouos","Discrete"))+scale_x_continuous("Proportion of deviance explained by CDR")+scale_y_continuous("Proportion of deviance explained by treatment")+scale_fill_distiller(type = "div",palette = 2,space = "Lab",trans="log",breaks=c(1,10,50,100,250,500))
pdf(file = "../inst/extdata/output/Supplementary_Figure_2.pdf",width = 8,height =
      5)
      print(p3)
      dev.off()

cat("Mean % Deviance")
kable(ddply(data.table(toplot)[cdr>0],.(dataset,component),function(x)mean(x$cdr)*100))
cat("% Deviance 90th Percentile")
kable(ddply(data.table(toplot)[cdr>0],.(dataset,component),function(x)quantile(x$cdr,0.9)*100))
print(p3)

Visualization of most differentially expressed genes

# Let's look at the top 40
genes_to_plot <- res_gene_hurdle[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),adj)][1:40,symbolid]
ggplot(dt_long[gene_id%in%genes_to_plot][,gene_id:=factor(gene_id,levels=genes_to_plot)],aes(x=condition, y=value,color=condition))+geom_violin()+geom_jitter()+facet_wrap(~gene_id)+ggtitle("Top 40 DE Genes in Activated MAIT Cells")

Heatmap of MAITs based on most differentially expressed genes

Top 100 DE genes between stimulated and non-stimulated MAIT cells.

fd<-fData(sca)
fd$primerid<-fd$symbolid
colnames(sca)<-fd$symbolid
sca@featureData@data<-fd
FCTHRESHOLD<-0.58
genes_to_plot <- res_gene_hurdle[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),adj)][1:100,symbolid]
COMPASS::pheatmap(exprs(sca[,genes_to_plot]),row_annotation=cData(sca)[,"condition",drop=FALSE],main="Top 100 DE genes",col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)))

Some of the activated MAITs have transcriptional profiles more similar to unactivated MAITS.

GSEA

#bootstrap
cache(boots <- bootVcov1(zlm, 50))

module_file <- list.files("../inst/extdata/", pattern = module, full.names = TRUE)
gene_set <- getGmt(module_file)
gene_ids <- geneIds(gene_set)
gene_ids=gene_ids[!names(gene_ids)%like%"TBA"&!names(gene_ids)%like%"B cell"]
saveRDS(gene_ids,file="../inst/extdata/MAIT_Modules_symbols.rds")
sets_indices <- ids2indices(gene_ids, fd$symbolid)
# Only keep modules with at least min_gene_in_module
sets_indices <- sets_indices[sapply(sets_indices, length) >= min_gene_in_module]
cache(gsea <- gseaAfterBoot(zlm, boots, sets_indices, CoefficientHypothesis("conditionStim")))
t_stat <- melt(calcZ(gsea, testType = "normal"))
t_stat_wide <- dcast(t_stat, set ~ comp + metric)
t_stat <- data.frame(calcZ(gsea, testType = "normal", combined = "stouffer"))
t_stat$set <- rownames(t_stat)
t_stat_comb <- data.table(merge(t_stat_wide, t_stat, by = "set"))
setorder(t_stat_comb,P)
t_stat_comb[,adj:=p.adjust(P,"fdr")]

##Filter by effect size: discrete odds ratio of 1
t_stat_comb=t_stat_comb[set%in%names(which(abs(gsea[,"disc","stat","test"]-gsea[,"disc","stat","null"])>log(2.5) ))|set%in%names(which(abs(gsea[,"cont","stat","test"]-gsea[,"cont","stat","null"])>log2(1.5)))]

for (i in c(1:10,12,13)) {
    genes_to_plot <- unlist(gene_ids[names(gene_ids) == as.character(t_stat_comb$set[i])])
    gp <- ggplot(dt_long[gene_id %in% genes_to_plot], aes(x = condition, y = value, 
        color = condition)) + geom_violin() + geom_jitter() + facet_wrap(~gene_id,scale="free") + 
        labs(title = paste0(t_stat_comb$set[i], ", Z=", round(t_stat_comb$Z[i], 
            2), ", P=", round(t_stat_comb$P[i], 2)))
    print(gp)
}

options("mc.cores" = LARGE_MEMORY_CORES)

score_4_hook <- function(x) {
  if (all(x@fitted)) {
    class(x@fitC) <- c("glm","lm")
    class(x@fitD) <- c("bayesglm","glm","lm")
    wh <- !colnames(x@modelMatrix) %like% "cngeneson"
    wh2 <- !names(coef(x@fitC)) %like% "cngeneson"
    fc <-
      x@modelMatrix[,!wh,drop = FALSE] %*% coef(x@fitC)[!wh2,drop = FALSE] #continuous CDR effect
    fd <-
      arm::invlogit(x@modelMatrix[,c("(Intercept)","cngeneson"),drop = FALSE] %*%
                      coef(x@fitD)[c("(Intercept)","cngeneson"),drop = FALSE]) #discrete CDR effect
    R <-
      matrix((x@response - fc * fd),nrow = 1) #residuals corrected for ngeneson in the continuous part
    colnames(R) <- names(residuals(x@fitD))
    R[x@response == 0] <- 0 - fd
    R
  }
}

zlm.scores <-
  zlm.SingleCellAssay(
    ~condition + cngeneson, sca, method = "bayesglm",
    ebayes = TRUE, ebayesControl = list(method = "MLE", model = "H1"),hook =
      score_4_hook
  )

scores <- do.call(rbind,zlm.scores@hookOut)
rownames(scores) <- names(zlm.scores@hookOut)

genes_to_plot <- sapply(1:nrow(t_stat_comb),function(i) {
unlist(gene_ids[names(gene_ids) == as.character(t_stat_comb$set[i])])
})
names(genes_to_plot) <- t_stat_comb$set

##Include the DE MAIT signature
mait_signature_genes_up = res_gene_hurdle[abs(logFC) > FCTHRESHOLD &
adj < 0.01][order(-abs(logFC),adj)][1:50,][logFC > 0,symbolid]
mait_signature_genes_down = res_gene_hurdle[abs(logFC) > FCTHRESHOLD &
adj < 0.01][order(-abs(logFC),adj)][1:50,][logFC < 0,symbolid]

genes_to_plot$"MAIT Activation [Up]" <- mait_signature_genes_up
genes_to_plot$"MAIT Activation [Down]" <- mait_signature_genes_down

genes_to_plot <- ldply(lapply(genes_to_plot,function(x)
data.frame(x)))
colnames(genes_to_plot) <- c("set","gene")
scores <- (melt(scores))

colnames(scores) <- c("gene","cell","score")

scores <- (merge(genes_to_plot,scores,by = "gene",all.y = TRUE))
condition <- cData(sca)[,"condition",drop = FALSE]
condition$cell <- rownames(condition)
scores <- (merge(scores,condition,by = "cell",all.x = TRUE))
scores <- na.omit(scores)
scores <- data.table(scores)
scores = scores[,set := factor(set)]
#t_stat_comb = t_stat_comb[,adj := p.adjust(P)]
scores$set <-
sapply(as.character(scores$set),function(x)
Kmisc::wrap(x,width = 25))
sets <-
sapply(t_stat_comb[order(-abs(t_stat_comb$Z),t_stat_comb$adj),set][c(1:9)],function(x)
Kmisc::wrap(x,25))
#sets <- gsub("\n$","",gsub(".\\(M.+$","",sets[!sets %like% "TBA"]))
gsea_scores <-
ggplot(scores[set %in% sets,list(score =
mean(score)),list(cell,set,condition)][,set := factor(set,levels = sets)]) +  geom_violin(show_guide =
FALSE) + geom_jitter(show_guide = FALSE) + aes(x = condition,y = score,col =
condition) + facet_wrap( ~set) + theme_linedraw() + theme(
plot.background = element_blank(),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.line =
element_line(colour = "black"),legend.position = "bottom"
) + scale_color_manual("",values = c(brewer.pal(name = "PiYG",n = 3)[c(1)],"#5729EE"))

Heatmap and PCA of module scores

#heatmap of module scores
mat = dcast(unique(scores[,list(cell,set,score)]),cell~set,value.var="score",fun.aggregate=mean)
#colnames(mat) = gsub("\\n"," ",colnames(mat))
rownames(mat)=mat[,1]
mat=mat[,-1L]
COMPASS::pheatmap(mat,clustering_distance_rows="manhattan",clustering_distance_cols="manhattan")

library(stringr)
plot(prcomp(mat)$x[,1:2],col=factor(str_split_fixed(rownames(mat),"-",4)[,3]),pch=20)
legend("bottom",c("Stim","Non-stim"),col=c("black","red"),pch=20)

plot(prcomp(mat))

#sort(prcomp(mat)$rotation[,1])

setorder(t_stat_comb,Z,adj)
ylevels=as.character(t_stat_comb[adj<0.01&!set %like% "TBA",set])
gsea_heatmap <-
ggplot(data.table:::melt.data.table(t_stat_comb[order(Z,adj)][adj < 0.01 &
!set %like% "TBA"],id = "set")[variable %like% "Z"]) + aes(y = set,x = variable,fill =
value) + geom_tile() + scale_y_discrete(limits = ylevels) +
  scale_fill_gradient2(
"GSEA Z-score",low = ("#E9A3C9"),mid = ("#F7F7F7"),high = ("#A1D76A"),space =
"Lab"
) + theme(legend.position = "bottom",axis.text.x = element_text(angle =
45,hjust = 1))

pdf("../inst/extdata/output/Figure3.pdf",width=15,height=7.5)
layout<-cbind(cbind(matrix(2,ncol=1,nrow=10),matrix(1,nrow=10,ncol=9)),cbind(matrix(3,nrow=10,ncol=9),matrix(2,ncol=1,nrow=10)))
lay<-grid.layout(nrow = nrow(layout),ncol=ncol(layout))
grid.newpage()
vp<-viewport(layout=lay)
pushViewport(vp)
ind<-which(layout==1,T)
vp<-viewport(layout.pos.row =ind[,1],layout.pos.col=ind[,2])
print(gsea_scores,vp=vp)
grid.text("A",x = 0.05,y=0.95)
ind<-which(layout==3,T)
vp<-viewport(layout.pos.row =ind[,1],layout.pos.col=ind[,2])
print(gsea_heatmap,vp=vp)
grid.text("B",x = 0.5,y=0.95)
dev.off()

png("../inst/extdata/output/Figure3.png",width=15*300,height=7.5*300,res=300)
layout<-cbind(cbind(matrix(2,ncol=1,nrow=10),matrix(1,nrow=10,ncol=9)),cbind(matrix(3,nrow=10,ncol=9),matrix(2,ncol=1,nrow=10)))
lay<-grid.layout(nrow = nrow(layout),ncol=ncol(layout))
grid.newpage()
vp<-viewport(layout=lay)
pushViewport(vp)
ind<-which(layout==1,T)
vp<-viewport(layout.pos.row =ind[,1],layout.pos.col=ind[,2])
print(gsea_scores,vp=vp)
grid.text("A",x = 0.05,y=0.95)
ind<-which(layout==3,T)
vp<-viewport(layout.pos.row =ind[,1],layout.pos.col=ind[,2])
print(gsea_heatmap,vp=vp)
grid.text("B",x = 0.5,y=0.95)
dev.off()

Figure 3 GSEA Results for MAIT Data Set

layout<-cbind(cbind(matrix(2,ncol=1,nrow=10),matrix(1,nrow=10,ncol=9)),cbind(matrix(3,nrow=10,ncol=9),matrix(2,ncol=1,nrow=10)))
lay<-grid.layout(nrow = nrow(layout),ncol=ncol(layout))
grid.newpage()
vp<-viewport(layout=lay)
pushViewport(vp)
ind<-which(layout==1,T)
vp<-viewport(layout.pos.row =ind[,1],layout.pos.col=ind[,2])
print(gsea_scores,vp=vp)
grid.text("A",x = 0.05,y=0.95)
ind<-which(layout==3,T)
vp<-viewport(layout.pos.row =ind[,1],layout.pos.col=ind[,2])
print(gsea_heatmap,vp=vp)
grid.text("B",x = 0.5,y=0.95)

genes_to_plot <- res_gene_hurdle[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),adj)][1:100,symbolid]

E<-exprs(sca[,genes_to_plot])
R<-cData(sca)[,"condition",drop=FALSE]
rownames(E)<-gsub("_.+","",gsub(".+Unstim-","U",gsub(".+Stim-","S",rownames(E))))
rownames(R)<-gsub("_.+","",gsub(".+Unstim-","U",gsub(".+Stim-","S",rownames(R))))

COMPASS::pheatmap(E,row_annotation=R,col=(colorRampPalette(colors = brewer.pal(name="Blues",n=10))(10)),cluster_rows=TRUE,clustering_method="ward",scale="none",row_annotation_colors = list(condition=c(Unstim=brewer.pal(name="PiYG",n=3)[1],Stim="#5729EE")),row_annotation_legend = FALSE,fontsize = 4)

gtr<-grid.grab(wrap=TRUE)
fig_sets<-"MAIT"

mait_scores<-ggplot(scores[set%like%fig_sets,list(score=mean(score)),list(cell,set,condition)])+geom_violin()+geom_jitter()+aes(x=condition,y=score,col=condition)+facet_wrap(~set,nrow = 2)+#geom_text(size=2.5,aes(label=gsub("Unstim-","U",gsub("Stim-","S",gsub("_.*","",gsub("1-MAIT-","",cell))))),data=scores[(cell%like%c("Stim-C15")|cell%like%c("Stim-C46")|cell%like%"Stim-C54"|cell%like%"Stim-C58"|cell%like%"Stim-C94"|cell%like%"Stim-C61")&set%like%fig_sets,list(score=mean(score)),list(cell,set,condition)],show_guide=FALSE)+
  theme_linedraw(base_size = 9)+theme(plot.background = element_blank(),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.line=element_line(colour="black"),legend.position="bottom")+scale_color_manual(values = c(brewer.pal(name = "PiYG",n=3)[c(1)],"#5729EE"))

lout<-matrix(c(1,1,1,1,1,1,1,1,2,2,2,2,
               1,1,1,1,1,1,1,1,3,3,3,3,
               1,1,1,1,1,1,1,1,3,3,3,3,
               1,1,1,1,1,1,1,1,3,3,3,3,
               1,1,1,1,1,1,1,1,3,3,3,3,
               1,1,1,1,1,1,1,1,3,3,3,3,
               1,1,1,1,1,1,1,1,3,3,3,3,
               1,1,1,1,1,1,1,1,3,3,3,3,
               c(rep(5,8),3,3,3,3)),ncol=12,byrow=TRUE)
pdf(file="../inst/extdata/output/Figure2.pdf",width=8*1.25,height=4*1.25)
grid.newpage()
lay<-grid.layout(nrow(lout), ncol(lout))
pushViewport(viewport(layout=lay,name = "top"))
matchind<-which(lout==c(1,1),T)
vp<-viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2],name="hmap")
pushViewport(vp)
grid.draw(gtr)
seekViewport("top")
matchind<-which(lout==c(3),T)
vp=viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2])
print(mait_scores,vp=vp)
grid.draw(textGrob("A",0.05,0.95))
grid.draw(textGrob("B",0.65,0.95))
dev.off()

png(file="../inst/extdata/output/Figure2.png",width=300*8*1.25,height=300*4*1.25,res = 300)
grid.newpage()
lay<-grid.layout(nrow(lout), ncol(lout))
pushViewport(viewport(layout=lay,name = "top"))
matchind<-which(lout==c(1,1),T)
vp<-viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2],name="hmap")
pushViewport(vp)
grid.draw(gtr)
grid.ls(viewport=TRUE,grob=FALSE)
seekViewport("top")
matchind<-which(lout==c(3),T)
vp=viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2])
print(mait_scores,vp=vp)
grid.draw(textGrob("A",0.05,0.95))
grid.draw(textGrob("B",0.65,0.95))
dev.off()

Figure 2 - Heatmap of Differentially Expressed Genes and Activated MAIT Module Scores

grid.newpage()
lay<-grid.layout(nrow(lout), ncol(lout))
pushViewport(viewport(layout=lay,name = "top"))
matchind<-which(lout==c(1,1),T)
vp<-viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2],name="hmap")
pushViewport(vp)
grid.draw(gtr)
seekViewport("top")
matchind<-which(lout==c(3),T)
vp=viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2])
print(mait_scores,vp=vp)
grid.draw(textGrob("A",0.05,0.95))
grid.draw(textGrob("B",0.65,0.95))

zlm.resid <- zlm.SingleCellAssay(~condition + cngeneson, sca, method = "bayesglm", 
    ebayes = TRUE, ebayesControl = list(method = "MLE", model = "H1"),hook=deviance_residuals_hook)

residuals<-do.call(rbind,zlm.resid@hookOut)
rownames(residuals)<-fData(sca)[,"symbolid"]
colnames(residuals)<-rownames(cData(sca))

genes_to_plot <- res_gene_hurdle[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),value)][1:20,symbolid]

#correlation amongst cells
C.s<-cor(residuals[genes_to_plot,cData(sca)[,"condition"]%in%"Stim"])
C.u<-cor(residuals[genes_to_plot,cData(sca)[,"condition"]%in%"Unstim"])
C<-cor(residuals[genes_to_plot,])

COMPASS::pheatmap(C,row_annotation=cData(sca)[,"condition",drop=FALSE],show_rownames=TRUE,show_colnames=FALSE,main="Correlation of residuals across all cells using top 20 DE genes",col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21))

options("mc.cores"=LARGE_MEMORY_CORES)
zlm.resid.nong <- zlm.SingleCellAssay(~condition , sca, method = "bayesglm", 
    ebayes = TRUE, ebayesControl = list(method = "MLE", model = "H1"),hook=deviance_residuals_hook)
lrt.nong <- lrTest(zlm.resid.nong, "condition")

# Still need to merge on logFC
res_gene_nong <- data.table(melt(lrt.nong))
res_gene_nong <- merge(res_gene_nong, fd, by="primerid")
res_gene_hurdle_nong <- res_gene_nong[metric=="Pr(>Chisq)" & test.type=="hurdle"]
res_gene_hurdle_nong[,adj:=p.adjust(`value`,"fdr")]
lfc_nong<-getLogFC(zlm.resid.nong)[contrast=="conditionStim"]
setkey(lfc_nong,primerid)
setkey(res_gene_hurdle_nong,primerid)
res_gene_hurdle_nong<-merge(lfc_nong,res_gene_hurdle_nong)


residuals.nong<-do.call(rbind,zlm.resid.nong@hookOut)
rownames(residuals.nong)<-fData(sca)[,"symbolid"]
colnames(residuals.nong)<-rownames(cData(sca))

genes_to_plot <- res_gene_hurdle_nong[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),value)][1:20,symbolid]
nrow(res_gene_hurdle_nong[abs(logFC)>FCTHRESHOLD&adj<0.01])
#correlation amongst cells
C<-cor(residuals.nong[genes_to_plot,])

COMPASS::pheatmap(C,row_annotation=cData(sca)[,"condition",drop=FALSE],show_rownames=FALSE,show_colnames=FALSE,clustering_method="ward",main="Correlation of residuals across all cells using top 20 DE genes (without ngeneson)",col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21))

residuals<-do.call(rbind,zlm.resid@hookOut)
rownames(residuals)<-fData(sca)[,"symbolid"]
colnames(residuals)<-rownames(cData(sca))
genes_to_plot <- res_gene_hurdle[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),value)][1:50,symbolid]
CC.s<-cor(t(residuals)[cData(sca)[,"condition"]%in%"Stim",genes_to_plot])
CC.u<-cor(t(residuals)[cData(sca)[,"condition"]%in%"Unstim",genes_to_plot])
CC.s<-data.table(melt(CC.s))
CC.u<-data.table(melt(CC.u))
setnames(CC.s,c("row","col","cor.s"))
setnames(CC.u,c("row","col","cor.u"))
CC.s<-CC.s[!is.na(cor)]
CC.u<-CC.u[!is.na(cor)]
setkey(CC.s,row,col)
setkey(CC.u,row,col)
CC.su<-merge(CC.s,CC.u)
CC.su[,cor.d:=cor.s-cor.u]
setorder(CC.su,cor.d)

toplot<-dcast(CC.su,row~col,value.var ="cor.s")[,-c(1)]
rownames(toplot)<-colnames(toplot)
o<-hclust(as.dist(1-cor(toplot)),method="ward")$order


toplot.ns<-dcast(CC.su,row~col,value.var ="cor.u")[,-c(1)]
rownames(toplot.ns)<-colnames(toplot.ns)
COMPASS::pheatmap(toplot.ns[o,o],col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21),cluster_rows=FALSE,cluster_cols=FALSE,legend=FALSE,fontsize = 7)
nostim<-grid.grab()
COMPASS::pheatmap(toplot[o,o],col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21),cluster_rows=FALSE,cluster_cols=FALSE,fontsize = 7)
stim<-grid.grab()
cutree(hclust(as.dist(1-cor(toplot)),method="ward"),3)
genes<-c("IFNG","GZMB","JUND","CD69","FOS","GZMA","TXN")

# genes_to_plot <- res_gene_hurdle_nong[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),value)][1:50,symbolid]
# CC.s<-cor(t(residuals.nong)[cData(sca)[,"condition"]%in%"Stim",genes_to_plot])
# CC.u<-cor(t(residuals.nong)[cData(sca)[,"condition"]%in%"Unstim",genes_to_plot])
# CC.s<-data.table(melt(CC.s))
# CC.u<-data.table(melt(CC.u))
# setnames(CC.s,c("row","col","cor.s"))
# setnames(CC.u,c("row","col","cor.u"))
# CC.s<-CC.s[!is.na(cor)]
# CC.u<-CC.u[!is.na(cor)]
# setkey(CC.s,row,col)
# setkey(CC.u,row,col)
# CC.su<-merge(CC.s,CC.u)
# CC.su[,cor.d:=cor.s-cor.u]
# setorder(CC.su,cor.d)
# 
# 
# toplot<-dcast(CC.su,row~col,value.var ="cor.u")[,-c(1)]
# rownames(toplot)<-colnames(toplot)
# COMPASS::pheatmap(toplot,main="Correlation among top DE genes in unstimulated cells (nongneson)",col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21))
# 
# toplot<-dcast(CC.su,row~col,value.var ="cor.s")[,-c(1)]
# rownames(toplot)<-colnames(toplot)
# COMPASS::pheatmap(toplot,main="Correlation among top DE genes in stimulated cells (no ngeneson)",col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21))

library(grid)
PCbiplot <- function(PC, x="PC1", y="PC2", colors=c('black', 'black', 'red', 'red'),arrow_colors="red",point_colors="black",arrow_alpha=0.75,gene_subset=NULL) {
    # PC being a prcomp object
    data <- data.frame(obsnames=row.names(PC$x), PC$x)
    plot <- ggplot(data, aes_string(x=x, y=y)) +geom_point(color=point_colors)#+ geom_text(alpha=.8, size=3, aes(label=obsnames), color=point_colors)
    plot <- plot + geom_hline(aes(0), size=.2) + geom_vline(aes(0), size=.2, color=colors[2])
    datapc <- data.frame(varnames=rownames(PC$rotation), PC$rotation)
    mult <- min(
        (max(data[,y]) - min(data[,y])/(max(datapc[,y])-min(datapc[,y]))),
        (max(data[,x]) - min(data[,x])/(max(datapc[,x])-min(datapc[,x])))
        )
    datapc <- transform(datapc,
            v1 = .7 * mult * (get(x)),
            v2 = .7 * mult * (get(y))
            )
    plot <- plot + coord_equal() + geom_text(data=datapc[gene_subset,], aes(x=v1, y=v2, label=varnames), size = 4, vjust=1, color=arrow_colors,alpha=arrow_alpha)
    plot <- plot + geom_segment(data=datapc[gene_subset,], aes(x=0, y=0, xend=v1, yend=v2), arrow=arrow(length=unit(0.2,"cm")), alpha=arrow_alpha, color=arrow_colors)
    plot
}
residuals<-do.call(rbind,zlm.resid@hookOut)
rownames(residuals)<-fData(sca)[,"symbolid"]
colnames(residuals)<-rownames(cData(sca))
point_colors<-factor(cData(sca)[,"condition"],labels=c(Unstim=brewer.pal(name="PiYG",n=1)[1],Stim="#5729EE"))

genes_to_plot <- res_gene_hurdle[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),value)][1:50,symbolid]
pcaplot<-PCbiplot(prcomp(t(residuals[genes_to_plot,])),point_colors=point_colors,arrow_colors="blue",arrow_alpha=0.75,gene_subset=genes)+theme_bw()

lout<-matrix(rep(c(1,2,3),each=100),ncol=30)
lout<-rbind(rep(4,30),lout,rep(4,30))
pdf(file="../inst/extdata/output/Figure4.pdf",width=15,height=5)
grid.newpage()
lay<-grid.layout(nrow(lout), ncol(lout))
pushViewport(viewport(layout=lay,name = "top"))
matchind<-which(lout==c(1),T)
vp<-viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2],name="hmap")
pushViewport(vp)
grid.draw(nostim)
seekViewport("top")
matchind<-which(lout==c(2),T)
vp=viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2])
pushViewport(vp)
grid.draw(stim)
seekViewport("top")
matchind<-which(lout==c(3),T)
vp=viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2])
print(pcaplot,vp=vp)
grid.draw(textGrob("A",0.01,0.95))
grid.draw(textGrob("B",0.3,0.95))
grid.draw(textGrob("C",0.65,0.95))
dev.off()

lout<-matrix(rep(c(1,2,3),each=100),ncol=30)
lout<-rbind(rep(4,30),lout,rep(4,30))
png(file="../inst/extdata/output/Figure4.png",width=300*15,height=300*5,res=300)
grid.newpage()
lay<-grid.layout(nrow(lout), ncol(lout))
pushViewport(viewport(layout=lay,name = "top"))
matchind<-which(lout==c(1),T)
vp<-viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2],name="hmap")
pushViewport(vp)
grid.draw(nostim)
seekViewport("top")
matchind<-which(lout==c(2),T)
vp=viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2])
pushViewport(vp)
grid.draw(stim)
seekViewport("top")
matchind<-which(lout==c(3),T)
vp=viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2])
print(pcaplot,vp=vp)
grid.draw(textGrob("A",0.01,0.95))
grid.draw(textGrob("B",0.3,0.95))
grid.draw(textGrob("C",0.65,0.95))
dev.off()

Figure 4 - Residual Correlations and PCA

grid.newpage()
lay<-grid.layout(nrow(lout), ncol(lout))
pushViewport(viewport(layout=lay,name = "top"))
matchind<-which(lout==c(1),T)
vp<-viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2],name="hmap")
pushViewport(vp)
grid.draw(nostim)
seekViewport("top")
matchind<-which(lout==c(2),T)
vp=viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2])
pushViewport(vp)
grid.draw(stim)
seekViewport("top")
matchind<-which(lout==c(3),T)
vp=viewport(layout.pos.row = matchind[,1],layout.pos.col = matchind[,2])
print(pcaplot,vp=vp)
grid.draw(textGrob("A",0.01,0.95))
grid.draw(textGrob("B",0.3,0.95))
grid.draw(textGrob("C",0.65,0.95))

# library(Matrix)
# library(igraph)
# library(glasso)
# residuals<-do.call(rbind,zlm.resid@hookOut)
# rownames(residuals)<-fData(sca)[,"symbolid"]
# colnames(residuals)<-rownames(cData(sca))
# 
# getNetwork<-function(residuals=NULL,genes=NULL,rho=0.35,main=NULL,plot=TRUE){
#   g<-glasso(cor(t(residuals[genes,]),method="kendall"),rho=rho)
#   colnames(g$w)<-genes
#   rownames(g$w)<-genes
#   g<-graph.adjacency(g$w,mode="undirected",diag=FALSE,weighted=TRUE)
#   sort(igraph:::degree(g))
#   G<-subgraph(g,which(igraph:::degree(g)>0))
#   set.seed(10)
#   lay<-layout.fruchterman.reingold(G,area=vcount(G)^2.3,repulserad=vcount(G)^2)
#   if(plot){
#     plot(G,layout=lay,edge.width=3,vertex.size=18,vertex.label.cex=0.5,main=main)
#   }
#   G
# }
# 
# genes_to_plot <- res_gene_hurdle[,adj:=p.adjust(value,"fdr")][abs(logFC)>3&adj<0.01][order(-abs(logFC),adj)][,symbolid]
# G<-getNetwork(residuals[,cData(sca)[,"condition"]%in%c("Stim","Unstim")],genes_to_plot,0.35,main="Residuals network in MAIT")
# 
# 
# COMPASS::pheatmap(exprs(sca[,vertex.attributes(G)$name]),row_annotation=cData(sca)[,"condition",FALSE],main="Genes from residuals network in MAIT (with ngeneson)",col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21))
# COMPASS::pheatmap(cor(t(residuals[vertex.attributes(G)$name,])),main="gene-gene correlation from network estimation (with ngeneson)",col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21))
# gene_colors<-c("red","blue")[as.numeric(factor(sign(res_gene_hurdle[symbolid%in%vertex.attributes(G)$name,][order(symbolid),logFC])))]
# PCbiplot(prcomp((residuals[sort(vertex.attributes(G)$name),])),arrow_color=NA,point_colors=gene_colors)+ggtitle("PCA of residuals of genes from network estimate with ngeneson")

# residuals.nong<-do.call(rbind,zlm.resid.nong@hookOut)
# rownames(residuals.nong)<-fData(sca)[,"symbolid"]
# colnames(residuals.nong)<-rownames(cData(sca))
# genes_to_plot <- res_gene_hurdle_nong[,adj:=p.adjust(value,"fdr")][abs(logFC)>3&adj<0.01][order(-abs(logFC),adj)][,symbolid]
# G<-getNetwork(residuals.nong[,cData(sca)[,"condition"]%in%c("Stim","Unstim")],genes_to_plot,0.35,main="Residuals network in MAIT (no ngeneson correction)")
# 
# 
# gene_colors<-c("red","blue")[as.numeric(factor(sign(res_gene_hurdle_nong[symbolid%in%vertex.attributes(G)$name,][order(symbolid),logFC])))]
# 
# COMPASS::pheatmap(exprs(sca[,vertex.attributes(G)$name]),row_annotation=cData(sca)[,"condition",FALSE],main="Genes from residuals network in MAIT (without ngeneson)",col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21))
# COMPASS::pheatmap(cor(t(residuals.nong[vertex.attributes(G)$name,])),main="Gene-gene correlation from network estimation (without ngeneson)",col=rev(colorRampPalette(colors = brewer.pal(name="PiYG",n=10))(20)),breaks = seq(-1,1,l=21))
# PCbiplot(prcomp((residuals.nong[vertex.attributes(G)$name,])),arrow_color=NA,point_color=gene_colors)+ggtitle("PCA of residuals of genes from network estimate without ngeneson")

library(stringr)
unstim<-colnames(residuals.nong)%like%"Unstim"

S1<-cor(t(residuals.nong[,unstim]))
S1.s<-cor(t(residuals.nong[,!unstim]))

S2<-cor(t(residuals[,unstim]))
S2.s<-cor(t(residuals[,!unstim]))

extractCor<-function(x,genes,class,model){
  S<-x[genes,genes]
  S<-S[upper.tri(S,diag=FALSE)]
  gn<-matrix(rep(genes,length(genes)),ncol=length(genes),byrow=TRUE)
  gn1<-gn[upper.tri(gn,diag=FALSE)]
  gn<-matrix(rep(genes,length(genes)),ncol=length(genes),byrow=FALSE)
  gn2<-gn[upper.tri(gn,diag=FALSE)]
  data.table(rho=S,gene1=gn1,gene2=gn2,class,model)
}

C<-extractCor(S1,rownames(S1),class = "no_ngeneson",model = "no_ngeneson")
C.s<-extractCor(S1.s,rownames(S1.s),class = "no_ngeneson",model = "no_ngeneson")
C2<-extractCor(S2,rownames(S2),class = "ngeneson",model = "ngeneson")
C2.s<-extractCor(S2.s,rownames(S2.s),class = "ngeneson",model = "ngeneson")

C3<-rbind(C,C2)
C3.s<-rbind(C.s,C2.s)

stopifnot(all.equal(C[,list(gene1, gene2)],C2[,list(gene1, gene2)]))
stopifnot(all.equal(C.s[,list(gene1, gene2)],C2.s[,list(gene1, gene2)]))
Cboth <- rbind(data.table("Rho (background)"=C[,rho], "Rho (CDR)"=C2[,rho], Condition='Non-stimulated'),
               data.table("Rho (background)"=C.s[,rho], "Rho (CDR)"=C2.s[,rho], Condition='Stimulated'))
Cboth <- Cboth[sample(nrow(Cboth), 1e5),]



C3 = C3[,delrho := rho[class == "no_ngeneson"] - rho[class == "ngeneson"],list(gene1,gene2)]
C3.s[,delrho := rho[class == "no_ngeneson"] - rho[class == "ngeneson"],list(gene1,gene2)]
setkey(C3,delrho)
setkey(C3.s,delrho)

C3[,bin:=cut(delrho,breaks = 1000)]
C3 = C3[!is.na(bin)]
C3sumry=C3[,.N,bin]
C3sumry[,cdf:=cumsum(N)/sum(N)]
C3sumry[,bin_c:=mean(as.numeric(str_split_fixed(pattern=",",gsub("]","",gsub("\\(","",bin)),2))),bin]

# p = ggplot(C3sumry) + geom_line(aes(x = bin_c,y = cdf,group = 1)) + geom_vline(x =
#                                                                                  0) + theme_linedraw() + theme(axis.text.x = element_text(angle = 90,hjust =
#                                                                                  1),aspect.ratio = 1) + scale_x_continuous("Difference in correlation (without CDR - with CDR)",breaks =
#                                                                                  seq(-1,1,l = 21)) + scale_y_continuous("fraction of gene pairs")+ggtitle("Unstimulated")
#  

C3.s[,bin:=cut(delrho,breaks = 1000)]
C3.s = C3.s[!is.na(bin)]
C3s.sumry=C3.s[,.N,bin]
C3s.sumry[,cdf:=cumsum(N)/sum(N)]
C3s.sumry[,bin_c:=mean(as.numeric(str_split_fixed(pattern=",",gsub("]","",gsub("\\(","",bin)),2))),bin]

p=ggplot(data.table(rbind(C3s.sumry,C3sumry),Stim=rep(c("Stim","Unstim"),c(nrow(C3s.sumry),nrow(C3sumry)))))+geom_line(aes(x=bin_c,y=cdf,col=Stim))+ geom_vline(x =
                                                                                 0) + theme_linedraw() + theme(axis.text.x = element_text(angle = 90,hjust =
                                                                                 1),aspect.ratio = 1) + scale_x_continuous("Difference in correlation (without CDR - with CDR)",breaks =
                                                                                 seq(-1,1,l = 21)) + scale_y_continuous("fraction of gene pairs")
# 
# p.s = ggplot(C3s.sumry) + geom_line(aes(x = bin_c,y = cdf,group = 1)) + geom_vline(x =
#                                                                                  0) + theme_linedraw() + theme(axis.text.x = element_text(angle = 90,hjust =
#                                                                                  1),aspect.ratio = 1) + scale_x_continuous("Difference in correlation (without CDR - with CDR)",breaks =
#                                                                                  seq(-1,1,l = 21)) + scale_y_continuous("fraction of gene pairs")+ggtitle("Stimulated")
#  
#                                                                          
# p<-ggplot(C3)+aes(x=delrho)+stat_ecdf()+theme_linedraw()+ggtitle("Unstimulated")+theme(aspect.ratio=1)
# p.s<-ggplot(C3.s)+aes(x=delrho)+stat_ecdf()+theme_linedraw()+ggtitle("Stimulated")+theme(aspect.ratio=1)

library(gridExtra)

colors<-c(brewer.pal(name = "PiYG",n=3)[c(1)],"#5729EE")
layout<-matrix(1,ncol=1)
pdf("../inst/extdata/output/Supplementary_Figure_10.pdf",width=8,height=4)
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow=1,ncol=1)))
print(p+plotheme+scale_color_discrete("Condition",labels=c("Stimulated","Non-stimulated")),vp=viewport(layout.pos.row = 1,layout.pos.col = 1))
#print(p.s+plotheme+scale_color_discrete("Model",labels=c("CDR adjusted","No CDR")),vp=viewport(layout.pos.row = 1,layout.pos.col = 2))
dev.off()

cond_de_genes_ng<-res_gene_hurdle[adj<0.01&abs(logFC)>FCTHRESHOLD,symbolid]
cond_de_genes_nong<-res_gene_hurdle_nong[adj<0.01&abs(logFC)>FCTHRESHOLD,symbolid]
foo<-data.table(S2[cond_de_genes_ng,cond_de_genes_ng][upper.tri(S2[cond_de_genes_ng,cond_de_genes_ng],diag=FALSE)])
ng.signif.cors.unstim<-nrow(foo)*foo[,mean(p.adjust(r.test(n=sum(unstim),r12=V1)$p,"fdr")<0.01)]
foo<-data.table(S2.s[cond_de_genes_ng,cond_de_genes_ng][upper.tri(S2.s[cond_de_genes_ng,cond_de_genes_ng],diag=FALSE)])
ng.signif.cors.stim<-nrow(foo)*foo[,mean(p.adjust(r.test(n=sum(!unstim),r12=V1)$p,"fdr")<0.01)]

foo<-data.table(S1[cond_de_genes_nong,cond_de_genes_nong][upper.tri(S1[cond_de_genes_nong,cond_de_genes_nong],diag=FALSE)])
nong.signif.cors.unstim<-nrow(foo)*foo[,mean(p.adjust(r.test(n=sum(unstim),r12=V1)$p,"fdr")<0.01,na.rm=TRUE)]
foo<-data.table(S1.s[cond_de_genes_nong,cond_de_genes_nong][upper.tri(S1.s[cond_de_genes_nong,cond_de_genes_nong],diag=FALSE)])
nong.signif.cors.stim<-nrow(foo)*foo[,mean(p.adjust(r.test(n=sum(!unstim),r12=V1)$p,"fdr")<0.01)]


C<-extractCor(S1[cond_de_genes_nong,cond_de_genes_nong],rownames(S1[cond_de_genes_nong,cond_de_genes_nong]),class = "no_ngeneson",model = "no_ngeneson")
C.s<-extractCor(S1.s[cond_de_genes_nong,cond_de_genes_nong],rownames(S1.s[cond_de_genes_nong,cond_de_genes_nong]),class = "no_ngeneson",model = "no_ngeneson")
C2<-extractCor(S2[cond_de_genes_ng,cond_de_genes_ng],rownames(S2[cond_de_genes_ng,cond_de_genes_ng]),class = "ngeneson",model = "ngeneson")
C2.s<-extractCor(S2.s[cond_de_genes_ng,cond_de_genes_ng],rownames(S2.s[cond_de_genes_ng,cond_de_genes_ng]),class = "ngeneson",model = "ngeneson")


# #################
highly.correlated.genes.ng.unstim<-unique(unlist(C2[data.table(S2[cond_de_genes_ng,cond_de_genes_ng][upper.tri(S2[cond_de_genes_ng,cond_de_genes_ng],diag=FALSE)])[,which((p.adjust(r.test(n=sum(unstim),r12=V1)$p,"fdr")<0.01))],list(gene1,gene2)],use.names=FALSE))
highly.correlated.genes.ng.stim<-unique(unlist(C2.s[data.table(S2.s[cond_de_genes_ng,cond_de_genes_ng][upper.tri(S2.s[cond_de_genes_ng,cond_de_genes_ng],diag=FALSE)])[,which((p.adjust(r.test(n=sum(!unstim),r12=V1)$p,"fdr")<0.01))],list(gene1,gene2)],use.names=FALSE))
highly.correlated.genes.nong.unstim<-unique(unlist(C[data.table(S1[cond_de_genes_nong,cond_de_genes_nong][upper.tri(S1[cond_de_genes_nong,cond_de_genes_nong],diag=FALSE)])[,which((p.adjust(r.test(n=sum(unstim),r12=V1)$p,"fdr")<0.01))],list(gene1,gene2)],use.names=FALSE))
highly.correlated.genes.nong.stim<-unique(unlist(C.s[data.table(S1.s[cond_de_genes_nong,cond_de_genes_nong][upper.tri(S1.s[cond_de_genes_nong,cond_de_genes_nong],diag=FALSE)])[,which((p.adjust(r.test(n=sum(!unstim),r12=V1)$p,"fdr")<0.01))],list(gene1,gene2)],use.names=FALSE))

sum(highly.correlated.genes.ng.stim%in%c(mait_signature_genes_up,mait_signature_genes_down))
sum(highly.correlated.genes.nong.stim%in%c(mait_signature_genes_up,mait_signature_genes_down))

sum(highly.correlated.genes.ng.unstim%in%c(mait_signature_genes_up,mait_signature_genes_down))
sum(highly.correlated.genes.nong.unstim%in%c(mait_signature_genes_up,mait_signature_genes_down))

Supplementary Figure 10 - Change in residual background correlation after adjusting for CDR.

grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow=1,ncol=1)))
print(p+plotheme+scale_color_discrete("Condition",labels=c("Stimulated","Non-Stimulated")),vp=viewport(layout.pos.row = 1,layout.pos.col = 1))
#print(p.s+plotheme+scale_color_discrete("Model",labels=c("CDR adjusted","No CDR")),vp=viewport(layout.pos.row = 1,layout.pos.col = 2))

library(class)
genes_to_plot <- res_gene_hurdle[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),adj)][1:100,symbolid]
qplot(x=GZMB,y=IFNG,data=(data.frame(cData(sca),exprs(subset(sca[,c("IFNG","GZMB")],condition%in%c("Stim","Unstim"))))),col=condition)+theme_linedraw()+geom_hline(y=12,lty=2)+geom_vline(x=10,lty=2)
#qplot(x=IFNG,y=SIN3A,data=(data.frame(cData(sca_stim),exprs((sca_stim[,c("SIN3A","RREB1","CD8A","IFNG","GZMB","PIM2","PIM1","NCOR1")])))),col=block)+theme_linedraw()

#IFNG < 12
#GZMB <10
blocked_cells<-data.table(cData(subset(sca,condition%in%"Stim")),exprs(subset(sca[,c("IFNG","GZMB")],condition%in%"Stim")))
blocked_cells[,block:=GZMB<10&IFNG < 12]
sca_stim<-subset(sca,condition%in%"Stim")
cData(sca_stim)<-merge(cData(sca_stim),blocked_cells[,c("wellKey","block"),with=FALSE],by="wellKey",all=TRUE)
#fit zlm
zlm_block <- zlm.SingleCellAssay(~block + cngeneson, sca_stim, method = "bayesglm", ebayes = TRUE, ebayesControl = list(method = "MLE", model = "H1"),hook  =deviance_residuals_hook)
block_logfc<-getLogFC(zlm_block)
block_logfc<-na.omit(block_logfc[contrast=="blockTRUE"])
block_logfc[order(-abs(logFC),2*(1-pnorm(abs(z)))),list(primerid,contrast,logFC,varLogFC,z,p=2*(1-pnorm(abs(z))),adj=p.adjust(2*(1-pnorm(abs(z))),"fdr"))][adj<0.05&abs(logFC)>FCTHRESHOLD]

lr_block<-lrTest(object = zlm_block,CoefficientHypothesis("blockTRUE"))
lr_block<-dcast(melt(lr_block[,"hurdle",]),primerid~metric)
lr_block<-data.table(lr_block)
lr_block[,adj:=p.adjust(`Pr(>Chisq)`,"fdr")]
lr_block[adj<0.05][order(adj)]

res_gene_block<-merge(lr_block,block_logfc,by="primerid")
genes_to_plot_block<-res_gene_block[adj<0.05&abs(logFC)>log2(2)][order(-abs(logFC),adj)][,primerid]
colors<-colorRampPalette(brewer.pal(n=9,"PiYG"))(21)
colors[11]<-"#FFFFFF"
COMPASS:::pheatmap(exprs(sca_stim[,genes_to_plot_block]),scale="column",col=colors,breaks=seq(-4,4,l=21))
hmap_stim_nonresponders <- grid.grab(wrap=TRUE)

library(class)
genes_to_plot <- res_gene_hurdle[abs(logFC)>FCTHRESHOLD&adj<0.01][order(-abs(logFC),adj)][1:100,symbolid]
qplot(x=GZMB,y=IFNG,data=(data.frame(cData(sca),exprs(subset(sca[,c("IFNG","GZMB")],condition%in%c("Stim","Unstim"))))),col=condition)+theme_linedraw()+geom_hline(y=12,lty=2)+geom_vline(x=10,lty=2)
#qplot(x=IFNG,y=SIN3A,data=(data.frame(cData(sca_stim),exprs((sca_stim[,c("SIN3A","RREB1","CD8A","IFNG","GZMB","PIM2","PIM1","NCOR1")])))),col=block)+theme_linedraw()

#IFNG < 12
#GZMB <10
blocked_cells<-data.table(cData(subset(sca,condition%in%c("Stim","Unstim"))),exprs(subset(sca[,c("IFNG","GZMB")],condition%in%c("Stim","Unstim"))))
blocked_cells[,block:=condition%in%"Stim"&GZMB<10&IFNG < 12]
sca_stim<-subset(sca,condition%in%"Stim")
cData(sca)<-merge(cData(sca),blocked_cells[,c("block","wellKey"),with=FALSE],by="wellKey",all=TRUE)
sca_noact<-subset(sca,condition%in%"Unstim"|block==TRUE)
#fit zlm
zlm_noact <- zlm.SingleCellAssay(~block + cngeneson, sca_noact, method = "bayesglm", ebayes = TRUE, ebayesControl = list(method = "MLE", model = "H1"),hook  =deviance_residuals_hook)
noact_logfc<-getLogFC(zlm_noact)
noact_logfc<-na.omit(noact_logfc[contrast=="blockTRUE"])
noact_logfc[order(-abs(logFC),2*(1-pnorm(abs(z)))),list(primerid,contrast,logFC,varLogFC,z,p=2*(1-pnorm(abs(z))),adj=p.adjust(2*(1-pnorm(abs(z))),"fdr"))][adj<0.05&abs(logFC)>FCTHRESHOLD]

lr_noact<-lrTest(object = zlm_noact,CoefficientHypothesis("blockTRUE"))
lr_noact<-dcast(melt(lr_noact[,"hurdle",]),primerid~metric)
lr_noact<-data.table(lr_noact)
lr_noact[,adj:=p.adjust(`Pr(>Chisq)`,"fdr")]
lr_noact[adj<0.05][order(adj)]

res_gene_noact<-merge(lr_noact,noact_logfc,by="primerid")
genes_to_plot_noact<-res_gene_noact[adj<0.15&abs(logFC)>log2(2)][order(-abs(logFC),adj)][,primerid]
colors<-colorRampPalette(brewer.pal(n=9,"PiYG"))(21)
colors[11]<-"#FFFFFF"
COMPASS:::pheatmap(exprs(sca_noact[,c(genes_to_plot_noact)]),scale="column",col=colors,breaks=seq(-4,4,l=21))
hmap_unstim_nonresponders<-grid.grab(wrap=TRUE)

nonresponders_biplot<-PCbiplot(prcomp(exprs(sca[,c(genes_to_plot_noact,genes_to_plot_block)])),x="PC1",y="PC2",point_colors = c("#E9A3C9","#5729EE")[cData(sca)$block+1],gene_subset=c(1:6,26,29:34))+theme_linedraw()+theme(aspect.ratio=1)+plotheme

pdf("../inst/extdata/output/Supplementary_Figure_14.pdf",width=24,height=7)
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow=1,ncol=3)))
vp<-viewport(layout.pos.row = 1,layout.pos.col = 1)
pushViewport(vp)
grid.draw(hmap_stim_nonresponders)
popViewport(n=2)
vp<-viewport(layout.pos.row = 1,layout.pos.col = 2)
pushViewport(vp)
grid.draw(hmap_unstim_nonresponders)
popViewport(n=2)
vp<-viewport(layout.pos.row = 1,layout.pos.col = 3)
pushViewport(vp)
print(nonresponders_biplot,vp=vp)
popViewport(n=2)
dev.off()

rowann<-cData(sca)[,"block",FALSE]
rowann$block<-factor(rowann$block)
rownames(rowann)<-rownames(exprs(sca[,c(genes_to_plot_noact,genes_to_plot_block)]))
#COMPASS:::pheatmap(exprs(sca[,unique(c(genes_to_plot_noact,genes_to_plot_block))]),row_annotation=rowann[,"block",FALSE])

library(org.Hs.eg.db)
symbols<-intersect(mappedkeys(org.Hs.egSYMBOL2EG),unique(c(genes_to_plot_noact,genes_to_plot_block)))
symbols<-melt(as.list(org.Hs.egSYMBOL2EG[symbols]))
setnames(symbols,c("egid","symbol"))
map<-melt(as.list(org.Hs.egMAP[as.character(symbols$egid)]))
setnames(map,c("chrmap","egid"))
symbols<-merge(map,symbols)
symbols<-data.table(symbols)
kable(symbols[chrmap%in%(symbols[,.N,chrmap][N>1]$chrmap)][order(chrmap)])

Supplementary Figure 14 - Heatmaps of DE genes between non-responding stimulated cells and stimulated / non-stimulated cells

grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow=1,ncol=3)))
vp<-viewport(layout.pos.row = 1,layout.pos.col = 1)
pushViewport(vp)
grid.draw(hmap_stim_nonresponders)
popViewport(n=2)
vp<-viewport(layout.pos.row = 1,layout.pos.col = 2)
pushViewport(vp)
grid.draw(hmap_unstim_nonresponders)
popViewport(n=2)
vp<-viewport(layout.pos.row = 1,layout.pos.col = 3)
pushViewport(vp)
print(nonresponders_biplot,vp=vp)
popViewport(n=2)

GSEA between non-stimulated and non-responding cells

module<-"bp"
module_file <- list.files("../inst/extdata/", pattern = module, full.names = TRUE)
gene_set <- getGmt(module_file)
gene_ids <- geneIds(gene_set)
sets_indices <- ids2indices(gene_ids, fd$symbolid)
# Only keep modules with at least min_gene_in_module
sets_indices <- sets_indices[sapply(sets_indices, length) >= min_gene_in_module]

cache(boot_noact_vs_unstim<-bootVcov1(zlm_noact,R=50))
cache(boot_noact_vs_stim<-bootVcov1(zlm_block,R=50))
cache(gsea_noact_unstim<-gseaAfterBoot(zFit = zlm_noact,boots=boot_noact_vs_unstim,hypothesis = CoefficientHypothesis("blockTRUE"),sets=sets_indices))
cache(gsea_noact_stim<-gseaAfterBoot(zFit = zlm_block,boots=boot_noact_vs_stim,hypothesis = CoefficientHypothesis("blockTRUE"),sets=sets_indices))

t_stat_unstim <- melt(calcZ(gsea_noact_unstim, testType = "normal"))
t_stat_unstim_wide <- dcast(t_stat_unstim, set ~ comp + metric)
t_stat_unstim <- data.frame(calcZ(gsea_noact_unstim, testType = "normal", combined = "stouffer"))
t_stat_unstim$set <- rownames(t_stat_unstim)
t_stat_unstim_comb <- data.table(merge(t_stat_unstim_wide, t_stat_unstim, by = "set"))
setorder(t_stat_unstim_comb,P)
t_stat_unstim_comb[p.adjust(disc_P,"fdr")<0.01, ][order(disc_Z)]
filt_unstim_filt=unique(c(names(which(abs(gsea_noact_unstim[,"disc","stat","test"]-gsea_noact_unstim[,"disc","stat","null"])> log(2.5) )),names(which(abs(gsea_noact_unstim[,"cont","stat","test"]-gsea_noact_unstim[,"cont","stat","null"])>log2(1.5)))))
t_stat_unstim_comb[set%in%filt_unstim_filt&p.adjust(disc_P,"fdr")<0.01][order(disc_Z)]

t_stat_stim <- melt(calcZ(gsea_noact_stim, testType = "normal"))
t_stat_stim_wide <- dcast(t_stat_stim, set ~ comp + metric)
t_stat_stim <- data.frame(calcZ(gsea_noact_stim, testType = "normal", combined = "stouffer"))
t_stat_stim$set <- rownames(t_stat_stim)
t_stat_stim_comb <- data.table(merge(t_stat_stim_wide, t_stat_stim, by = "set"))
setorder(t_stat_stim_comb,P)
filt_stim_filt=unique(c(names(which(abs(gsea_noact_stim[,"disc","stat","test"]-gsea_noact_stim[,"disc","stat","null"])> log(2.5) )),names(which(abs(gsea_noact_stim[,"cont","stat","test"]-gsea_noact_stim[,"cont","stat","null"])>log2(1.5)))))
t_stat_stim_comb[set%in%filt_stim_filt&p.adjust(disc_P,"fdr")<0.01][order(disc_Z)]

lim=t_stat_stim_comb[set%in%filt_stim_filt&p.adjust(P,"fdr")<0.01][order(Z),set]

ggplot(setDT(melt(t_stat_stim_comb[set%in%filt_stim_filt&p.adjust(P,"fdr")<0.01][order(Z),]))[variable%like%"Z"])+geom_tile(aes(x=variable,y=set,fill=value))+scale_y_discrete(limits=lim)+scale_fill_gradient2("GSEA Z-score",low = ("#E9A3C9"),mid=("#F7F7F7"),high=("#A1D76A"),space="Lab")+scale_x_discrete("Component",labels=c("Continous","Discrete","Combined"))+theme_linedraw()+plotheme+theme(axis.text.x=element_text(angle=90,hjust=1))+ggtitle("GSEA of non-responding MAITs vs stimulated MAITs")
ggsave(file="../inst/extdata/output/GSEA_nonresponsive_vs_stim_gobp.pdf",width=15,height=10)
lim2=t_stat_unstim_comb[set%in%filt_unstim_filt&p.adjust(P,"fdr")<0.01][order(Z),set]

ggplot(setDT(melt(t_stat_unstim_comb[set%in%filt_unstim_filt&p.adjust(P,"fdr")<0.01][order(Z)]))[variable%like%"Z"])+geom_tile(aes(x=variable,y=set,fill=value))+scale_y_discrete(limits=lim2)+scale_fill_gradient2("GSEA Z-score",low = ("#E9A3C9"),mid=("#F7F7F7"),high=("#A1D76A"),space="Lab")+scale_x_discrete("Component",labels=c("Continous","Discrete","Combined"))+theme_linedraw()+plotheme+theme(axis.text.x=element_text(angle=90,hjust=1))+ggtitle("GSEA of non-responding MAITs vs unstimulated MAITs")
ggsave(file="../inst/extdata/output/GSEA_nonresponsive_vs_unstim_gobp.pdf",width=15,height=10)

layer(sca_mait)<-"tpm"
melted<-data.table(MAST:::melt.SingleCellAssay(sca_mait))
melted[,ngeneson:=mean(tpm>0),wellKey]
melted[,cut_ngeneson:=cut(.SD[,(ngeneson)],quantile(.SD[,(ngeneson)],seq(0,1,length.out=11))),condition]
ggplot(melted[!is.na(cut_ngeneson)][,list(cond_mean=ifelse(sum(tpm)==0,NA_real_,log2(mean(tpm[tpm>0])))),list(cut_ngeneson,primerid,condition)])+geom_boxplot()+aes(x=cut_ngeneson,y=cond_mean)+theme_linedraw()+theme(axis.text.x=element_text(angle=90,hjust=1))+geom_smooth(aes(group=condition))+facet_wrap(~condition)

PCA plots for ngeneson effect

sca=sca_mait
sca=subset(sca,nGeneOn>4000)
pcsm<-prcomp(exprs(sca))$x
scree_mait<-lapply(list(prcomp(exprs(sca))$sdev^2),function(x)x/sum(x))[[1]][1:8]
foo=data.table(MAST:::melt.SingleCellAssay(sca))

Figure 1

mait_pc<-data.frame(pcsm[,1:2],ngeneson=foo[,list(ngeneson=mean(tpm>0)),`wellKey`][,ngeneson],experiment="MAIT",stim=cData(sca)[,"condition"])

pcs<-prcomp(exprs(sca_alex))$x
scree_shalek<-lapply(list(prcomp(exprs(sca_alex))$sdev^2),function(x)x/sum(x))[[1]][1:8]
shalek_pc<-data.frame(pcs[,1:2],ngeneson=cData(sca_alex)[,"ngeneson",FALSE],experiment="DCs",stim=cData(sca_alex)[,"Stim",FALSE])
setnames(shalek_pc,"Stim","stim")
colors2<-brewer.pal(name="PiYG",n=4)
colors<-c(brewer.pal(name = "PiYG",n=3)[c(1)],"#5729EE")

scree_shalek[1:2]*100

f2<-ggplot(data.table(mait_pc))+geom_point(aes(y=PC2,x=ngeneson,col=stim),show_guide=FALSE)+facet_wrap(~experiment,scales="free")+theme_linedraw()+scale_color_manual(values=colors)+scale_y_continuous(paste0("PC2 (",sprintf("%2.1f",100*scree_mait[2]),"%)"))+scale_x_continuous("Cellular Detection Rate")
f1<-ggplot(data.table(mait_pc))+geom_point(aes(y=PC1,x=ngeneson,col=stim),show_guide=FALSE)+facet_wrap(~experiment,scales="free")+theme_linedraw()+scale_color_manual(values=colors)+scale_y_continuous(paste0("PC1 (",sprintf("%2.1f",100*scree_mait[1]),"%)"))+scale_x_continuous("Cellular Detection Rate")
f3<-ggplot(data.table(shalek_pc))+geom_point(aes(y=PC2,x=ngeneson,col=stim),show_guide=FALSE)+facet_wrap(~experiment,scales="free")+theme_linedraw()+scale_color_manual(values=colors2)+scale_y_continuous(paste0("PC2 (",sprintf("%2.1f",100*scree_shalek[2]),"%)"))+scale_x_continuous("Cellular Detection Rate")
f4<-ggplot(data.table(shalek_pc))+geom_point(aes(y=PC1,x=ngeneson,col=stim),show_guide=FALSE)+facet_wrap(~experiment,scales="free")+theme_linedraw()+scale_color_manual(values=colors2)+scale_y_continuous(paste0("PC1 (",sprintf("%2.1f",100*scree_shalek[1]),"%)"))+scale_x_continuous("Cellular Detection Rate")

scree<-rbind(data.frame(pvar=100*scree_mait,PC=paste0("PC",1:8),experiment="MAIT"),data.frame(pvar=100*scree_shalek,PC=paste0("PC",1:8),experiment="DCs"))

ggplot(scree)+geom_bar(aes(x=PC,y=pvar,fill=experiment),stat="identity",position="dodge")+theme_linedraw()

layout=matrix(1:4,ncol=2)
grid.newpage()
pushViewport(viewport(layout=grid.layout(nrow=nrow(layout),ncol=ncol(layout))))
vp=viewport(layout.pos.row = 1,layout.pos.col = 2)
print(f1,vp=vp)
vp=viewport(layout.pos.row = 1,layout.pos.col = 1)
print(f4,vp=vp)
vp=viewport(layout.pos.row = 2,layout.pos.col = 2)
print(f2,vp=vp)
vp=viewport(layout.pos.row = 2,layout.pos.col = 1)
print(f3,vp=vp)
grid.text("A",x=0.02,y=0.96)
grid.text("B",x=0.52,y=0.96)
grid.text("C",x=0.02,y=0.52)
grid.text("D",x=0.52,y=0.52)

pdf(file="../inst/extdata/output/Figure1.pdf",width=6.5,height=6.5)
grid.newpage()
pushViewport(viewport(layout=grid.layout(nrow=nrow(layout),ncol=ncol(layout))))
vp=viewport(layout.pos.row = 1,layout.pos.col = 2)
print(f1+plotheme,vp=vp)
vp=viewport(layout.pos.row = 1,layout.pos.col = 1)
print(f4+plotheme,vp=vp)
vp=viewport(layout.pos.row = 2,layout.pos.col = 2)
print(f2+plotheme,vp=vp)
vp=viewport(layout.pos.row = 2,layout.pos.col = 1)
print(f3+plotheme,vp=vp)
grid.text("A",x=0.02,y=0.96)
grid.text("B",x=0.52,y=0.96)
grid.text("C",x=0.02,y=0.52)
grid.text("D",x=0.52,y=0.52)
dev.off()

grid.newpage()
pushViewport(viewport(layout=grid.layout(nrow=nrow(layout),ncol=ncol(layout))))
vp=viewport(layout.pos.row = 1,layout.pos.col = 2)
print(f1+plotheme,vp=vp)
vp=viewport(layout.pos.row = 1,layout.pos.col = 1)
print(f4+plotheme,vp=vp)
vp=viewport(layout.pos.row = 2,layout.pos.col = 2)
print(f2+plotheme,vp=vp)
vp=viewport(layout.pos.row = 2,layout.pos.col = 1)
print(f3+plotheme,vp=vp)
grid.text("A",x=0.02,y=0.96)
grid.text("B",x=0.52,y=0.96)
grid.text("C",x=0.02,y=0.52)
grid.text("D",x=0.52,y=0.52)

Flow cytometry data

library(flowCore)
library(flowWorkspace)
library(gridExtra)
ws = openWorkspace("../inst/extdata/(2)140702-BRI MAIT ICS-HM.xml")
g = parseWorkspace(ws,name=4,subset=c(2,3))
nds = getNodes(g)[13:16]
pData(g)$Stim=c("Unstim","Stim")
p1=plotGate(g,"L/Live/MAITs",cond="Stim")
p3=plotGate(g,nds,cond="Stim")
toplot=data.table(fsApply(getData(g,"L/Live/MAITs"),function(x)exprs(x)[,c("FSC-A","SSC-A")]),Stim=rep(c("Unstim","Stim"),fsApply(getData(g,"L/Live/MAITs"),dim)[,1]))
p2<-ggplot(toplot)+aes(x=`FSC-A`,y=`SSC-A`)+geom_point(alpha=0.25)+geom_density2d(aes(col=Stim))+theme_linedraw()+plotheme+ scale_color_manual("",values = rev(c(brewer.pal(name = "PiYG",n = 3)[c(1)],"#5729EE")))
pdf(file="../inst/extdata/output/Supplementary_Figure_6.pdf",width=15,height=5)
grid.arrange(p1,p2,p3,ncol=3)
grid.text(c("A","B","C"),x = c(0.05,0.35,0.65),y=c(0.95,0.95,0.95))
dev.off()
png(file="../inst/extdata/output/Supplementary_Figure_6.png",width=15*150,height=5*150,res=150)
grid.arrange(p1,p2,p3,ncol=3)
grid.text(c("A","B","C"),x = c(0.05,0.35,0.65),y=c(0.95,0.95,0.95))
dev.off()
grid.arrange(p1,p2,p3,ncol=3)
grid.text(c("A","B","C"),x = c(0.05,0.35,0.65),y=c(0.95,0.95,0.95))

RGLab/MASTdata documentation built on May 8, 2019, 5:52 a.m.