Go back to desc homepage

Data Summary

This dataset was generated by Kang et al. (2018) Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nature Biotechnology 36(1):89-94, which can be downloaded in the Gene Expression Omnibus under the accession number GSE96583. The raw gene data, which include count matrix, meta.data( tsne coordinates, the ClusterID in original paper, celltype, and BatchID) were download from Click here.

options(warn = -1)
.libPaths(c("/home/xiaoxiang/R/x86_64-pc-linux-gnu-library/3.5",.libPaths()))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(cowplot))
suppressPackageStartupMessages(library(gridExtra))
suppressPackageStartupMessages(library(png))
suppressPackageStartupMessages(library(RANN))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(kableExtra))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(vioplot))
suppressPackageStartupMessages(library(RColorBrewer))
datadirpath="/media/xiaoxiang/D/raw_data/DECv.1/batchremove_result/"
fig_path="/home/xiaoxiang/Documents/DESC_paper_prepare/DESC_paper_final/formal_revised/figures_sep/"
# get DE genes 
suppressPackageStartupMessages(library(Seurat))# version 2.3.4
datadir="/media/xiaoxiang/D/raw_data/batchcorrect/GSE96583/"
mtx=readMM(paste0(datadir,"matrix.mtx"))
barcodes=read.table(paste0(datadir,"barcodes.tsv"),sep="\t",header=T,stringsAsFactors=F)
meta.data=read.table("GSE96583_cell_info_adata.tsv",sep="\t",header=T,stringsAsFactors=F)
rownames(meta.data)=meta.data$adata_obs_names
gene.info=read.table("GSE96583_gene_info_adata.tsv",sep="\t",header=T,stringsAsFactors=F)
rownames(gene.info)=gene.info$adata_var_names
rownames(mtx)=gene.info$adata_var_names
colnames(mtx)=barcodes$cellname

#
mtx=mtx[,meta.data$adata_obs_names]
cells_id=rep(T,length=dim(mtx)[2])
cells_id[barcodes$use_cells_id=="False"]=F
mtx=mtx[,cells_id]

obj <- CreateSeuratObject(raw.data = mtx,meta.data=meta.data, project = "IMMUNE", min.cells = 10)
#obj=FilterCells(obj,subset.names="nGene",low.thresholds=1,high.thresholds=Inf)
obj=NormalizeData(obj)
obj=ScaleData(obj,display.progress = F)
#
#obj@ident=factor(obj@meta.data$BatchID)
#names(obj@ident)=obj@cell.names
obj@meta.data$BatchID=toupper(obj@meta.data$BatchID)

obj=SetAllIdent(obj,id="BatchID")
#other function FindMarker, FindAllMarkers
###
#data=FindConservedMarkers(obj, ident.1="CTRL", ident.2 = "STIM", grouping.var="celltype",assay.type = "RNA")
#write.table(data,file="DE_CTRL_STIM_groupby_celtype.marker.txt",sep="\t",row.names=F,col.names=T,quote=F)
# within each group, 
# for i,1,2,3,
#
zz=log1p(AverageExpression(obj, show.progress = FALSE,add.ident="celltype"))
zz$genes=rownames(zz)
write.table(zz,file="aveExpression.txt",col.names=T,row.names=F,quote=F,sep="\t")


obj_ctrl=SubsetData(obj,ident.use="CTRL",subset.raw=T)
obj_ctrl=SetAllIdent(obj_ctrl,id="celltype")


obj_stim=SubsetData(obj,ident.use="STIM",subset.raw=T)
obj_stim=SetAllIdent(obj_stim,id="celltype")
#obj_stim@ident=factor(obj_stim@meta.data$celltype)
#names(obj_stim@ident)=obj_stim@cell.names

#obj@ident=factor(obj@meta.data$celltype)
#names(obj@ident)=obj@cell.names
obj=SetAllIdent(obj,id="celltype")

cell_type= names(table(obj@meta.data$celltype))
gde_all_ctrl=data.frame()
for (i in 1:length(cell_type)){
    message(paste("Calculating cluster", cell_type[i]))
    for (j in 1:length(cell_type)){
        gde=FindMarkers(
                  object = obj_ctrl,
                  ident.1 =cell_type[i],
                  ident.2 = cell_type[j],
                  test.use = "MAST")
        if (is.null(x=unlist(gde))){
            next
        }   
            if (nrow(x = gde) > 0) {
            gde=gde[order(gde$p_val,-gde$avg_logFC),]
            gde$clusteri=cell_type[i]
            gde$clusterj=cell_type[j]
                gde$gene <- rownames(x = gde)
                gde_all_ctrl <- rbind(gde_all_ctrl, gde)
            }   
        
    }
}


gde_all_stim=data.frame()
for (i in 1:length(cell_type)){
    message(paste("Calculating 2 cluster", cell_type[i]))
    for (j in 1:length(cell_type)){
        gde=FindMarkers(
                  object = obj_stim,
                  ident.1 =cell_type[i],
                  ident.2 = cell_type[j],
                  test.use = "MAST")
        if (is.null(x=unlist(gde))){
            next
        }   
            if (nrow(x = gde) > 0) {
            gde=gde[order(gde$p_val,-gde$avg_logFC),]
            gde$clusteri=cell_type[i]
            gde$clusterj=cell_type[j]
                gde$gene <- rownames(x = gde)
                gde_all_stim <- rbind(gde_all_stim, gde)
            }   
        
    }
}


gde_cluster_ctrl_stim=data.frame()
obj=SetAllIdent(obj,id="BatchID")
meta.data=obj@meta.data
for (i in 1:length(cell_type)){
    message(paste("Calculating 3 cluster", cell_type[i]))
    obj_tmp=SubsetData(obj,cells.use=rownames(meta.data[meta.data$celltype==cell_type[i],]))
     gde=FindMarkers(
                  object = obj_tmp,
                  ident.1 ="CTRL",
                  ident.2 = "STIM",
                  test.use = "MAST")
         if (is.null(x=unlist(gde))){
                    next
         }
         if (nrow(x = gde) > 0) {
                 gde=gde[order(gde$p_val,-gde$avg_logFC),]
                 gde$cluster=cell_type[i]
                 gde$gene <- rownames(x = gde)
                 gde_cluster_ctrl_stim<- rbind(gde_cluster_ctrl_stim, gde)
         }
}
save(gde_all_ctrl,gde_all_stim,gde_cluster_ctrl_stim,file="resnew.RData")
datadir="/media/xiaoxiang/D/raw_data/DECv.1/batchremove_result"
load(paste0(datadir,"/resnew.RData"))
celltype=c("B cells","CD14+ Monocytes","CD4 T cells","CD8 T cells","Dendritic cells","FCGR3A+ Monocytes","Megakaryocytes","NK cells")
# load necessay
source("/media/xiaoxiang/D/DESC_reproducible_file/helpfunc_new.R")
old=theme_set(theme_bw()+theme(strip.background = element_rect(fill="white"),
                                         panel.background = element_blank(),
                                         panel.grid =element_blank()))
x=gde_all_stim[gde_all_stim$p_val_adj<=0.01,]
x=as.data.frame(table(x$clusteri,x$clusterj))
p_left <- ggplot2::ggplot( data = x, ggplot2::aes( x = Var1,y = Var2, fill = Freq))+ggplot2::theme_bw() +
       ggplot2::geom_tile(color = "grey", size = 0.3) +
       ggplot2::ggtitle("") +
       ggplot2::scale_x_discrete(expand = c(0, 0)) +
       ggplot2::scale_y_discrete(expand = c(0, 0)) +
     geom_text(aes(label=round(Freq)),size=5)+coord_fixed(ratio = 1) +
      ggplot2::scale_fill_gradient(name="# DE genes",
          low = "white",
          high = "steelblue",
      guide = guide_colorbar(title.theme =element_text(size = 12, angle = 0),
                title.vjust = 1.9, barheight = 6, barwidth = 1.6,
                label.theme = ggplot2::element_text(size = 9, angle = 0),
                label.hjust = 0.2)
          )+theme_new
#p_left
x=gde_all_ctrl[gde_all_ctrl$p_val_adj<=0.01,]
x=as.data.frame(table(x$clusteri,x$clusterj))
p_right <- ggplot2::ggplot( data = x, ggplot2::aes( x = Var1,y = Var2, fill = Freq) ) +
       ggplot2::theme_bw() +
    ggplot2::geom_tile(color = "grey", size = 0.3) +
       ggplot2::ggtitle("") +
       ggplot2::scale_x_discrete(expand = c(0, 0)) +
       ggplot2::scale_y_discrete(expand = c(0, 0)) +
     geom_text(aes(label=round(Freq)),size=5)+
       ggplot2::coord_fixed(ratio = 1) +
     ggplot2::scale_fill_gradient(name="# DE genes",
        low = "white",
        high = "steelblue",
            guide = ggplot2::guide_colorbar(title.theme =
                ggplot2::element_text(size = 12, angle = 0),
                title.vjust = 1.9, barheight = 6, barwidth = 1.6,
                label.theme = ggplot2::element_text(size = 9, angle = 0),
                label.hjust = 0.2)) +theme_new
#p_right

Here we only focused on genes with adjusted pvalue <=0.01

x=gde_cluster_ctrl_stim[gde_cluster_ctrl_stim$p_val_adj<=0.01,]
x=as.data.frame(table(x$cluster))
x$x_row="CTRL vs STIM"

p_mid <- ggplot2::ggplot( data = x, ggplot2::aes(y = Var1, x = x_row, fill = Freq) ) +
       ggplot2::theme_bw() +
       ggplot2::geom_tile(aes(width=1),color = "grey", size = 0.3) +
       ggplot2::ggtitle("") +
       ggplot2::scale_x_discrete(expand = c(0, 0)) +
       ggplot2::scale_y_discrete(expand = c(0, 0)) +
     geom_text(aes(label=round(Freq)),size=5)+
       ggplot2::coord_fixed(ratio = 1) +
     ggplot2::scale_fill_gradient(name="# DE genes",
        low = "white",
        high = "steelblue",
            guide = ggplot2::guide_colorbar(title.theme =
                ggplot2::element_text(size = 15, angle = 0),
                title.vjust = 1.9, barheight = 6, barwidth = 1.6,
                label.theme = ggplot2::element_text(size = 9, angle = 0),
                label.hjust = 0.2))+theme_new+
  geom_tile(aes(x=1,y=2),alpha=0,color="red",size=1)
#
#p_mid+coord_flip()+ggplot2::scale_x_discrete(expand = c(0, 0)) +
#       ggplot2::scale_y_discrete(expand = c(0, 0))
gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}
ad=import("anndata",convert = FALSE)
adata_group=ad$read(paste0(datadir,"/GSE96583_bygroup/adata_desc.h5ad"))
adata_group=change_obs(adata_group,exclude = c("n_gene","n_counts","tsne1","tsne2"))
obj=Convert_from_anndata_to_seurat(adata_group,raw.X.slot = "count.data")
obj=NormalizeData(obj,display.progress = F)
obj_group=obj

#adata=ad$read(paste0(datadir,"/GSE96583_bygroup/adata_desc.h5ad"))
adata=ad$read(paste0(datadir,"/GSE96583_prefinal/adata_desc.h5ad"))# nobatch
adata=change_obs(adata,exclude = c("n_gene","n_counts","tsne1","tsne2"))
obj=Convert_from_anndata_to_seurat(adata,raw.X.slot = "count.data")
obj=NormalizeData(obj,display.progress = F)
obj0=obj
tmp=table(obj@meta.data$celltype,obj@meta.data$BatchID)
cellnum=data.frame(CTRL=tmp[,1],STIM=tmp[,2],"Sum"=tmp[,1]+tmp[,2],cluster=names(tmp[,1]))
#cellnum=read.table("cellnum.txt",sep="\t",header=T)
df_cellnum=reshape2::melt(cellnum,id.vars=c("cluster"))
colnames(df_cellnum)=c("Cluster","Condition","Num")
p_num=ggplot2::ggplot( data = df_cellnum, ggplot2::aes( x =Condition ,y = Cluster, fill = Num) ) +
       ggplot2::theme_bw() +
    ggplot2::geom_tile(color = "grey", size = 0.3) +
       ggplot2::ggtitle("") +
       ggplot2::scale_x_discrete(expand = c(0, 0)) +
       ggplot2::scale_y_discrete(expand = c(0, 0)) +
     geom_text(aes(label=round(Num)),size=5)+
       ggplot2::coord_fixed(ratio = 1) +
     ggplot2::scale_fill_gradient(name="# Cells",
        low = "white",
        high = "steelblue",
            guide = ggplot2::guide_colorbar(title.theme =
                ggplot2::element_text(size = 15, angle = 0),
                title.vjust = 1.9, barheight = 6, barwidth = 1.6,
                label.theme = ggplot2::element_text(size = 9, angle = 0),
                label.hjust = 0.2)) +theme_new
p_num=p_num+theme(axis.text.x=element_text(size=15,angle=60),
                  axis.text.y=element_text(size=15),
                  legend.text = element_text(size=15),
                  legend.title = element_text(size=18),
                  plot.title = element_text(size=18))
#save_plot(filename = "~/figure_supp/pbmc_supp/fig.s5.1a.pdf",p_num,base_height = 5,base_width = 10)
#p_num

Supplementary Figure 7

text_tmp=theme(axis.text.x=element_text(size=15),
               axis.text.y=element_text(size=15),
               legend.text = element_text(size=15),
               legend.position="none",
               legend.title = element_text(size=18),
               plot.title = element_text(size=18))

prow=egg::ggarrange(p_left +ggtitle("CTRL")+text_tmp,
              p_mid+ylab("sdfsd")+ggtitle("CTRL vs STIM")+text_tmp+
                theme(axis.text.y=element_blank()), 
              p_right+ggtitle("STIM")+text_tmp+theme(axis.text.y=element_blank(),legend.position = "right"),
              p_num,
              ncol=4,draw=F,widths = c(8,1,8,2.8),
              labels = c("a","","","b"),
              label.args = list(gp = grid::gpar(cex =3.5,color="black",face="bold"),just="bottom"))
obj0=SetAllIdent(obj0,id="BatchID")
avg=log1p(AverageExpression(obj0, show.progress = FALSE,add.ident="celltype"))    
avg$genes=rownames(avg)   
tmp_names=gsub(colnames(avg),pattern = "\\.",replacement = " ")
celltype_tmp=c("B cells","CD4 T cells","CD8 T cells","CD14+ Monocytes","Dendritic cells","FCGR3A+ Monocytes","Megakaryocytes","NK cells","B cells","CD4 T cells","CD8 T cells","CD14+ Monocytes","Dendritic cells","FCGR3A+ Monocytes","Megakaryocytes","NA","NK cells") #ctrl vs stim
pList=rep(list(),length=length(celltype))
names(pList)=celltype
for (i in 1:length(celltype)){
    tmp1=gde_cluster_ctrl_stim[gde_cluster_ctrl_stim$cluster==celltype[i],]
    tmp2=avg[,c(which(celltype_tmp==celltype[i]),ncol(avg))]
    colnames(tmp2)=c("CTRL","STIM","gene")
    tmp=merge(tmp1,tmp2,by="gene",sort=F)
    tmp$color=factor(tmp$p_val_adj<1e-50)
    #pList[[i]]=plot_ctrl_stim(df=tmp,title0 = celltype[i])
    pList[[i]]=ggplot()+geom_point(data=tmp,aes(x=CTRL, y=STIM,color=color),size=2,alpha=0.7)+
      scale_colour_manual(values = c("TRUE"="red","FALSE" = "grey"))+
      ggtitle(celltype[i])+
      theme(legend.position="none",plot.title = element_text(colour = "black"))+
      annotate("text",x=-Inf,y=Inf,label=paste0("# significant genes=",sum(tmp$p_val_adj<=1e-50)),hjust=-0.1,vjust=2,color="blue",fontface="bold",size=8)

}
#plot_grid(plotlist=pList,ncol=4)
pList_new=lapply(pList,function(x) x+theme22(15)+theme(panel.border = element_blank(),
                                                       plot.title = element_text(hjust=0.5,face="bold"),
                                                       axis.line = element_line(size=1,color="black")))
p=egg::ggarrange(plots=pList_new,ncol=4,draw=F)
p0=ggdraw()+
  draw_plot(prow,x=0,y=10.5/18.5,width=1,height=8/18.5)+
  draw_plot(p,x=0,y=0,width =1,height=10.5/18.5)+
  draw_label("c",x=0.01,y=10.5/18.5,hjust=0.5,vjust=1,size=40)#fontface = "bold"
## # Supplementary Figure 7
ggsave(file.path(fig_path,"Figure.S7.tiff"),p0,width=21,height = 18.5,dpi=300,compression="lzw")
ggsave(file.path(fig_path,"Figure.S7.pdf"),p0,width=21,height = 18.5,dpi=300)
p0

#p_adj_val*  >= 50 (that is p_adj_val<=1e-50)
df_fold=gde_cluster_ctrl_stim[,c(1,2,5,6,7)]
colnames(df_fold)=c("pvalue","logFC","p_val_adj","cluster","genes")
df_fold$pvalue=-log(df_fold$pvalue,base=10)
df_fold$p_val_adj=-log(df_fold$p_val_adj,base=10)
# replace Inf by max value non-infiniate
tmp=df_fold$pvalue
tmp[is.infinite(tmp)]=max(tmp[!is.infinite(tmp)])
df_fold$pvalue=tmp

tmp=df_fold$p_val_adj
tmp[is.infinite(tmp)]=max(tmp[!is.infinite(tmp)])
df_fold$p_val_adj=tmp

df_fold1=df_fold[df_fold$logFC>0,]
df_fold1$color=apply(df_fold1,1,function(x) ifelse(as.numeric(x[3])>=50,"blue","grey"))

df_fold2=df_fold[df_fold$logFC<0,]
df_fold2$color=apply(df_fold2,1,function(x) ifelse(as.numeric(x[3])>=50,"red","grey"))

df_fold=rbind(df_fold1,df_fold2)
table(df_fold$cluster,df_fold$color)
##                    
##                     blue grey red
##   B cells              6  364  82
##   CD14+ Monocytes    204  531 281
##   CD4 T cells         21   93 129
##   CD8 T cells          0  256  37
##   Dendritic cells      0  831  28
##   FCGR3A+ Monocytes   14  723  82
##   Megakaryocytes       0  493   0
##   NK cells             1  387  39
celltype=unique(df_fold$cluster)
pList=rep(list(),length=length(celltype))
names(pList)=celltype

label_size=7
for (i in 1:length(pList)){
    df=df_fold[df_fold$cluster==celltype[i],]
    xtab=table(factor(df$color,levels=c("grey","blue","red")))
    p1<-ggplot()
    range0=range(df$pvalue)
    p1<-p1+geom_point(data=df,aes(x=logFC, y=pvalue,color=color),size=2)+ylim(range0[1],range0[2]+25)+scale_colour_manual(values = c("blue" = "blue","grey" = "grey","red" = "red"))
    #subset1=subset(df,color=="red")
    #p1<-p1+geom_text_repel(data=subset1,aes(x=logFC,y=pvalue,label=genes),hjust=0,vjust=0)
    #subset2=subset(df,color=="blue")
    #p1<-p1+geom_text_repel(data=subset2,aes(x=logFC,y=pvalue,label=genes),hjust=0,vjust=0)
    #p1<-p1+labs(x=expression(Log[2]~FoldChange),y=expression(-Log[10]~Pvalue))
    title0=strsplit(celltype[i],"\\+")[[1]]
    p1<-p1+ylab(expression(paste("-", log[10],"(", p,"-",value,")")))+ggtitle(celltype[i])+theme(legend.position="none",plot.title = element_text(colour = "black"))+
      annotate("text",x=-Inf,y=Inf,label=paste0("#Genes=",xtab[names(xtab)=="red"]),hjust=-0.1,vjust=1.3,color="red",fontface="bold",size=label_size)+
      annotate("text",x=Inf,y=Inf,label=paste0("#Genes=",xtab[names(xtab)=="blue"]),hjust=1.0,vjust=1.3,color="blue",fontface="bold",size=label_size)
    #yrange0 <- ggplot_build(p1)$layout$panel_scales_y[[1]]$range$range
  #xrange0 <- ggplot_build(p1)$layout$panel_scales_x[[1]]$range$range
    #p1=p1+annotate("rect", xmin = xrange0[1], xmax =xrange0[1]+0.2*(xrange0[2]-xrange0[1]) , ymin = xrange0[1]+0.8*(yrange0[2]-yrange0[1]), ymax = yrange0[2],alpha = .3)
    pList[[celltype[i]]]=p1
}

#plot_grid(plotlist=pList,ncol=4)

pList_new=lapply(pList,function(x) x+theme22(15)+theme(panel.border = element_blank(),
                                                       plot.title = element_text(hjust=0.5,face="bold"),
                                                       axis.line = element_line(size=1,color="black")))
pList_new_volcano=pList_new
p_volcano=egg::ggarrange(plots=pList_new,ncol=4,draw=F)
pp=egg::ggarrange(plots=lapply(pList_new_volcano,function(x) x+theme22(25)),ncol=8,draw=F)

The Clustering Result

In this section we compared the results obtained from different methods.

markers=list("CD4 T cells"=c("IL7R"),
        "CD14+ Monocytes"=c('CD14',"LYZ"),
        "B cells"=c("MS4A1"),
        "CD8 T cells"=c("CD8A"),
        "FCGR3A+ Monocytes"=c("FCGR3A", "MS4A7"),
        "NK cells"=c("GNLY","NKG7"),
        "Dendritic Cells"=c("FCER1A", "CST3"),
        "Megakaryocytes"=c("PPBP")
)

Note: The above genes are marker genes for each celltype.

DESC result

df0=list()
df0$desc=list()
df0$cca=list()
df0$cca3.0=list()
df0$mnn=list()
df0$scVI=list()
df0$BERMUDA=list()
df0$scanorama=list()
KL_list=list()
KL_list$desc=list()
KL_list$cca=list()
KL_list$cca3.0=list()
KL_list$mnn=list()
KL_list$scvi=list()
KL_list$BERMUDA=list()
KL_list$scanorama=list()
#KL_list$saucie=list()
#if (file.exists("df0.RData")) load("df0.RData")

gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}
resolution_use="0.6"
desc_use=paste0("desc_",resolution_use)
tsne_use=paste0("tsne",resolution_use)
obsm_use=paste0("X_Embeded_z",resolution_use)
markers_update=lapply(markers,function(x) x[x%in%rownames(obj@raw.data)])
df_desc=getplotdata(obj,reduction.use = tsne_use)
df_desc[,desc_use]=factor(as.numeric(df_desc[,desc_use])) # we use resolution=0.8

DESC scaleing across all cells

No batch information provided.

#obj=Convert_from_anndata_to_seurat(adata,raw.X.slot="count.data")
df_desc_nobatch=getplotdata(adata,reduction.use = tsne_use)
df0$desc$df_desc_nobatch=df_desc_nobatch
df0$desc$dr_desc_nobatch=py_to_r(adata_group$obsm[[obsm_use]]) #just for KL
if (! "nobatch" %in% names(KL_list$desc)){
  id14=df_desc_nobatch$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_desc_nobatch,df0$desc$dr_desc_nobatch,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_desc_nobatch[!id14,],df0$desc$dr_desc_nobatch[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="DESC-nobatch"
  KL_list$desc[["nobatch"]]=tmp
}
df_desc_nobatch_feature=getfeature.data(obj0,genes.use=unlist(markers_update),reduction.use = tsne_use)
df_desc_nobatch[,desc_use]=factor(as.numeric(as.character(df_desc_nobatch[,desc_use])))

pdesc_nobatch_1=getplot(df_desc_nobatch,by.group = desc_use,ggtitle0 = "DESC",pt.size = 0.2)+
  theme(plot.title = element_text(color="red"))

pdesc_nobatch_2=getplot(df_desc_nobatch,by.group = "BatchID",ggtitle0 = "DESC",pt.size = 0.2)+scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_desc_nobatch[,'BatchID']))))

pdesc_nobatch_3=getplot(df_desc_nobatch,by.group = "celltype",ggtitle0 = "DESC",pt.size = 0.2)+scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_desc_nobatch[,'celltype']))))

x=cbind(df_desc_nobatch[,c(1,2)],"max.prob"=apply(py_to_r(adata$uns[paste0("prob_matrix",resolution_use)]),1,max))
pdesc_nobatch_4=getfeature.plot(x,pt.size = 0.2,cols.use = c("blue","grey"))[[1]]
pdesc_nobatch_list=getfeature.plot(df_desc_nobatch_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pdesc_nobatch_1,pdesc_nobatch_2,pdesc_nobatch_3,pdesc_nobatch_4,pdesc_nobatch_list[[1]],pdesc_nobatch_list[[2]],ncol=3,draw = F)
p0

DESC scaleing by condition

#obj=Convert_from_anndata_to_seurat(adata,raw.X.slot="count.data")
df_desc=getplotdata(adata_group,reduction.use = tsne_use)
df0$desc=list()
df0$desc$df_desc=df_desc
df0$desc$dr_desc=py_to_r(adata_group$obsm[[obsm_use]]) #just for KL
if (! desc_use %in% names(KL_list$desc)){
  id14=df_desc$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_desc,df0$desc$dr_desc,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_desc[!id14,],df0$desc$dr_desc[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="DESC"
  KL_list$desc[[desc_use]]=tmp
}
df_desc_feature=getfeature.data(obj_group,genes.use=unlist(markers_update),reduction.use = tsne_use)
df_desc[,desc_use]=factor(as.numeric(as.character(df_desc[,desc_use])))
pdesc_1=getplot(df_desc,by.group = desc_use,ggtitle0 = "DESC",pt.size = 0.2)+theme(plot.title = element_text(color="red"))
pdesc_2=getplot(df_desc,by.group = "BatchID",ggtitle0 = "DESC",pt.size = 0.2)+scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_desc[,'BatchID']))))
pdesc_3=getplot(df_desc,by.group = "celltype",ggtitle0 = "DESC",pt.size = 0.2)+scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_desc[,'celltype']))))
x=cbind(df_desc[,c(1,2)],"max.prob"=apply(py_to_r(adata$uns[paste0("prob_matrix",resolution_use)]),1,max))
pdesc_4=getfeature.plot(x,pt.size = 0.2,cols.use = c("blue","grey"))[[1]]
pdesc_list=getfeature.plot(df_desc_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pdesc_1,pdesc_2,pdesc_3,pdesc_4,pdesc_list[[1]],pdesc_list[[2]],ncol=3,draw = F)
p0

CCA3.0

batch_id="BatchID"
detach("package:Seurat",unload = T)
suppressPackageStartupMessages(library("Seurat",lib.loc = "/usr/lib/R/self_library/"))
filename=paste0(datadirpath,"/cca3.0_result/","obj_cca3.0_",batch_id,"_2000.rds")
obj=readRDS(filename)
df_cca3.0=cbind(obj[["tsne"]]@cell.embeddings[,c(1,2)],obj@meta.data)
df0$cca3.0$df_cca3.0=df_cca3.0
var_genes_ori=obj@assays$integrated@var.features
df0$cca3.0$dr_cca3.0=obj[["pca"]]@cell.embeddings[,1:20]
detach("package:Seurat",unload = T)
suppressPackageStartupMessages(library("Seurat"))
#print(mclust::adjustedRandIndex(df_cca$cluster,df_cca$integrated_snn_res.0.2))
#getplot(df_cca,by.group = "sample",pt.size = 0.1)
cca3.0_use="BatchIDtop2000"
if (! cca3.0_use %in% names(KL_list$cca3.0)){
  id14=df_cca3.0$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_cca3.0,df0$cca3.0$dr_cca3.0,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_cca3.0[id14,],df0$cca3.0$dr_cca3.0[id14,],n_cells=100,batch="BatchID")
                 )
  tmp$Method="Seurat3.0"
  KL_list$cca3.0[[cca3.0_use]]=tmp
}
stopifnot(any(rownames(df_desc)==rownames(df_cca3.0)))
df_cca3.0_feature=df_desc_nobatch_feature
df_cca3.0_feature[,c(1,2)]=df_cca3.0[,c(1,2)]
df_cca3.0[,"integrated_snn_res.0.1"]=factor(as.numeric(as.character(df_cca3.0[,"integrated_snn_res.0.1"])))
df_cca3.0[,"res"]=df_cca3.0[,"integrated_snn_res.0.1"]
pcca3.0_1=getplot(df_cca3.0,by.group = "res",ggtitle0 = "Seurat3.0",pt.size = 0.2)+theme(plot.title = element_text(color="red"))
pcca3.0_2=getplot(df_cca3.0,by.group = "BatchID",ggtitle0 = "Seurat3.0",pt.size = 0.1)+scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_cca3.0[,'BatchID']))))
df_cca3.0$celltype[is.na(df_cca3.0$celltype)]="NA"
pcca3.0_3=getplot(df_cca3.0,by.group = "celltype",ggtitle0 = "Seurat3.0",pt.size = 0.1)+scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_cca3.0[,'celltype']))))
pcca3.0_list=getfeature.plot(df_cca3.0_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pcca3.0_1,pcca3.0_2,pcca3.0_3,pcca3.0_list[[1]],pcca3.0_list[[2]],ncol=3,draw = F)
## adding dummy grobs
p0

CCA

batch_id="BatchID"
filename=paste0(datadirpath,"/cca_result/","obj_cca_",batch_id,"_2000_False.rds")
obj=readRDS(filename)
obj=readRDS(filename)
df_cca=cbind(obj@dr$tsne@cell.embeddings[,c(1,2)],obj@meta.data)
df0$cca$df_cca=df_cca
df0$cca$dr_cca=obj@dr$cca@cell.embeddings
#print(mclust::adjustedRandIndex(df_cca$cluster,df_cca$integrated_snn_res.0.2))
#getplot(df_cca,by.group = "sample",pt.size = 0.1)
cca_use="BatchIDtop2000"
if (! cca_use %in% names(KL_list$cca)){
  id14=df_cca$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_cca,df0$cca$dr_cca,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_cca[!id14,],df0$cca$dr_cca[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="CCA"
  KL_list$cca[[cca_use]]=tmp
}
#stopifnot(any(rownames(df_desc)==rownames(df_cca)))
id0=match(rownames(df_desc),df_cca$cellname)
df_cca_feature=df_desc_feature
df_cca_feature[,c(1,2)]=df_cca[id0,c(1,2)]
#df_cca_feature=df_desc_feature
#df_cca_feature[,c(1,2)]=df_cca[,c(1,2)]
df_cca[,"res.0.2"]=factor(as.numeric(as.character(df_cca[,"res.0.2"])))
df_cca[,"res"]=df_cca[,"res.0.2"]
pcca_1=getplot(df_cca,by.group = "res",ggtitle0 = "CCA",pt.size = 0.2)+theme(plot.title = element_text(color="red"))
pcca_2=getplot(df_cca,by.group = "BatchID",ggtitle0 = "CCA",pt.size = 0.1)+scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_cca[,'BatchID']))))
df_cca$celltype[is.na(df_cca$celltype)]="NA"
pcca_3=getplot(df_cca,by.group = "celltype",ggtitle0 = "CCA",pt.size = 0.1)+scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_cca[,'celltype']))))
pcca_list=getfeature.plot(df_cca_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pcca_1,pcca_2,pcca_3,pcca_list[[1]],pcca_list[[2]],ncol=3,draw = F)
## adding dummy grobs
p0

MNN

batch_id="BatchID"
filename=paste0(datadirpath,"/mnn_result/","adata_mnn_",batch_id,"_2000_False.h5ad")
ad=import("anndata",convert = FALSE)
adata=ad$read_h5ad(filename)
adata=change_obs(adata,exclude=c("n_genes","n_counts","sample"))# colum
#print(mclust::adjustedRandIndex(df_cca$cluster,df_cca$integrated_snn_res.0.2))
df_mnn=getplotdata(adata,reduction.use = "tsne")
df0$mnn$df_mnn=df_mnn
df0$mnn$dr_mnn=py_to_r(adata$obsm["X_pca"])[,c(1:20)]
mnn_use="BatchIDtop2000"
if (! mnn_use %in% names(KL_list$mnn)){
  id14=df_mnn$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_mnn,df0$mnn$dr_mnn,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_mnn[!id14,],df0$mnn$dr_mnn[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="MNN"
  KL_list$mnn[[mnn_use]]=tmp
}
#stopifnot(any(rownames(df_desc)==df_mnn$cellname))
id0=match(rownames(df_desc),df_mnn$cellname)
df_mnn_feature=df_desc_feature
df_mnn_feature[,c(1,2)]=df_mnn[id0,c(1,2)]
df_mnn[,"louvain_ori0.3"]=factor(as.numeric(as.character(df_mnn[,"louvain_ori0.3"])))
df_mnn[,"res"]=df_mnn[,"louvain_ori0.3"]
pmnn_1=getplot(df_mnn,by.group = "res",ggtitle0 = "MNN",pt.size = 0.2)+theme(plot.title = element_text(color="red"))
pmnn_2=getplot(df_mnn,by.group = "BatchID",ggtitle0 = "MNN",pt.size = 0.1)+scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_mnn[,'BatchID']))))
pmnn_3=getplot(df_mnn,by.group = "celltype",ggtitle0 = "MNN",pt.size = 0.1)+scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_mnn[,'celltype']))))
pmnn_list=getfeature.plot(df_mnn_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pmnn_1,pmnn_2,pmnn_3,pmnn_list[[1]],pmnn_list[[2]],ncol=3,draw = F)
## adding dummy grobs
p0

scVI

scVI no batch

batch_id="nobatch"
filename=paste0(datadirpath,"/scvi_result/","adata_scvi_",batch_id,"_2000.h5ad")
ad=import("anndata",convert = FALSE)
adata=ad$read_h5ad(filename)
adata=change_obs(adata,exclude=c("n_genes","n_counts","sample"))# colum
df_scvi_nobatch=getplotdata(adata,reduction.use = "tsne")
df0$scvi$df_scvi_nobatch=df_scvi_nobatch
df0$scvi$dr_scvi_nobatch=py_to_r(adata$X)
if (! "nobatch" %in% names(KL_list$scvi)){
  id14=df_scvi_nobatch$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_scvi_nobatch,df0$scvi$dr_scvi_nobatch,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_scvi_nobatch[!id14,],df0$scvi$dr_scvi_nobatch[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="scvi-nobatch"
  KL_list$scvi[["nobatch"]]=tmp
}
stopifnot(any(rownames(df_desc)==df_scvi_nobatch$cellname))
df_scvi_nobatch_feature=df_desc_nobatch_feature
df_scvi_nobatch_feature[,c(1,2)]=df_scvi_nobatch[,c(1,2)]
df_scvi_nobatch[,"louvain_ori0.4"]=factor(as.numeric(as.character(df_scvi_nobatch[,"louvain_ori0.4"])))
df_scvi_nobatch[,"res"]=df_scvi_nobatch[,"louvain_ori0.4"]
pscvi_nobatch_1=getplot(df_scvi_nobatch,by.group = "res",ggtitle0 = "scVI",pt.size = 0.2)+theme(plot.title = element_text(color="red"))
pscvi_nobatch_2=getplot(df_scvi_nobatch,by.group = "BatchID",ggtitle0 = "scVI",pt.size = 0.1)+scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_scvi_nobatch[,'BatchID']))))
pscvi_nobatch_3=getplot(df_scvi_nobatch,by.group = "celltype",ggtitle0 = "scVI",pt.size = 0.1)+scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_scvi_nobatch[,'celltype']))))
pscvi_nobatch_list=getfeature.plot(df_scvi_nobatch_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pscvi_nobatch_1,pscvi_nobatch_2,pscvi_nobatch_3,pscvi_nobatch_list[[1]],pscvi_nobatch_list[[2]],ncol=3,draw = F)
## adding dummy grobs
p0

scVI using batch

batch_id="BatchID"
filename=paste0(datadirpath,"/scvi_result/","adata_scvi_",batch_id,"_2000_False.h5ad")
ad=import("anndata",convert = FALSE)
adata=ad$read_h5ad(filename)
adata=change_obs(adata,exclude=c("n_genes","n_counts","sample"))# colum
df_scvi=getplotdata(adata,reduction.use = "tsne")
df0$scvi$df_scvi=df_scvi
df0$scvi$dr_scvi=py_to_r(adata$X)
scvi_use="BatchIDtop2000"
if (! scvi_use %in% names(KL_list$scvi)){
  id14=df_scvi$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_scvi,df0$scvi$dr_scvi,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_scvi[!id14,],df0$scvi$dr_scvi[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="scvi"
  KL_list$scvi[[scvi_use]]=tmp
}
stopifnot(any(rownames(df_desc)==df_scvi$cellname))
df_scvi_feature=df_desc_feature
df_scvi_feature[,c(1,2)]=df_scvi[,c(1,2)]
df_scvi[,"louvain_ori0.4"]=factor(as.numeric(as.character(df_scvi[,"louvain_ori0.4"])))
df_scvi[,"res"]=df_scvi[,"louvain_ori0.4"]
pscvi_1=getplot(df_scvi,by.group = "res",ggtitle0 = "scVI",pt.size = 0.2)+theme(plot.title = element_text(color="red"))
pscvi_2=getplot(df_scvi,by.group = "BatchID",ggtitle0 = "scVI",pt.size = 0.1)+scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_scvi[,'BatchID']))))
pscvi_3=getplot(df_scvi,by.group = "celltype",ggtitle0 = "scVI",pt.size = 0.1)+scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_scvi[,'celltype']))))
pscvi_list=getfeature.plot(df_scvi_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pscvi_1,pscvi_2,pscvi_3,pscvi_list[[1]],pscvi_list[[2]],ncol=3,draw = F)
## adding dummy grobs
p0

BERMUDA

BERMUDA threshold=0.9

batch_id="BatchID"
filename=paste0(datadirpath,"/BERMUDA_result/adata_bermuda0.9_BatchID.h5ad")
ad=import("anndata",convert = FALSE)
adata=ad$read_h5ad(filename)
adata=change_obs(adata,exclude=c("n_genes","n_counts","sample"))# colum
#print(mclust::adjustedRandIndex(df_cca$cluster,df_cca$integrated_snn_res.0.2))
df_bermuda=getplotdata(adata,reduction.use = "tsne")
df0$BERMUDA$df_bermuda=df_bermuda
df0$BERMUDA$dr_bermuda=py_to_r(adata$obsm["X_bermuda"])[,c(1:20)]
BERMUDA_use="BatchID0.9"
if (! BERMUDA_use %in% names(KL_list$BERMUDA)){
  id14=df_bermuda$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_bermuda,df0$BERMUDA$dr_bermuda,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_bermuda[!id14,],df0$BERMUDA$dr_bermuda[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="BERMUDA"
  KL_list$BERMUDA[[BERMUDA_use]]=tmp
}
#stopifnot(any(rownames(df_desc)==df_mnn$cellname))
id0=match(rownames(df_desc),df_bermuda$cellname)
df_bermuda0.90_feature=df_desc_feature
df_bermuda0.90_feature[,c(1,2)]=df_bermuda[id0,c(1,2)]
df_bermuda[,"louvain_0.4"]=factor(as.numeric(as.character(df_bermuda[,"louvain_0.4"])))
df_bermuda[,"res"]=df_bermuda[,"louvain_0.4"]
df_bermuda_0.90=df_bermuda
pbermuda0.90_1=getplot(df_bermuda,by.group = "res",ggtitle0 = "BERMUDA",pt.size = 0.2)+theme(plot.title = element_text(color="red"))
pbermuda0.90_2=getplot(df_bermuda,by.group = "BatchID",ggtitle0 = "BERMUDA",pt.size = 0.1)+
  scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_bermuda[,'BatchID']))))
pbermuda0.90_3=getplot(df_bermuda,by.group = "celltype",ggtitle0 = "BERMUDA",pt.size = 0.1)+
  scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_bermuda[,'celltype']))))
pbermuda0.90_list=getfeature.plot(df_bermuda0.90_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pbermuda0.90_1,pbermuda0.90_2,pbermuda0.90_3,pbermuda0.90_list[[1]],pbermuda0.90_list[[2]],ncol=3,draw = F)
## adding dummy grobs
p0

BERMUDA threshold=0.85

batch_id="BatchID"
filename=paste0(datadirpath,"/BERMUDA_result/adata_bermuda0.85_BatchID.h5ad")
ad=import("anndata",convert = FALSE)
adata=ad$read_h5ad(filename)
adata=change_obs(adata,exclude=c("n_genes","n_counts","sample"))# colum
#print(mclust::adjustedRandIndex(df_cca$cluster,df_cca$integrated_snn_res.0.2))
df_bermuda=getplotdata(adata,reduction.use = "tsne")
df0$BERMUDA$df_bermuda=df_bermuda
df0$BERMUDA$dr_bermuda=py_to_r(adata$obsm["X_bermuda"])[,c(1:20)]
BERMUDA_use="BatchID0.85"
if (! BERMUDA_use %in% names(KL_list$BERMUDA)){
  id14=df_bermuda$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_bermuda,df0$BERMUDA$dr_bermuda,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_bermuda[!id14,],df0$BERMUDA$dr_bermuda[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="BERMUDA0.85"
  KL_list$BERMUDA[[BERMUDA_use]]=tmp
}
#stopifnot(any(rownames(df_desc)==df_mnn$cellname))
id0=match(rownames(df_desc),df_bermuda$cellname)
df_bermuda0.85_feature=df_desc_feature
df_bermuda0.85_feature[,c(1,2)]=df_bermuda[id0,c(1,2)]
df_bermuda[,"louvain_0.4"]=factor(as.numeric(as.character(df_bermuda[,"louvain_0.4"])))
df_bermuda[,"res"]=df_bermuda[,"louvain_0.4"]
df_bermuda_0.85=df_bermuda
pbermuda0.85_1=getplot(df_bermuda,by.group = "res",ggtitle0 = "BERMUDA0.85",pt.size = 0.2)+theme(plot.title = element_text(color="red"))
pbermuda0.85_2=getplot(df_bermuda,by.group = "BatchID",ggtitle0 = "BERMUDA0.85",pt.size = 0.1)+
  scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_bermuda[,'BatchID']))))
pbermuda0.85_3=getplot(df_bermuda,by.group = "celltype",ggtitle0 = "BERMUDA0.85",pt.size = 0.1)+
  scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_bermuda[,'celltype']))))
pbermuda0.85_list=getfeature.plot(df_bermuda0.85_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pbermuda0.85_1,pbermuda0.85_2,pbermuda0.85_3,pbermuda0.85_list[[1]],pbermuda0.85_list[[2]],ncol=3,draw = F)
## adding dummy grobs
p0

scanorama

Taking scanorama embeded representation as dimension reduction

batch_id="BatchID"
filename=paste0(datadirpath,"/scanorama_result/adata_scanorama_BatchID_2000_False.h5ad")
ad=import("anndata",convert = FALSE)
adata=ad$read_h5ad(filename)
adata=change_obs(adata,exclude=c("n_genes","n_counts","sample"))# colum
#print(mclust::adjustedRandIndex(df_cca$cluster,df_cca$integrated_snn_res.0.2))
df_scanorama=getplotdata(adata,reduction.use = "tsne")
df0$scanorama$df_scanorama=df_scanorama
df0$scanorama$dr_scanorama=py_to_r(adata$obsm["X_scanorama"])[,c(1:100)]
scanorama_use="BatchID"
if (! scanorama_use %in% names(KL_list$scanorama)){
  id14=df_scanorama$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_scanorama,df0$scanorama$dr_scanorama,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_scanorama[!id14,],df0$scanorama$dr_scanorama[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="scanorama"
  KL_list$scanorama[[scanorama_use]]=tmp
}
#stopifnot(any(rownames(df_desc)==df_mnn$cellname))
id0=match(rownames(df_desc),df_scanorama$cellname)
df_scanorama_feature=df_desc_feature
df_scanorama_feature[,c(1,2)]=df_scanorama[id0,c(1,2)]
df_scanorama[,"louvain_0.4"]=factor(as.numeric(as.character(df_scanorama[,"louvain_0.4"])))
df_scanorama[,"res"]=df_scanorama[,"louvain_0.4"]
pscanorama_1=getplot(df_scanorama,by.group = "res",ggtitle0 = "scanorama",pt.size = 0.2)+theme(plot.title = element_text(color="red"))
pscanorama_2=getplot(df_scanorama,by.group = "BatchID",ggtitle0 = "scanorama",pt.size = 0.1)+scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_scanorama[,'BatchID']))))
pscanorama_3=getplot(df_scanorama,by.group = "celltype",ggtitle0 = "scanorama",pt.size = 0.1)+scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_scanorama[,'celltype']))))
pscanorama_list=getfeature.plot(df_scanorama_feature,pt.size = 0.1,cols.use = c("red","grey"))
p0=egg::ggarrange(pscanorama_1,pscanorama_2,pscanorama_3,pscanorama_list[[1]],pscanorama_list[[2]],ncol=3,draw = F)
## adding dummy grobs
p0

Original paper

This show the the results of original paper for Kang et al (2018) 1

df_ori=obj0@meta.data
x=df_ori[,c("tsne1","tsne2")]
colnames(x)=c("tSNE_1","tSNE_2")
df_ori=cbind(x,df_ori)
df0$ori$df_ori=df_ori
#
obj0=ScaleData(obj0,genes.use = var_genes_ori,display.progress = F)
obj0@var.genes=var_genes_ori
obj0=RunPCA(obj0,pcs.compute = 50,do.print = F,genes.print = F)
dr_ori=obj0@dr$pca@cell.embeddings[,1:20]
ori_use="BatchIDtop2000"
if (! ori_use %in% names(KL_list$ori)){
  id14=df_ori$celltype%in%"CD14+ Monocytes"
  tmp=data.frame(BatchID=BatchKL(df_ori,n_cells=100,batch="BatchID"),
                 BatchID_no14=BatchKL(df_ori[!id14,],n_cells=100,batch="BatchID"))
  tmp$Method="ori"
KL_list$ori[[ori_use]]=tmp
}
df_ori_feature=df_desc_feature
df_ori_feature[,c(1,2)]=df_ori[,c(1,2)]
pori_2=getplot(df_ori,by.group = "BatchID",ggtitle0 = "Original Paper",pt.size = 0.1)+scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_ori[,'BatchID']))))
pori_3=getplot(df_ori,by.group = "celltype",ggtitle0 = "Original Paper",pt.size = 0.1)+scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_ori[,'celltype']))))
pori_list=getfeature.plot(df_ori_feature,pt.size = 0.1)
p0=egg::ggarrange(pori_2,pori_3,pori_list[[1]],pori_list[[2]],ncol=3,draw = F)
## adding dummy grobs
p0

Figure 5

theme_use=theme(legend.key.size = unit(1,"cm"),
                legend.key.width = unit(0,'cm'),
                legend.text = element_text(size=16),
                legend.title = element_text(size=20),
                legend.key.height = unit(1.0,"cm"),
                axis.text.x = element_blank(),
                axis.text.y = element_blank())

theme_use2=theme(legend.key.size = unit(0.5,"cm"),
                  legend.text = element_text(size=8),
                  legend.key.height = unit(0.5,"cm"),
                  legend.title = element_text(size=12),
                 #axis.title.x = element_blank(),
                 axis.text.x = element_blank(),
                axis.text.y = element_blank())

theme_use_feature=theme(legend.justification=0.5,
                legend.key.size = unit(1,"cm"),
                legend.text = element_text(size=25),
                legend.title = element_text(size=30),
                legend.key.height = unit(1.3,"cm"),
                axis.text.x = element_blank(),
                axis.text.y=element_blank())
#Figure 2a and #Figure 2b
pp=egg::ggarrange(pdesc_nobatch_2+theme22(25)+theme_use+
                    scale_color_discrete(name="BatchID        ")+
                    theme(legend.title = element_text(size=25),
                          legend.text = element_text(size=20),
                          legend.key.height = unit(1.3,"cm"))+ggtitle("DESC nobatch"),
                  
                  pdesc_nobatch_3+theme22(25)+theme_use+
                        theme(legend.text = element_text(size=20),
                              legend.title = element_text(size=25),
                              legend.key.height = unit(1.0,"cm"))+ggtitle("DESC nobatch"),ncol=2,draw=F,
                  labels = c("a","b"),
                 label.args = list(gp = grid::gpar(cex =3.5,color="black",face="bold"),just="bottom"))
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
#Figure 2c
pp1=egg::ggarrange(plots=lapply(pList_new_volcano,function(x) x+theme22(15)),ncol=4,draw=F,
                  labels = c("c","","","","","","",""),
                 label.args = list(gp = grid::gpar(cex =3.5,color="black",face="bold"),just="bottom"))
data_KL=do.call("rbind",list(KL_list$desc[[1]],
                             KL_list$desc[[2]],
                             KL_list$cca3.0[[1]],
                             KL_list$cca[[1]],
                             KL_list$mnn[[1]],
                             KL_list$scvi[[1]],
                             KL_list$scvi[[2]],
                             KL_list$BERMUDA[[1]],
                             KL_list$BERMUDA[[2]],
                             KL_list$scanorama[[1]],
                             KL_list$ori[[1]]))
data_rect=data.frame(xmin=seq(0.5,9.5,by=1),
                     xmax=seq(1.5,10.5,by=1),
                     ymin=-Inf,
                     color=rep(c("#E8F2F4","#8FB0FF"),length=10),
                     ymax=Inf,stringsAsFactors = F)
#Figure 2d
data_KL=subset(data_KL,!Method%in%c("BERMUDA0.85"))
df=reshape2::melt(data_KL,id.vars="Method")
df$Method=factor(df$Method,levels=c("DESC-nobatch","DESC","Seurat3.0","CCA","MNN","scvi-nobatch","scvi","BERMUDA","scanorama","ori"),
                 labels=c("DESC-nobatch","DESC","Seurat3.0","CCA","MNN","scVI-nobatch","scVI","BERMUDA","scanorama","Original Paper"))
colnames(df)=c("Method","Group","value")
df$Group=factor(df$Group,levels=c("BatchID","BatchID_no14"),labels=c("With CD14+ Monocytes","Without CD14+ Monocytes"))
p_KL=ggplot()+geom_boxplot(data=df,aes(Method,value,color=Group))+
    xlab("")+
   ylab("KL divergence")+
   theme_self+theme(
   #aspect.ratio = 1,
   panel.border =element_blank(),
   axis.line = element_line())+
  scale_x_discrete(expand = c(0,0,0,0))+
  theme(legend.title = element_blank(),
        legend.position=c(0.1,0.9))+
   ggtitle("")+
  theme(legend.position = "right",
                                 legend.justification = "center",
                                 legend.title  = element_blank(),
                                 legend.text = element_text(size=12),
                                 legend.spacing.x = unit(0.5,"cm"),
                                 legend.key.width = unit(1.2,"cm"),
                                 legend.key.height = unit(1.6,"cm"),
                                 axis.text.x   = element_text(angle = 45,hjust=1,size=18,vjust=1))+xlab("")+
                            guides(color=guide_legend(nrow=1,label.position = "bottom",label.theme = element_text(angle=-90,hjust=0,vjust=0.5,size=18)))+
  geom_rect(data=data_rect,aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,fill=color),alpha=0.05)+
    scale_fill_discrete(guide=F)+
  coord_cartesian(clip = "off")
p_KL=egg::ggarrange(p_KL, labels = c("d"),
                 label.args = list(gp = grid::gpar(cex =3.5,color="black",face="bold"),just="bottom"),draw = F)
# Figure 5
p_f5=ggdraw()+
  draw_plot(pp,x=0,y=16/22,width=1,height=6/22)+
  draw_plot(pp1,x=0,y=7/22,width=1,height=9/22)+
  draw_plot(p_KL,x=0,y=0,width=1,height=7/22)
# Figure 5
ggsave(file.path(fig_path,"Figure.5.tiff"),p_f5,width =18,height = 22,dpi=300,compression="lzw")
ggsave(file.path(fig_path,"Figure.5.pdf"),p_f5,width =18,height = 22,dpi=300)
p_f5

pp=egg::ggarrange(pdesc_2+theme22(25)+theme_use+
                    scale_color_discrete(name="BatchID         ")+
                    theme(legend.title = element_text(size=25),
                          legend.text = element_text(size=20),
                          legend.key.height = unit(1.3,"cm"))+ggtitle("DESC"),
                  pdesc_3+theme22(25)+theme_use+
                        theme(legend.text = element_text(size=20),
                              legend.title = element_text(size=25),
                              legend.key.height = unit(1.0,"cm"))+ggtitle("DESC"),ncol=2,draw=F)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Supplementary Figure 8

#supp figure
theme_use_batch=theme_use+theme(legend.text = element_text(size=25),
                                legend.title = element_text(size=30),
                                legend.key.height = unit(1.3,"cm"),
                                axis.text.x = element_blank(),
                                axis.text.y = element_blank(),
                                axis.line.x = element_blank(),
                                axis.line.y=element_blank(),
                                axis.title.y=element_blank(),
                                axis.title.x =element_blank())
theme_use_celltype=theme_use+theme(legend.text = element_text(size=25),
                                   legend.title = element_text(size=30),
                                   legend.key.height = unit(1.0,"cm"),
                                   axis.text.x = element_blank(),
                                   axis.text.y = element_blank(),
                                   axis.line.x = element_blank(),
                                   axis.line.y=element_blank(),
                                   axis.title.y=element_blank(),
                                   axis.title.x =element_blank() )
p_legend=get_legend(pdesc_2+theme(legend.position = "top",
                                          legend.title = element_blank(),
                                          legend.spacing.x = unit(2,"cm"),
                                          legend.text = element_text(size=25))+
                      guides(colour = guide_legend(override.aes = list(alpha = 1, size=15),nrow=1)))
p00=list(pdesc_2+theme(legend.position = "none")+theme_use_batch+ggtitle("DESC"),
                     pcca3.0_2+theme(legend.position = "none")+theme_use_batch,
                     pcca_2+theme(legend.position = "none")+theme_use_batch,
                     pmnn_2+theme(legend.position = "none")+theme_use_batch,
                     pscvi_2+theme(legend.position = "none")+theme_use_batch,
                     pscvi_nobatch_2+theme(legend.position = "none")+theme_use_batch+ggtitle("scVI-nobatch"),
                     pbermuda0.90_2+theme(legend.position = "none")+theme_use_batch,
                     #pbermuda0.85_2+theme(legend.position = "none")+theme_use_batch,
                     pscanorama_2+theme(legend.position = "none")+theme_use_batch,
                  plot_grid(p_legend))
#p01=grid.arrange(grobs=p00,layout_matrix=matrix(1:9,3,3,byrow=T))
p01=egg::ggarrange(plots=p00,ncol = 3,draw = F)
#p01=egg::ggarrange(plots=p00,ncol = 3,draw = F,
#                   labels = c("a","","","","","","","",""),
#                 label.args = list(gp = grid::gpar(cex =3.5,color="black",face="bold"),just="left"))
#p0
p00=egg::ggarrange(pdesc_3+theme(legend.position = "none")+theme_use_batch+ggtitle("DESC"),
                     pcca3.0_3+theme(legend.position = "none")+theme_use_batch,
                     pcca_3+theme(legend.position = "none")+theme_use_batch,
                     pmnn_3+theme(legend.position = "none")+theme_use_batch,
                     pscvi_3+theme(legend.position = "none")+theme_use_batch,
                     pscvi_nobatch_3+theme(legend.position = "none")+theme_use_batch+ggtitle("scVI-nobatch"),
                     pbermuda0.90_3+theme(legend.position = "none")+theme_use_batch,
                     #pbermuda0.85_3+theme(legend.position = "none")+theme_use_batch,
                     pscanorama_3+theme(legend.position = "none")+theme_use_batch,
                   ncol=3,
                   draw=F)
## adding dummy grobs
p_legend=get_legend(pdesc_3+theme(legend.position = "top",
                                          legend.title = element_blank(),
                                          legend.spacing.x = unit(1.3,"cm"),
                                          legend.text = element_text(size=25))+
                      guides(colour = guide_legend(override.aes = list(alpha = 1, size=15),nrow=2)))
p02=grid.arrange(p00,p_legend,ncol=1,heights=c(8,1),newpage = F)

p=ggdraw()+
  draw_label("a",x=0.005,y=0.995,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p01,x=0.005,y=0.54,width=0.995,height=0.45)+
  draw_label("b",x=0.005,y=0.5,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p02,x=0.005,y=0,width=0.995,height=0.5)
# Supplementary Figure 8
ggsave(file.path(fig_path,"Figure.S8.tiff"),p,width =21,height = 30,dpi=300,compression="lzw")
ggsave(file.path(fig_path,"Figure.S8.pdf"),p,width =21,height = 30,dpi=300)
p

Supplementary Figure 9

#save.image(file="all.RData")
#load("all.RData")
theme_use_feature=theme(legend.justification=0.5,
                legend.key.size = unit(1,"cm"),
                legend.text = element_text(size=20),
                legend.title = element_text(size=25),
                legend.key.height = unit(1.0,"cm"),
                axis.text.x = element_blank(),
                axis.text.y=element_blank())
p1=egg::ggarrange(plots=rlist::list.append(lapply(pdesc_nobatch_list,function(x) x+theme_use_feature+coord_cartesian()),
                                                       pdesc_nobatch_1+
                                                 theme_use+
                                                 theme(legend.title = element_blank())+
                                                 ggtitle("DESC-nobatch")),
                      ncol=5,draw=F)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## adding dummy grobs
#
p2=egg::ggarrange(plots=rlist::list.append(lapply(pdesc_list,function(x) x+theme_use_feature+coord_cartesian()),
                                                       pdesc_1+
                                                 theme_use+
                                                 theme(legend.title = element_blank())+
                                                 ggtitle("DESC")), ncol=5,draw=F)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## adding dummy grobs
p=ggdraw()+
  draw_label("a",x=0.005,y=1,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p1,x=0.005,y=0.515,width=0.995,height=0.485)+
  draw_label("b",x=0.005,y=0.485,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p2,x=0.005,y=0,width=0.995,height=0.485)
# Supplementary Figure 8
p

ggsave(file.path(fig_path,"Figure.S9.tiff"),p,width =29,height = 28,dpi=300,compression="lzw")
ggsave(file.path(fig_path,"Figure.S9.pdf"),p,width =29,height = 28,dpi=300)

Supplementary Figure 10

p1=egg::ggarrange(plots=rlist::list.append(lapply(pcca3.0_list,function(x) x+theme_use_feature+coord_cartesian()),
                                                       pcca3.0_1+
                                                 theme_use+
                                                 theme(legend.title = element_blank())+
                                                 ggtitle("Seurat3.0")),
                      ncol=5,draw=F)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## adding dummy grobs
#
p2=egg::ggarrange(plots=rlist::list.append(lapply(pcca_list,function(x) x+theme_use_feature+coord_cartesian()),
                                                       pcca_1+
                                                 theme_use+
                                                 theme(legend.title = element_blank())+
                                                 ggtitle("CCA")),
                      ncol=5,draw=F)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## adding dummy grobs
p=ggdraw()+
  draw_label("a",x=0.005,y=1,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p1,x=0.005,y=0.515,width=0.995,height=0.485)+
  draw_label("b",x=0.005,y=0.485,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p2,x=0.005,y=0,width=0.995,height=0.485)
# Supplementary Figure 10
p

ggsave(file.path(fig_path,"Figure.S10.tiff"),p,width =29,height = 28,dpi=300,compression="lzw")
ggsave(file.path(fig_path,"Figure.S10.pdf"),p,width =29,height = 28,dpi=300)

Supplementary Figure 11

p1=  egg::ggarrange(plots=rlist::list.append(lapply(pmnn_list,function(x) x+theme_use_feature+coord_cartesian()),
                                                       pmnn_1+
                                                 theme_use+
                                                 theme(legend.title = element_blank())+
                                                 ggtitle("MNN")),
                      ncol=5,draw=F)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## adding dummy grobs
#
p2=egg::ggarrange(plots=rlist::list.append(lapply(pscanorama_list,function(x) x+theme_use_feature+coord_cartesian()),
                                                       pscanorama_1+
                                                 theme_use+
                                                 theme(legend.title = element_blank())+
                                                 ggtitle("scanorama")),
                      ncol=5,draw=F)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## adding dummy grobs
p=ggdraw()+
  draw_label("a",x=0.005,y=1,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p1,x=0.005,y=0.515,width=0.995,height=0.485)+
  draw_label("b",x=0.005,y=0.485,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p2,x=0.005,y=0,width=0.995,height=0.485)
# Supplementary Figure 11
p

ggsave(file.path(fig_path,"Figure.S11.tiff"),p,width =29,height = 28,dpi=300,compression="lzw")
ggsave(file.path(fig_path,"Figure.S11.pdf"),p,width =29,height = 28,dpi=300)

Supplementary Figure 12

p1=  egg::ggarrange(plots=rlist::list.append(lapply(pscvi_list,function(x) x+theme_use_feature+coord_cartesian()),
                                                       pscvi_1+
                                                 theme_use+
                                                 theme(legend.title = element_blank())+
                                                 ggtitle("scVI")),
                      ncol=5,draw=F)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## adding dummy grobs
#
p2= egg::ggarrange(plots=rlist::list.append(lapply(pscvi_nobatch_list,function(x) x+theme_use_feature+coord_cartesian()),
                                                       pscvi_nobatch_1+
                                                 theme_use+
                                                 theme(legend.title = element_blank())+
                                                 ggtitle("scVI nobatch")),
                      ncol=5,draw=F)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## adding dummy grobs
p=ggdraw()+
  draw_label("a",x=0.005,y=1,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p1,x=0.005,y=0.515,width=0.995,height=0.485)+
  draw_label("b",x=0.005,y=0.485,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_plot(p2,x=0.005,y=0,width=0.995,height=0.485)
# Supplementary Figure 12
p

ggsave(file.path(fig_path,"Figure.S12.tiff"),p,width =29,height = 28,dpi=300,compression="lzw")
ggsave(file.path(fig_path,"Figure.S12.pdf"),p,width =29,height = 28,dpi=300)

Supplementary Figure 13

p1=  egg::ggarrange(plots=rlist::list.append(lapply(pbermuda0.90_list,function(x) x+theme_use_feature+coord_cartesian()),
                                                       pbermuda0.90_1+
                                                 theme_use+
                                                 theme(legend.title = element_blank())+
                                                 ggtitle("BERMUDA")),
                      ncol=5,draw=F)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## adding dummy grobs
#
pori_2=getplot(df_ori,by.group = "BatchID",ggtitle0 = "Original Paper",pt.size = 1)+
  scale_color_manual(name="BatchID",values=gg_color_hue(length(unique(df_ori[,'BatchID']))))
pori_3=getplot(df_ori,by.group = "celltype",ggtitle0 = "Original Paper",pt.size = 1)+
  scale_color_manual(name="Celltype",values=gg_color_hue(length(unique(df_ori[,'celltype']))))

p2=egg::ggarrange(pori_2+theme22(25)+theme_use+
                    guides(colour = guide_legend(override.aes = list(alpha = 1, size=10)))+
                    theme(legend.title = element_text(size=25),
                          legend.text = element_text(size=20),
                          legend.key.height = unit(1.5,"cm"))+ggtitle("Original Paper"),
                  ggplot()+theme_void(),
                  pori_3+theme22(25)+theme_use+
  guides(colour = guide_legend(override.aes = list(alpha = 1, size=10)))+
                        theme(legend.text = element_text(size=20),
                              legend.title = element_text(size=25),
                              legend.key.height = unit(1.2,"cm"))+ggtitle("Original Paper"),ncol=3,draw=F,widths = c(4,1.1,4))


p=ggdraw()+
  draw_plot(p1,x=0.005,y=11/25,width=0.995,height=14/25)+
  draw_plot(p2,x=0.005,y=0,width=0.995,height=10/25)+
  draw_label("a",x=0.005,y=1,fontface="bold",hjust=0,vjust=1,size=40)+
  draw_label("b",x=0.005,y=10/25,fontface="bold",hjust=0,vjust=1,size=40)
# Supplementary Figure 13
p

ggsave(file.path(fig_path,"Figure.S13.tiff"),p,width =29,height = 25,dpi=300,compression="lzw")
ggsave(file.path(fig_path,"Figure.S13.pdf"),p,width =29,height = 25,dpi=300)

  1. Kang, H.M., Subramaniam, M., Targ, S., Nguyen, M., Maliskova, L., McCarthy, E., Wan, E., Wong, S., Byrnes, L., Lanata, C.M., et al. (2018). Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nat Biotechnol 36, 89-94.↩