---
title: "The results for monocytes dataset"
subtitle: "raw count matrix, CarDEC, DCA+combat, scVI and scanorama"
author:  Xiangjie Li
date: "`r format(Sys.time(), '%m/%d/%Y')`"
output:
  html_notebook:
    number_sections: yes
    toc: yes
  jekyllthat::jekylldown:
  html_document:
    df_print: paged
    toc: yes
    number_sections: yes
  prettydoc::html_pretty:
    theme: cayman
    highlight: github
    math: katex
    toc: yes
---

<style>
pre {
  max-height: 200px;
  float: left;
  width: 910px;
  overflow-y: auto;
}
pre.r {
  max-height: none;
}
</style>


**Data Summary:** 

This dataset was generated by our group, which can be downloaded from [GEO (GSE146974)](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE146974) or or [https://drive.google.com/file/d/1kR8Hhufoo2h2OtomW8n3kM0gaQhVS564/view?usp=sharing](https://drive.google.com/file/d/1kR8Hhufoo2h2OtomW8n3kM0gaQhVS564/view?usp=sharing). This dataset was generated from human peripheral blood mononuclear clear cells by Ficoll Separation followed by CD14 and CD16 positive cell selection. Since the CD14 and CD16 antibodies are not 100% specific, some T cells were also present in the scRNA-seq data. We performed clustering analysis using leiden’s algorithm for each batch and identified 288 T cells in total based on the T cell marker genes CD3D, CD3E and CD3G. Aftering removing these 288 T cells, there are 10,878 cells and 21,289 genes, which was processed and sequenced at three different days, resulting in three batches (3,640 cells in T1, 4,833 cells in T2 and 2,405 cells in T3) left in the remaining analysis. 

__***Human monocyte preparation***__: Monocyte preparation uses a modification of published protocols. Briefly, ~20 ml blood drawn in sodium heparin was processed immediately in the lab in the Clinical Research Center at Columbia University. PBMCs were isolated by gradient Ficoll paque centrifugation, which maintains cell viability and prevents ex vivo activation during cell recovery. Cells were stained with antibodies against human HLADR, CD14 and CD16 and monocyte subsets defined as HLADR+CD14++CD16-(classical), HLADR+CD14++CD16+ (intermediate), HLADR+CD14dim/CD16++ (nonclassical, patrolling monocyte). DAPI staining was used to exclude dead cells. Monocytes were sorted by a BD Influx Sorter into tubes for real-time 10x Genomics analysis.

# Summary 

Here I used monocle3 (monocle3_0.2.1) to conduct the pseudotime analysis. 

> <span style="color:red"> CarDEC, scVI and DCA are both deep learning based methods. For each method, we used all genes as the input, the way of `Using latetn`One is the standard pipline for monocle3 (`denosied count` -> `normalization` -> `scaling` -> `pca` dimension reduction -> `umap visualization` based on `pca` dimension reduction) and the other method replces the `pca` by the latent representation and then umap visualization based on latent representation. </span>


============================================================================================

- __***HVGs raw***__: Only using highly variable genes (HVG) as the inputs for Monocle3
- __***All genes raw***__: Using all genes (HVG) as the inputs for Monocle3
- __***Using latent***__: Replceing the `pca` by the latent representation and then umap visualization based on latent representation (CarDEC, DCA+combat, scVI and scanorama).
- __***HVG denoised***__: Only using denosined expression count from those highly variable genes (HVGs) as the input for monocle3 and then conducted standard pipline of monocle3 (CarDEC, DCA, scVI and scanorama).
- __***All genes denoised***__: Used all denoised expression values as the input for monocle3 and then conducted the standard pipline of monocle3 (CarDEC, DCA, and scVI).

```{r}
options(warn=-1) # turn off warning message globally
.libPaths(c("/home/xiaoxiang/R/x86_64-pc-linux-gnu-library/3.5",.libPaths()))
Sys.setenv(RETICULATE_PYTHON_ENV="/home/xiaoxiang/anaconda3/envs/py36")#="/home/xiaoxiang/.conda/envs/DESCVIR"
Sys.setenv(RETICULATE_PYTHON="/usr/bin/python3")
#RETICULATE_PYTHON="/home/xiaoxiang/anaconda3/bin/python3",
if ("Seurat" %in% loadedNamespaces()) detach("package:Seurat",unload = T)
dyn.load("/home/xiaoxiang/R/x86_64-pc-linux-gnu-library/3.5/sf/libs/sf.so")
#suppressPackageStartupMessages(library(monocle,lib.loc = "/usr/lib/R/monocle_alpha"))# devtools::install_github("")
#devtools::install_github("cole-trapnell-lab/DDRTree", ref="simple-ppt-like",lib="/usr/lib/R/monocle_alpha")
#devtools::install_github("r-spatial/sf") if 
#install.packages("~/Downloads/monocle-release-monocle3_alpha/", repos = NULL,lib = "/usr/lib/R/monocle_alpha")
suppressPackageStartupMessages(library(reticulate))
#suppressPackageStartupMessages(library(devtools))
suppressPackageStartupMessages(library(monocle3))

#suppressPackageStartupMessages(library(flexclust))
#suppressPackageStartupMessages(library(mcclust))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggjoy))
suppressPackageStartupMessages(library(VGAM))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(kableExtra))
suppressPackageStartupMessages(library(cowplot))
#py_install('umap-learn', pip = T, pip_ignore_installed = T)
#import("leiden")
#fig_path="/home/xiaoxiang/Documents/DESC_paper_prepare/DESC_paper_final/formal_revised/figures_sep/"
datadirpath="./"
knitr::opts_chunk$set(echo=T)
```


```{r}
df_pseudotime_list=list()
```

```{r,echo=F,include=T,eval=T}
get_p_new =function(x,x2=NULL) {
 res1=sapply(x,function(i0) ifelse(i0<=0,"<2.2e-16",scales::scientific(i0,digits = 3))) 
 res2=rep(c(""),length=length(res1))
 if(!is.null(x2)){
   res2=round(x2,3) 
 }
 return(paste0(res2," (",res1,")"))
}

num_methods=16#original 13
Stable5=data.frame(Method=c("Raw count HVGs","Raw count All","CarDEC (latent)","CarDEC (denoised HVGs)","CarDEC (denoised All)","scVI (latent)","scVI (denoised HVGs)","scVI (denoised All)","DCA (latent)","DCA (denoised HVGs)","DCA (denoised All)","MNN (denoised HVGs)","MNN (denoised All)","Scanorama (latent)","Scanorama (denoised HVGs)","Scanorama (denoised All)"),
                  `T1 v.s. T2`=rep("<2.2e-16",length=num_methods),
                  `T1 v.s. T3`=rep("<2.2e-16",length=num_methods),
                  `T2 v.s. T3`=rep("<2.2e-16",length=num_methods),stringsAsFactors = F)
rownames(Stable5)=Stable5$Method
#table5 %>%
#  kable() %>%
#  kable_styling()
```


```{r}
# load necessay function
#source("/media/xiaoxiang/D/DESC_reproducible_file/helpfunc_new.R")
#source("/media/xiaoxiang/D/Upenn_computer_backup/Documents/Human_Heart_Project/heart/Heart_result_updated/helpfunc_new.R")
old=theme_set(theme_bw()+theme(strip.background = element_rect(fill="white"),
                                         panel.background = element_blank(),
                               legend.background = element_blank(),
                                         panel.grid =element_blank()))

BatchKL=function(df,dimensionData=NULL,replicates=200,n_neighbors=100,n_cells=100,batch="BatchID"){
  #entropy of batch mixiing
  #replicates is the number of boostrap times
  #n_neighbors is the number of nearest neighbours of cell(from all batchs)
  #n_cells is the number of randomly picked cells
  if (is.null(dimensionData)){
        tsnedata=as.matrix(df[,c("tSNE_1","tSNE_2")])
  }else{
        tsnedata=as.matrix(dimensionData)
  }
  batchdata=factor(as.vector(df[,batch]))
  table.batchdata=as.matrix(table(batchdata))[,1]
  tmp00=table.batchdata/sum(table.batchdata)#proportation of population
  n=dim(df)[1]
  KL=sapply(1:replicates,function(x){
    bootsamples=sample(1:n,n_cells)
    #nearest=nn2(tsnedata,tsnedata[bootsamples,],k=n_neighbors)
    nearest=nabor::knn(tsnedata,tsnedata[bootsamples,],k=min(5*length(tmp00),n_neighbors))
    KL_x=sapply(1:length(bootsamples),function(y){
      id=nearest$nn.idx[y,]
      tmp=as.matrix(table(batchdata[id]))[,1]
      tmp=tmp/sum(tmp)
      return(sum(tmp*log2(tmp/tmp00),na.rm = T))
    })
    return(mean(KL_x,na.rm = T))
  })
  return(KL)
}
```

```{r}
Convert_to_seurat3=function(adata){
  suppressPackageStartupMessages(library("Seurat",lib.loc = "/usr/lib/R/self_library/"))
  mtx=py_to_r(adata$X$T$tocsc())
  cellinfo=py_to_r(adata$obs)
  geneinfo=py_to_r(adata$var)
  colnames(mtx)=cellinfo$cellname
  rownames(mtx)=rownames(geneinfo)
  obj=CreateSeuratObject(mtx,meta.data = cellinfo[,!colnames(cellinfo)%in%c("n_genes","n_counts"),drop=F],min.features  = 1)
  return(obj)
}
getwd()
```


```{r}
get_plot4=function(df00){
  p1=ggplot()+geom_point(data =df00,aes(x=UMAP_1,y=UMAP_2,color=FCGR3A),size=0.01)+
    scale_color_gradient(low="grey",high="red")+
    theme(legend.position = "top")+
    guides(color=guide_colorbar(title.vjust = 0.7))
  
  p2=ggplot()+geom_point(data =df00,aes(x=UMAP_1,y=UMAP_2,color=S100A8),size=0.01)+
    scale_color_gradient(low="grey",high="red")+
    theme(legend.position = "top")+
     guides(color=guide_colorbar(title.vjust = 0.7))
  
  p3=ggplot(data =df00,aes(x=pseudotime,y=FCGR3A))+
      geom_point(aes(color=BatchID),size=0.01)+
      guides(color=guide_legend(override.aes = list(size=5)))+
      geom_smooth(aes(color=BatchID),method="gam",formula = y ~ s(x, bs="cs"))+
       geom_smooth(color="black",method="gam",formula = y ~ s(x, bs="cs"),size=0.5)+
      ggtitle("")+xlab("Pseudotime")+theme(legend.position = "top",
                               plot.title = element_text(size=18,face="bold",hjust=0.5),
                               legend.text = element_text(size=15,face="bold"),
                               plot.margin = unit(c(0,1,0,0),"cm"),
                               legend.title = element_blank())+
    scale_color_brewer(palette = "Set2")
  
  p4=ggplot(data =df00,aes(x=pseudotime,y=S100A8))+
      geom_point(aes(color=BatchID),size=0.01)+
      guides(color=guide_legend(override.aes = list(size=5)))+
      geom_smooth(aes(color=BatchID),method="gam",formula = y ~ s(x, bs="cs"))+
       geom_smooth(color="black",method="gam",formula = y ~ s(x, bs="cs"),size=0.5)+
      ggtitle("")+xlab("Pseudotime")+theme(legend.position = "top",
                               plot.title = element_text(size=18,face="bold",hjust=0.5),
                               legend.text = element_text(size=15,face="bold"),
                               legend.title = element_blank())+scale_color_brewer(palette = "Set2")
  
  p=egg::ggarrange(p1,p3,p2,p4,ncol=4,draw=F)
  return(p)
}

get_plot4_sep=function(df00){
  p1=ggplot()+geom_point(data =df00,aes(x=UMAP_1,y=UMAP_2,color=FCGR3A),size=0.01)+
    scale_color_gradient(low="grey",high="red")+
    theme(legend.position = "top")+
    guides(color=guide_colorbar(title.vjust = 0.7))
  
  p2=ggplot()+geom_point(data =df00,aes(x=UMAP_1,y=UMAP_2,color=S100A8),size=0.01)+
    scale_color_gradient(low="grey",high="red")+
    theme(legend.position = "top")+
     guides(color=guide_colorbar(title.vjust = 0.7))
  
  p3=ggplot(data =df00,aes(x=pseudotime,y=FCGR3A))+
      geom_point(aes(color=BatchID),size=0.05)+
      guides(color=guide_legend(override.aes = list(size=5)))+
      geom_smooth(aes(color=BatchID),method="gam",formula = y ~ s(x, bs="cs"))+
       geom_smooth(color="black",method="gam",formula = y ~ s(x, bs="cs"),size=0.5)+
      ggtitle("")+xlab("Pseudotime")+theme(legend.position = "top",
                               plot.title = element_text(size=18,face="bold",hjust=0.5),
                               legend.text = element_text(size=15,face="bold"),
                               #plot.margin = unit(c(0,1,0,0),"cm"),
                               legend.title = element_blank())+
    scale_color_brewer(palette = "Set2")
  p4=ggplot(data =df00,aes(x=pseudotime,y=S100A8))+
      geom_point(aes(color=BatchID),size=0.05)+
      guides(color=guide_legend(override.aes = list(size=5)))+
      geom_smooth(aes(color=BatchID),method="gam",formula = y ~ s(x, bs="cs"))+
       geom_smooth(color="black",method="gam",formula = y ~ s(x, bs="cs"),size=0.5)+
      ggtitle("")+xlab("Pseudotime")+theme(legend.position = "top",
                               plot.title = element_text(size=18,face="bold",hjust=0.5),
                               legend.text = element_text(size=15,face="bold"),
                               legend.title = element_blank())+
    scale_color_brewer(palette = "Set2")
  return(list(p1,p2,p3,p4))
}

```




```{r}
ad=import("anndata",convert = FALSE)
adata=ad$read_h5ad("../../dca_test.h5ad")
obj0=Convert_to_seurat3(adata)
obj0=NormalizeData(obj0,verbose = F)
raw.data=obj0@assays$RNA@counts
```

```{r}
maprules=c("2017_0801"="T1","2017_1017"="T2","2017_1120"="T3")
maprules
```

Here we compared different methods, including `DCA` and `scVI`. 

```{r}
hvg_genes=read.table("../final_processed_results/CarDEC_hvg_used.tsv",header = T,sep="\t",stringsAsFactors = F)
hvg_genes=subset(hvg_genes,Variance.Type=="HVG") #top 2000 genes 
```

# Monocle3 for `raw data`
## HVGs raw 

```{r}
cell.meta.data=obj0@meta.data
cell.meta.data$dataset_batch=plyr::mapvalues(cell.meta.data$batch_label,names(maprules),maprules)
gene_ann=data.frame(gene_short_name = make.unique(rownames(raw.data)),row.names = make.unique(rownames(raw.data)))
#pd <- new("AnnotatedDataFrame",data=cell.meta.data)
#fd <- new("AnnotatedDataFrame",data=gene_ann)
cds <- new_cell_data_set(raw.data[rownames(raw.data)%in%hvg_genes$genename,], 
                         cell_metadata = cell.meta.data,
                         gene_metadata =gene_ann[gene_ann$gene_short_name%in%hvg_genes$genename,,drop=F])
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log",verbose = F)
## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)
# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, cluster=c("1","5")){
  root_pr_nodes=sapply(cluster,function(ii){
    cell_ids <- which(colData(cds)[, "clusters"] %in%ii)
  
  closest_vertex <-cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
  
  closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
  root_pr_nodes <-igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names(which.max(table(closest_vertex[cell_ids,]))))]
  })
  root_pr_nodes
}
# root cells
ids=get_earliest_principal_node(cds,cluster=c("1","3","4"))
cds <- order_cells(cds,root_pr_nodes = ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```

```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])

set.seed(10)
theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_ori_1=plot_cells(cds,color_cells_by = "dataset_batch",graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")+theme(legend.title = element_blank())
          
p_ori_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_ori_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_ori=egg::ggarrange(p_ori_1,p_ori_2,p_ori_3,ncol=3,draw=F)
```

```{r}
# printed how many cells with no pseudotime
table(as.numeric(is.infinite(pData(cds)[,c("pseudotime")]))) #0 mean normal pseudotime and 1 means infinity.
```

```{r fig.height=5,fig.width=18}
p_monocle_ori
```

```{r}
cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
#cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)/size_factors(cds))))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$raw=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1

tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2

tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[1,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```

## All genes raw 
```{r}
cell.meta.data=obj0@meta.data
cell.meta.data$dataset_batch=plyr::mapvalues(cell.meta.data$batch_label,names(maprules),maprules)
gene_ann=data.frame(gene_short_name = make.unique(rownames(raw.data)),row.names = make.unique(rownames(raw.data)))
#pd <- new("AnnotatedDataFrame",data=cell.meta.data)
#fd <- new("AnnotatedDataFrame",data=gene_ann)
cds <- new_cell_data_set(raw.data, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log",verbose = F)
## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)
# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```


```{r,fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, cluster=c("1","5")){
  root_pr_nodes=sapply(cluster,function(ii){
    cell_ids <- which(colData(cds)[, "clusters"] %in%ii)
  
  closest_vertex <-cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
  
  closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
  root_pr_nodes <-igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names(which.max(table(closest_vertex[cell_ids,]))))]
  })
  root_pr_nodes
}
# root cells
ids=get_earliest_principal_node(cds,cluster=c("4","5"))
cds <- order_cells(cds,root_pr_nodes = ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```

```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])

set.seed(10)
theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_ori_all_1=plot_cells(cds,color_cells_by = "dataset_batch",graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")+theme(legend.title = element_blank())
          
p_ori_all_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_ori_all_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_ori_all=egg::ggarrange(p_ori_all_1,p_ori_all_2,p_ori_all_3,ncol=3,draw=F)
```

```{r}
# printed how many cells with no pseudotime
table(as.numeric(is.infinite(pData(cds)[,c("pseudotime")]))) #0 mean normal pseudotime and 1 means infinity.
```

```{r fig.height=5,fig.width=18}
p_monocle_ori_all
```

```{r}
cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
#cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)/size_factors(cds))))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$raw_all=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r,fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[2,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```

# Monocle3 using carDEC

In this section, we will evalutate the performance of carDEC. 

Note that: carDEC used all genes and extracted HVG to evaluate.

```{r}
adata=ad$read_h5ad("../final_processed_results/CarDEC Results/adata_CarDEC.h5ad")
```

```{r}
cell.meta.data=py_to_r(adata$obs)
cell.meta.data$dataset_batch=plyr::mapvalues(cell.meta.data$batch_label,names(maprules),maprules)
gene_ann0=py_to_r(adata$var)
gene_ann=data.frame(gene_short_name = make.unique(rownames(gene_ann0)),
                    VarianceType=gene_ann0$`Variance Type`,
                    row.names = make.unique(rownames(gene_ann0)))
mtx=t(py_to_r(adata$layers['denoised counts']))
colnames(mtx)=cell.meta.data$cellname
rownames(mtx)=rownames(gene_ann)
mtx_sizefactor=1e4/colSums(mtx)
```


## Using latent 

```{r}
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log",verbose = F)

tmp0=py_to_r(adata$obsm["embedding"])
colnames(tmp0)=paste0("PC",1:ncol(tmp0))
reducedDims(cds)$PCA=tmp0

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("3"))
cds <- order_cells(cds, root_pr_nodes=ids)
plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_carDEC_latent_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_carDEC_latent_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_carDEC_latent_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_carDEC_latent=egg::ggarrange(p_carDEC_latent_1,p_carDEC_latent_2,p_carDEC_latent_3,ncol=3,draw=F)
```


```{r, fig.height=5,fig.width=18}
p_monocle_carDEC_latent
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))=1e4/rowSums(mtx)
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)*mtx_sizefactor)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$carDEC_latent=df0
```


- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[3,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```


## HVGs denoised

```{r}
cds <- new_cell_data_set(mtx[gene_ann$VarianceType=="HVG",], cell_metadata = cell.meta.data,gene_metadata =gene_ann[gene_ann$VarianceType=="HVG",,drop=F])
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log")

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("4"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_carDEC_denoised_hvg_1 = plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_carDEC_denoised_hvg_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_carDEC_denoised_hvg_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_carDEC_denoised_hvg=egg::ggarrange(p_carDEC_denoised_hvg_1,p_carDEC_denoised_hvg_2,p_carDEC_denoised_hvg_3,ncol=3,draw=F)
```


```{r, fig.height=5,fig.width=18}
p_monocle_carDEC_denoised_hvg
```


```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)*mtx_sizefactor)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$carDEC_denoised_hvg=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[4,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```

## All genes denoised

```{r}
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log",verbose = F)

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("4"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```

```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_carDEC_denoised_all_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_carDEC_denoised_all_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_carDEC_denoised_all_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_carDEC_denoised_all=egg::ggarrange(p_carDEC_denoised_all_1,p_carDEC_denoised_all_2,p_carDEC_denoised_all_3,ncol=3,draw=F)
```


```{r, fig.height=5,fig.width=18}
p_monocle_carDEC_denoised_all
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)*mtx_sizefactor)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$carDEC_denoised_all=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[5,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```

# Monocle3 using scVI

```{r}
#adata=ad$read_h5ad("../final_processed_results/scVI Results/monocytes_ALL/adata_all.h5ad")
adata=ad$read_h5ad("../final_processed_results/scVI Results New/monocytes_ALL/adata_all.h5ad")
```

```{r}
cell.meta.data=py_to_r(adata$obs)
cell.meta.data$dataset_batch=plyr::mapvalues(cell.meta.data$batch_label,names(maprules),maprules)
gene_ann0=py_to_r(adata$var)
gene_ann=data.frame(gene_short_name = make.unique(rownames(gene_ann0)),
                    row.names = make.unique(rownames(gene_ann0)))
mtx=t(py_to_r(adata$X))
colnames(mtx)=cell.meta.data$cellname
rownames(mtx)=rownames(gene_ann)
mtx_sizefactor=1e4/colSums(mtx)
```

## Using latent

```{r}
#mtx=mtx[gene_ann$VarianceType=="HVG",]
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log",verbose = F)
tmp0=py_to_r(adata$obsm["X_latent"])
colnames(tmp0)=paste0("PC",1:ncol(tmp0))
reducedDims(cds)$PCA=tmp0

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("3"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
#saveRDS(cds,file = "cds_scvi.rds")
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_scVI_latent_all_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_scVI_latent_all_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_scVI_latent_all_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_scVI_latent_all=egg::ggarrange(p_scVI_latent_all_1,p_scVI_latent_all_2,p_scVI_latent_all_3,ncol=3,draw=F)
```

```{r, fig.height=5,fig.width=18}
p_monocle_scVI_latent_all
```



```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)*mtx_sizefactor)))
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$scVI_latent_all=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[6,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```


## HVGs denoised


```{r}
cds <- new_cell_data_set(mtx[rownames(mtx)%in%hvg_genes$genename,], cell_metadata = cell.meta.data,gene_metadata =gene_ann[gene_ann$gene_short_name%in%hvg_genes$genename,,drop=F])
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log",verbose = F)

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("3"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_scvi_denoised_hvg_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_scvi_denoised_hvg_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_scvi_denoised_hvg_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_scvi_denoised_hvg=egg::ggarrange(p_scvi_denoised_hvg_1,p_scvi_denoised_hvg_2,p_scvi_denoised_hvg_3,ncol=3,draw=F)
```

```{r, fig.height=5,fig.width=18}
p_monocle_scvi_denoised_hvg
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)*mtx_sizefactor)))
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$scVI_denosied_hvg=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[7,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```

## All genes denoised

```{r}
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log",verbose = F)

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("3"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
x0=pseudotime(cds)
x0[is.infinite(x0)]=NA
colData(cds)$pseudotime=x0

colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_scVI_denoised_all_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_scVI_denoised_all_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_scVI_denoised_all_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_scVI_denoised_all=egg::ggarrange(p_scVI_denoised_all_1,p_scVI_denoised_all_2,p_scVI_denoised_all_3,ncol=3,draw=F)
```

```{r fig.height=5,fig.width=18}
p_monocle_scVI_denoised_all
```


```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)*mtx_sizefactor)))
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$scVI_denosied_all=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[8,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```

# Monocle3 using dca+combat

```{r}
#adata=ad$read_h5ad("../final_processed_results/dca Results/adata_all.h5ad")
adata=ad$read_h5ad("../final_processed_results/dca Results New/adata_all.h5ad")
```

```{r}
cell.meta.data=py_to_r(adata$obs)
cell.meta.data$dataset_batch=plyr::mapvalues(cell.meta.data$batch_label,names(maprules),maprules)
gene_ann0=py_to_r(adata$var)
gene_ann=data.frame(gene_short_name = make.unique(rownames(gene_ann0)),
                    row.names = make.unique(rownames(gene_ann0)))
mtx=t(py_to_r(adata$X))
colnames(mtx)=cell.meta.data$cellname
rownames(mtx)=rownames(gene_ann)
mtx_sizefactor=1e4/colSums(mtx)
```

## Using combated latent from dca

```{r}
#mtx=mtx[gene_ann$VarianceType=="HVG",]
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log")

tmp0=py_to_r(adata$obsm["X_dca_latent"]) #original dca
colnames(tmp0)=paste0("PC",1:ncol(tmp0))
reducedDims(cds)$PCA=tmp0

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA")

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden")

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("4"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
#saveRDS(cds,file="cds_dca.rds")
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_dca_latent_all_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_dca_latent_all_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_dca_latent_all_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_dca_latent_all=egg::ggarrange(p_dca_latent_all_1,p_dca_latent_all_2,p_dca_latent_all_3,ncol=3,draw=F)
```

```{r fig.height=5,fig.width=18}
p_monocle_dca_latent_all
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)*mtx_sizefactor)))
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$dca_latent_all=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[9,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```

## HVGs denoised

```{r}
cds <- new_cell_data_set(mtx[rownames(mtx)%in%hvg_genes$genename,], cell_metadata = cell.meta.data,gene_metadata =gene_ann[gene_ann$gene_short_name%in%hvg_genes$genename,,drop=F])
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log",verbose = F)

tmp0=py_to_r(adata$obsm["X_pcahvg"]) #original dca
colnames(tmp0)=paste0("PC",1:ncol(tmp0))
reducedDims(cds)$PCA=tmp0


## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("1","3","4"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_dca_denoised_hvg_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_dca_denoised_hvg_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_dca_denoised_hvg_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_dca_denoised_hvg=egg::ggarrange(p_dca_denoised_hvg_1,p_dca_denoised_hvg_2,p_dca_denoised_hvg_3,ncol=3,draw=F)
```

```{r fig.height=5,fig.width=18}
p_monocle_dca_denoised_hvg
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$dca_denoised_hvg=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[10,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```

## All genes denoised

```{r}
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log",verbose = F)

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

tmp0=py_to_r(adata$obsm["X_pcaall"]) #original dca
colnames(tmp0)=paste0("PC",1:ncol(tmp0))
reducedDims(cds)$PCA=tmp0


## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("2","4","5"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_dca_denoised_all_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_dca_denoised_all_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_dca_denoised_all_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_dca_denoised_all=egg::ggarrange(p_dca_denoised_all_1,p_dca_denoised_all_2,p_dca_denoised_all_3,ncol=3,draw=F)
```

```{r, fig.height=5,fig.width=18}
p_monocle_dca_denoised_all
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,log1p(t(cds_exprs)*mtx_sizefactor)))
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$dca_denoised_all=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[11,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```


# Monocle3 using MNN

```{r}
output=readRDS("../final_processed_results/MNN_corrected_all.rds")
mtx=output@assays$data$corrected
cell.meta.data=colData(output)
cell.meta.data$dataset_batch=plyr::mapvalues(cell.meta.data$batch_label,names(maprules),maprules)
gene_ann=data.frame(gene_short_name = make.unique(rownames(mtx)),
                    row.names = make.unique(rownames(mtx)))

```

### HVGs denoised

```{r}
#mtx=mtx[gene_ann$VarianceType=="HVG",]
cds <- new_cell_data_set(mtx[rownames(mtx)%in%hvg_genes$genename,], 
                         cell_metadata = cell.meta.data,
                         gene_metadata =gene_ann[gene_ann$gene_short_name%in%hvg_genes$genename,,drop=F])
## Step
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log")
#tmp0=py_to_r(adata$obsm["X_dca_latent"])
#colnames(tmp0)=paste0("PC",1:ncol(tmp0))
#reducedDims(cds)$PCA=tmp0

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA")

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden")

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("2","3","4"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
#saveRDS(cds,file="cds_mnn.rds")
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_mnn_denoised_hvg_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_mnn_denoised_hvg_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_mnn_denoised_hvg_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_mnn_denoised_hvg=egg::ggarrange(p_mnn_denoised_hvg_1,p_mnn_denoised_hvg_2,p_mnn_denoised_hvg_3,ncol=3,draw=F)
```

```{r, fig.height=5,fig.width=18}
p_monocle_mnn_denoised_hvg
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)*mtx_sizefactor))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$mnn_denoised_hvg=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[12,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```


### All genes denoised

```{r}
output=readRDS("../final_processed_results/MNN_corrected_all.rds")
mtx=output@assays$data$corrected

cell.meta.data=colData(output)
cell.meta.data$dataset_batch=plyr::mapvalues(cell.meta.data$batch_label,names(maprules),maprules)
gene_ann=data.frame(gene_short_name = make.unique(rownames(mtx)),
                    row.names = make.unique(rownames(mtx)))

#colnames(mtx)=colData(output)$cellname
#rownames(mtx)=rownames(gene_ann)
#mtx_sizefactor=1e4/colSums(mtx)
```

```{r}
#mtx=mtx[gene_ann$VarianceType=="HVG",]
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 32,method="PCA",norm_method="log")
#tmp0=py_to_r(adata$obsm["X_dca_latent"])
#colnames(tmp0)=paste0("PC",1:ncol(tmp0))
#reducedDims(cds)$PCA=tmp0

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA")

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden")

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("6","3","5"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_mnn_denoised_all_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_mnn_denoised_all_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_mnn_denoised_all_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_mnn_denoised_all=egg::ggarrange(p_mnn_denoised_all_1,p_mnn_denoised_all_2,p_mnn_denoised_all_3,ncol=3,draw=F)
```

```{r, fig.height=5,fig.width=18}
p_monocle_mnn_denoised_all
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)*mtx_sizefactor))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$mnn_denoised_all=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[13,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```




# Monocle3 using scanorama


```{r}
adata=ad$read_h5ad("../final_processed_results/scanorama Results/adata_ALL.h5ad")#
```

```{r}
cell.meta.data=py_to_r(adata$obs)
cell.meta.data$dataset_batch=plyr::mapvalues(cell.meta.data$batch_label,names(maprules),maprules)
gene_ann0=py_to_r(adata$raw$var)
gene_ann=data.frame(gene_short_name = make.unique(rownames(gene_ann0)),
                    row.names = make.unique(rownames(gene_ann0)))
mtx=t(py_to_r(adata$X$tocsc()))#adata$raw
colnames(mtx)=cell.meta.data$cellname
rownames(mtx)=rownames(gene_ann)
mtx_sizefactor=1e4/colSums(mtx)
```

## Using latent

```{r}
#mtx=mtx[gene_ann$VarianceType=="HVG",]
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 50,method="PCA",norm_method="log")

tmp0=py_to_r(adata$obsm["X_scanorama"])
colnames(tmp0)=paste0("PC",1:ncol(tmp0))
reducedDims(cds)$PCA=tmp0

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA")

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden")

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("2"))#need to specify
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```

```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
#saveRDS(cds,file="cds_scanorama.rds")
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_monocle_scanorama_latent_hvg_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_monocle_scanorama_latent_hvg_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_monocle_scanorama_latent_hvg_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_scanorama_latent_hvg=egg::ggarrange(p_monocle_scanorama_latent_hvg_1,p_monocle_scanorama_latent_hvg_2,p_monocle_scanorama_latent_hvg_3,ncol=3,draw=F)
```

```{r, fig.height=5,fig.width=18}
p_monocle_scanorama_latent_hvg
```


```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$scanorama_latent_hvg=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```



```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[14,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```


## HVGs denoised

```{r}
#cds <- new_cell_data_set(mtx[rownames(mtx)%in%hvg_genes$genename,], cell_metadata = cell.meta.data,gene_metadata =gene_ann[gene_ann$gene_short_name%in%hvg_genes$genename,,drop=F])
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 50,method="PCA",norm_method="log",verbose = F)

tmp0=py_to_r(adata$obsm["X_hvgpca"])
colnames(tmp0)=paste0("PC",1:ncol(tmp0))
reducedDims(cds)$PCA=tmp0

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("2"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_monocle_scanorama_denoised_hvg_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_monocle_scanorama_denoised_hvg_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_monocle_scanorama_denoised_hvg_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_scanorama_denoised_hvg=egg::ggarrange(p_monocle_scanorama_denoised_hvg_1,p_monocle_scanorama_denoised_hvg_2,p_monocle_scanorama_denoised_hvg_3,ncol=3,draw=F)
```

```{r, fig.height=5,fig.width=18}
p_monocle_scanorama_denoised_hvg
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)*mtx_sizefactor))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$scanorama_denoised_hvg=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[15,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```


## All genes denoised

```{r}
#cds <- new_cell_data_set(mtx[rownames(mtx)%in%hvg_genes$genename,], cell_metadata = cell.meta.data,gene_metadata =gene_ann[gene_ann$gene_short_name%in%hvg_genes$genename,,drop=F])
cds <- new_cell_data_set(mtx, cell_metadata = cell.meta.data,gene_metadata =gene_ann)
## Step 1: Normalize and pre-process the data
cds <- preprocess_cds(cds, num_dim = 30,method="PCA",norm_method="log",verbose = F)

tmp0=py_to_r(adata$obsm["X_pca"])
colnames(tmp0)=paste0("PC",1:ncol(tmp0))
reducedDims(cds)$PCA=tmp0[,1:30]

## Step 2: Remove batch effects with cell alignment
##cds <- align_cds(cds, alignment_group = "BatchID", residual_model_formula_str = NULL)
## Step 3: Reduce the dimensions using UMAP
cds <- reduce_dimension(cds,reduction_method = "UMAP",preprocess_method="PCA",verbose = F)

## Step 4: Cluster the cells
cds <- cluster_cells(cds,reduction_method ="UMAP",cluster_method = "leiden",verbose = F)

# Construct the graph
# Note that, for the rest of the code to run, the graph should be fully (partionly) connected
## Step 5: Learn a graph
cds <- learn_graph(cds, use_partition = T,verbose = F)
colData(cds)$clusters=cds@clusters$UMAP$clusters
p1=plot_cells(cds,color_cells_by = "partition",label_cell_groups = F)+theme(legend.position = "top")
p2=plot_cells(cds,color_cells_by = "clusters",label_cell_groups=F,graph_label_size=2, label_leaves=F,label_branch_points=F)+theme(legend.position = "top")
p=cowplot::plot_grid(p1,p2,align = "h",ncol = 3)
```

```{r, fig.height=5,fig.width=18}
p
```

```{r}
## Step 6: Order cells
# root cells
ids=get_earliest_principal_node(cds,cluster=c("1"))
cds <- order_cells(cds, root_pr_nodes=ids)
#plot_cells(cds,color_cells_by = "pseudotime")
```


```{r,fig.height=5,fig.width=18}
colData(cds)$pseudotime=pseudotime(cds)
colData(cds)$Pseudotime=colData(cds)$pseudotime/max(colData(cds)$pseudotime,na.rm = T)
df_den=pData(cds)[,c("Pseudotime","dataset_batch")]
df_den=as.data.frame(df_den[!is.infinite(df_den$Pseudotime),])
set.seed(10)

theme_use=theme(legend.text = element_text(size=16),
                legend.title = element_text(size=20))

p_monocle_scanorama_denoised_all_1=plot_cells(cds,color_cells_by = "dataset_batch",,graph_label_size=0,alpha=1,cell_size = 0.6)+
  guides(colour = guide_legend(override.aes = list(alpha=0.7, size=5)))+
  theme_use+
  theme(legend.position = "top")
          
p_monocle_scanorama_denoised_all_2=plot_cells(cds,color_cells_by = "Pseudotime",label_branch_points=T,graph_label_size=2,alpha=1,cell_size = 0.6)+
            theme(legend.position = "top",
              legend.title = element_text(vjust = 0.2),
              legend.text = element_text(angle=-50 ),
                  legend.key.height = unit(0.5,"cm"),
                  legend.key.width = unit(1,"cm"))+
            guides(color = guide_colourbar(label.position = "top"))+theme_use

p_monocle_scanorama_denoised_all_3=ggplot(data=df_den)+geom_density(aes(x=Pseudotime,fill=dataset_batch),alpha=0.7)+
            scale_y_continuous(expand = c(0,0))+
            scale_x_continuous(expand = c(0,0))+
            theme(legend.position="top")+theme_use

p_monocle_scanorama_denoised_all=egg::ggarrange(p_monocle_scanorama_denoised_all_1,p_monocle_scanorama_denoised_all_2,p_monocle_scanorama_denoised_all_3,ncol=3,draw=F)
```

```{r, fig.height=5,fig.width=18}
p_monocle_scanorama_denoised_all
```

```{r}
#cds_exprs=FetchData(obj0,vars = c("FCGR3A","S100A8"))
#df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,cds_exprs))
cds_exprs=as.matrix(SingleCellExperiment::counts(cds)[c("FCGR3A","S100A8"),])
df0=data.frame(cbind(pseudotime=pData(cds)$Pseudotime,t(cds_exprs)*mtx_sizefactor))
df0$UMAP_1=reducedDims(cds)$UMAP[,1]
df0$UMAP_2=reducedDims(cds)$UMAP[,2]
df0$BatchID=pData(cds)$dataset_batch
df0=df0[is.finite(df0$pseudotime),]
df0=df0[order(df0$pseudotime,decreasing = F),,drop=F]
df0$x=df0$pseudotime/max(df0$pseudotime)
df_pseudotime_list$scanorama_denoised_all=df0
```

- Feature plots of `FCGR3A` and `S100A8`

```{r}
p=get_plot4(df00 = df0)
```

```{r, fig.height=4,fig.width=18}
p
```

```{r}
df_den=colData(cds)
tt1=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt1
tt2=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T1"],df_den$Pseudotime[df_den$dataset_batch=="T2"])
tt2
tt3=ks.test(df_den$Pseudotime[df_den$dataset_batch=="T2"],df_den$Pseudotime[df_den$dataset_batch=="T3"])
tt3
```

```{r}
Stable5[16,2:4]=matrix(get_p_new(c(tt1$p.value,tt2$p.value,tt3$p.value),c(tt1$statistic,tt2$statistic,tt3$statistic)),1,3)
```

```{r}
suppressPackageStartupMessages(library(cowplot))
```

```{r}
#Supplementary Table 5
write.table(Stable5,file="KS_table.csv",sep=",")
openxlsx::write.xlsx(Stable5,file="KS_table.xlsx")

Stable5 %>%
  kable() %>%
  kable_styling()
```

# Figures

##  Main Figure (Figure 5)

```{r, fig.width=30,fig.height=25}
fig_width=30
fig_height=25
labels=letters[1:5]
Methods=c("Raw count HVGs","Raw count All","CarDEC (latent)","CarDEC (denoised HVGs)","CarDEC (denoised All)","scVI (latent)","scVI (denoised HVGs)","scVI (denoised All)","DCA (latent)","DCA (denoised HVGs)","DCA (denoised All)","MNN (denoised HVGs)","MNN (denoised All)","Scanorama (latent)","Scanorama (denoised HVGs)","Scanorama (denoised All)")
use_id=c(4,15,10,7,12)

get_draw_plot=function(plot_id=1,plist0){
  x=0.02
  y=1-plot_id/5
  width=0.98
  height=1/5-0.01 # total number of figures is 12
  pp=draw_plot(egg::ggarrange(plots=plist0,nrow = 1,draw = F,newpage = F),x = x,y = y,width = width,height = height)
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(pp)
}

get_label_pos=function(plot_id=1){
  x=0
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/5-1/30))
}

get_title_pos=function(plot_id=1){
  x=0.015
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/10-1/60))
}

get_plot_list=function(x,y){
  x0=rep(list(),length=length(x)+length(y))
  x0[1:length(x)]=x[1:3]
  x0[(length(x)+1):length(x0)]=y[1:2]
  return(x0)
}
p=ggdraw()+get_draw_plot(1,get_plot_list(list(p_carDEC_denoised_hvg_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                              p_carDEC_denoised_hvg_2,
                                              p_carDEC_denoised_hvg_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")),
                                         get_plot4_sep(df_pseudotime_list[[4]])[c(3,4)]))+
get_draw_plot(2,get_plot_list(list(p_monocle_scanorama_denoised_hvg_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                   p_monocle_scanorama_denoised_hvg_2,
                                   p_monocle_scanorama_denoised_hvg_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")),
                              get_plot4_sep(df_pseudotime_list[[15]])[c(3,4)]))+
  get_draw_plot(3,get_plot_list(list(p_dca_denoised_hvg_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                     p_dca_denoised_hvg_2,
                                     p_dca_denoised_hvg_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")), 
                                get_plot4_sep(df_pseudotime_list[[10]])[c(3,4)]))+
  get_draw_plot(4,get_plot_list(list(p_scvi_denoised_hvg_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                     p_scvi_denoised_hvg_2,
                                     p_scvi_denoised_hvg_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")),
                                get_plot4_sep(df_pseudotime_list[[7]])[c(3,4)]))+
get_draw_plot(5,get_plot_list(list(p_mnn_denoised_hvg_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                   p_mnn_denoised_hvg_2,
                                   p_mnn_denoised_hvg_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")),
                              get_plot4_sep(df_pseudotime_list[[12]])[c(3,4)]))

for (i in 1:length(use_id)){
  p=p+draw_label(labels[i],x=get_label_pos(i)[1],y=get_label_pos(i)[2],size=30,color="black",hjust = 0,vjust = 1)+
    draw_label(paste0(Methods[use_id[i]],collapse = ""),
               x=get_title_pos(i)[1],
               y=get_title_pos(i)[2],
               size=20,
               color = "black",angle = 90,hjust = 0.5,vjust = 0.5)
}
```

```{r, fig.width=30,fig.height=25}
p
```

```{r}
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Main_fig1.pdf",width = 30,height = 25)
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Main_fig1.tiff",width = 30,height = 25,compression="lzw")
```


##  Supplementary Figures about monocyte

1. Raw 

```{r}
plist0=rep(list(),length=9)
plist0[1:3]=list(p_ori_all_1+scale_color_brewer(name="dataset_batch",palette = "Set2")+theme_use,
                 p_ori_all_2,
                 p_ori_all_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")+theme(plot.margin = unit(c(0,0,1,0),"cm")))[1:3]
plist0[4:5]=get_plot4_sep(df_pseudotime_list[[2]])[1:2]
plist0[[6]]=ggplot()+theme_void()+theme(plot.margin = unit(c(0,0,1,0),"cm"))
plist0[7:8]=get_plot4_sep(df_pseudotime_list[[2]])[3:4]
plist0[[9]]=ggplot()+theme_void()
p=ggdraw()+draw_plot(egg::ggarrange(plots = plist0,ncol=3,draw = F),x=0,y=0,width=1,height=1)+
 draw_label("a",x=0,y=1-0.02,size=30,color="black",hjust = 0,vjust = 1)+
 draw_label("b",x=0,y=2/3-0.03,size=30,color="black",hjust = 0,vjust = 1)+
   draw_label("c",x=0,y=1/3-0.03,size=30,color="black",hjust = 0,vjust = 1)
```

```{r, fig.width=18,fig.height=16}
p
```

```{r}
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Supp_fig1.pdf",width = 18,height = 16)
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Supp_fig1.tiff",width = 18,height = 16,compression="lzw")
```

2. All genes

```{r,fig.width=36,fig.height=30}
fig_width=30
fig_height=25
labels=letters[1:5]
Methods=c("Raw count HVGs","Raw count All","CarDEC (latent)","CarDEC (denoised HVGs)","CarDEC (denoised All)","scVI (latent)","scVI (denoised HVGs)","scVI (denoised All)","DCA (latent)","DCA (denoised HVGs)","DCA (denoised All)","MNN (denoised HVGs)","MNN (denoised All)","Scanorama (latent)","Scanorama (denoised HVGs)","Scanorama (denoised All)")
use_id=c(5,16,11,8,13)

get_draw_plot=function(plot_id=1,plist0){
  x=0.02
  y=1-plot_id/5
  width=0.98
  height=1/5-0.01 # total number of figures is 12
  pp=draw_plot(egg::ggarrange(plots=plist0,nrow = 1,draw = F,newpage = F),x = x,y = y,width = width,height = height)
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(pp)
}

get_label_pos=function(plot_id=1){
  x=0
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/5-1/30))
}

get_title_pos=function(plot_id=1){
  x=0.015
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/10-1/60))
}
get_plot_list=function(x,y){
  x0=rep(list(),length=length(x)+length(y))
  x0[1:length(x)]=x[1:3]
  x0[(length(x)+1):length(x0)]=y[1:2]
  return(x0)
}

p=ggdraw()+
  get_draw_plot(1,get_plot_list(list(p_carDEC_denoised_all_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                     p_carDEC_denoised_all_2,
                                     p_carDEC_denoised_all_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")), 
                                get_plot4_sep(df_pseudotime_list[[5]])[c(3,4)]))+
 get_draw_plot(2,get_plot_list(list(p_monocle_scanorama_denoised_all_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                   p_monocle_scanorama_denoised_all_2,
                                   p_monocle_scanorama_denoised_all_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")),
                              get_plot4_sep(df_pseudotime_list[[16]])[c(3,4)]))+
  get_draw_plot(3,get_plot_list(list(p_dca_denoised_all_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                     p_dca_denoised_all_2,
                                     p_dca_denoised_all_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")), 
                                get_plot4_sep(df_pseudotime_list[[11]])[c(3,4)]))+
  get_draw_plot(4,get_plot_list(list(p_scVI_denoised_all_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                     p_scVI_denoised_all_2,
                                     p_scVI_denoised_all_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")), 
                                get_plot4_sep(df_pseudotime_list[[8]])[c(3,4)]))+
  get_draw_plot(5,get_plot_list(list(p_mnn_denoised_all_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                     p_mnn_denoised_all_2,
                                     p_mnn_denoised_all_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")), 
                                get_plot4_sep(df_pseudotime_list[[13]])[c(3,4)]))
 
for (i in 1:length(use_id)){
  p=p+draw_label(labels[i],x=get_label_pos(i)[1],y=get_label_pos(i)[2],size=30,color="black",hjust = 0,vjust = 1)+
    draw_label(paste0(Methods[use_id[i]],collapse = ""),
               x=get_title_pos(i)[1],
               y=get_title_pos(i)[2],
               size=20,
               color = "black",angle = 90,hjust = 0.5,vjust = 0.5)
}
```

```{r, fig.width=30,fig.height=25}
p
```

```{r}
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Supp_fig2.pdf",width = 30,height = 25,limitsize = F)
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Supp_fig2.tiff",width = 30,height = 25,limitsize = F,compression="lzw")
```


3. Latent

```{r,fig.width=36,fig.height=30}
fig_width=30
fig_height=25
labels=letters[1:5]
Methods=c("Raw count HVGs","Raw count All","CarDEC (latent)","CarDEC (denoised HVGs)","CarDEC (denoised All)","scVI (latent)","scVI (denoised HVGs)","scVI (denoised All)","DCA (latent)","DCA (denoised HVGs)","DCA (denoised All)","MNN (denoised HVGs)","MNN (denoised All)","Scanorama (latent)","Scanorama (denoised HVGs)","Scanorama (denoised All)")

use_id=c(3,14,9,6,1)

get_draw_plot=function(plot_id=1,plist0){
  x=0.02
  y=1-plot_id/5
  width=0.98
  height=1/5-0.01 # total number of figures is 12
  pp=draw_plot(egg::ggarrange(plots=plist0,nrow = 1,draw = F,newpage = F),x = x,y = y,width = width,height = height)
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(pp)
}

get_label_pos=function(plot_id=1){
  x=0
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/5-1/30))
}

get_title_pos=function(plot_id=1){
  x=0.015
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/10-1/60))
}
get_plot_list=function(x,y){
  x0=rep(list(),length=length(x)+length(y))
  x0[1:length(x)]=x[1:3]
  x0[(length(x)+1):length(x0)]=y[1:2]
  return(x0)
}


p=ggdraw()+get_draw_plot(1,get_plot_list(list(p_carDEC_latent_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                              p_carDEC_latent_2,
                                              p_carDEC_latent_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")),
                                         get_plot4_sep(df_pseudotime_list[[3]])[c(3,4)]))+
  get_draw_plot(2,get_plot_list(list(p_monocle_scanorama_latent_hvg_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                   p_monocle_scanorama_latent_hvg_2,
                                   p_monocle_scanorama_latent_hvg_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")),
                              get_plot4_sep(df_pseudotime_list[[14]])[c(3,4)]))+
  get_draw_plot(3,get_plot_list(list(p_dca_latent_all_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                     p_dca_latent_all_2,
                                     p_dca_latent_all_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")), 
                                get_plot4_sep(df_pseudotime_list[[9]])[c(3,4)]))+
  get_draw_plot(4,get_plot_list(list(p_scVI_latent_all_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                     p_scVI_latent_all_2,
                                     p_scVI_latent_all_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")), 
                                get_plot4_sep(df_pseudotime_list[[6]])[c(3,4)]))+
  get_draw_plot(5,get_plot_list(list(p_ori_1+scale_color_brewer(name="dataset_batch",palette = "Set2"),
                                     p_ori_2,
                                     p_ori_3+scale_fill_brewer(name="dataset_batch",palette = "Set2")), 
                                get_plot4_sep(df_pseudotime_list[[1]])[c(3,4)]))
for (i in 1:length(use_id)){
  p=p+draw_label(labels[i],x=get_label_pos(i)[1],y=get_label_pos(i)[2],size=30,color="black",hjust = 0,vjust = 1)+
    draw_label(paste0(Methods[use_id[i]],collapse = ""),
               x=get_title_pos(i)[1],
               y=get_title_pos(i)[2],
               size=20,
               color = "black",angle = 90,hjust = 0.5,vjust = 0.5)
}
```

```{r, fig.width=30,fig.height=25}
p
```

```{r}
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Supp_fig3.pdf",width = 30,height = 25,limitsize = F)
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Supp_fig3.tiff",width = 30,height = 25,limitsize = F,compression="lzw")
```

## S100A8 and FCGR3A's feature plots

1. monocles' UMAP of denoised counts from HVGs 

```{r}
fig_width=15
fig_height=25
labels=letters[1:13]
Methods=c("Raw count HVGs","Raw count All","CarDEC (latent)","CarDEC (denoised HVGs)","CarDEC (denoised All)","scVI (latent)","scVI (denoised HVGs)","scVI (denoised All)","DCA (latent)","DCA (denoised HVGs)","DCA (denoised All)","MNN (denoised HVGs)","MNN (denoised All)","Scanorama (latent)","Scanorama (denoised HVGs)","Scanorama (denoised All)")
use_id=c(4,15,10,7,12)
#use_id=c(3,9,6,1)
get_draw_plot=function(plot_id=1,plist0){
  x=0.02
  y=1-plot_id/5
  width=0.98
  height=1/5-0.01 # total number of figures is 12
  pp=draw_plot(egg::ggarrange(plots=plist0,nrow = 1,draw = F,newpage = F),x = x,y = y,width = width,height = height)
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(pp)
}

get_label_pos=function(plot_id=1){
  x=0
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/5-1/100))
}

get_title_pos=function(plot_id=1){
  x=0.015
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/10-1/100))
}
get_plot_list=function(x,y){
  x0=rep(list(),length=length(x)+length(y))
  x0[1:length(x)]=x[1:3]
  x0[(length(x)+1):length(x0)]=y[1:2]
  return(x0)
}

p=ggdraw()
for(i in 1:length(use_id)){
  p=p+get_draw_plot(i,get_plot4_sep(df_pseudotime_list[[use_id[i]]])[c(1,2)])
}

for (i in 1:length(use_id)){
  p=p+draw_label(labels[i],x=get_label_pos(i)[1],y=get_label_pos(i)[2],size=30,color="black",hjust = 0,vjust = 1)+
    draw_label(paste0(Methods[use_id[i]],collapse = ""),
               x=get_title_pos(i)[1],
               y=get_title_pos(i)[2],
               size=20,
               color = "black",angle = 90,hjust = 0.5,vjust = 0.5)
}

```

```{r fig.width=15,fig.height=25}
ggsave("./revised_figures/CarDEC_monocyte_Supp_fig4.pdf",p,width = 15,height = 25)
ggsave("./revised_figures/CarDEC_monocyte_Supp_fig4.tiff",p,width = 15,height = 25,compression="lzw")
p
```

2. monocles' UMAP of denoised counts from All genes 


```{r}
fig_width=15
fig_height=25
labels=letters[1:13]
Methods=c("Raw count HVGs","Raw count All","CarDEC (latent)","CarDEC (denoised HVGs)","CarDEC (denoised All)","scVI (latent)","scVI (denoised HVGs)","scVI (denoised All)","DCA (latent)","DCA (denoised HVGs)","DCA (denoised All)","MNN (denoised HVGs)","MNN (denoised All)","Scanorama (latent)","Scanorama (denoised HVGs)","Scanorama (denoised All)")
use_id=c(5,16,11,8,13)
get_draw_plot=function(plot_id=1,plist0){
  x=0.02
  y=1-plot_id/5
  width=0.98
  height=1/5-0.01 # total number of figures is 12
  pp=draw_plot(egg::ggarrange(plots=plist0,nrow = 1,draw = F,newpage = F),x = x,y = y,width = width,height = height)
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(pp)
}

get_label_pos=function(plot_id=1){
  x=0
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/5-1/100))
}

get_title_pos=function(plot_id=1){
  x=0.015
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/10-1/100))
}
get_plot_list=function(x,y){
  x0=rep(list(),length=length(x)+length(y))
  x0[1:length(x)]=x[1:3]
  x0[(length(x)+1):length(x0)]=y[1:2]
  return(x0)
}

p=ggdraw()
for(i in 1:length(use_id)){
  p=p+get_draw_plot(i,get_plot4_sep(df_pseudotime_list[[use_id[i]]])[c(1,2)])
}

for (i in 1:length(use_id)){
  p=p+draw_label(labels[i],x=get_label_pos(i)[1],y=get_label_pos(i)[2],size=30,color="black",hjust = 0,vjust = 1)+
    draw_label(paste0(Methods[use_id[i]],collapse = ""),
               x=get_title_pos(i)[1],
               y=get_title_pos(i)[2],
               size=20,
               color = "black",angle = 90,hjust = 0.5,vjust = 0.5)
}

```

```{r fig.width=15,fig.height=25}
ggsave("./revised_figures/CarDEC_monocyte_Supp_fig5.pdf",p,width = 15,height = 25)
ggsave("./revised_figures/CarDEC_monocyte_Supp_fig5.tiff",p,width = 15,height = 25,compression="lzw")
p
```


3. monocles' UMAP based on different methods' latent


```{r}
fig_width=15
fig_height=25
labels=letters[1:13]
Methods=c("Raw count HVGs","Raw count All","CarDEC (latent)","CarDEC (denoised HVGs)","CarDEC (denoised All)","scVI (latent)","scVI (denoised HVGs)","scVI (denoised All)","DCA (latent)","DCA (denoised HVGs)","DCA (denoised All)","MNN (denoised HVGs)","MNN (denoised All)","Scanorama (latent)","Scanorama (denoised HVGs)","Scanorama (denoised All)")
use_id=c(3,14,9,6,1)
get_draw_plot=function(plot_id=1,plist0){
  x=0.02
  y=1-plot_id/5
  width=0.98
  height=1/5-0.01 # total number of figures is 12
  pp=draw_plot(egg::ggarrange(plots=plist0,nrow = 1,draw = F,newpage = F),x = x,y = y,width = width,height = height)
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(pp)
}

get_label_pos=function(plot_id=1){
  x=0
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/5-1/100))
}

get_title_pos=function(plot_id=1){
  x=0.015
  y=1-plot_id/5
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/10-1/100))
}
get_plot_list=function(x,y){
  x0=rep(list(),length=length(x)+length(y))
  x0[1:length(x)]=x[1:3]
  x0[(length(x)+1):length(x0)]=y[1:2]
  return(x0)
}

p=ggdraw()
for(i in 1:length(use_id)){
  p=p+get_draw_plot(i,get_plot4_sep(df_pseudotime_list[[use_id[i]]])[c(1,2)])
}

for (i in 1:length(use_id)){
  p=p+draw_label(labels[i],x=get_label_pos(i)[1],y=get_label_pos(i)[2],size=30,color="black",hjust = 0,vjust = 1)+
    draw_label(paste0(Methods[use_id[i]],collapse = ""),
               x=get_title_pos(i)[1],
               y=get_title_pos(i)[2],
               size=20,
               color = "black",angle = 90,hjust = 0.5,vjust = 0.5)
}

```

```{r fig.width=15,fig.height=25}
ggsave("./revised_figures/CarDEC_monocyte_Supp_fig6.pdf",p,width = 15,height = 25)
ggsave("./revised_figures/CarDEC_monocyte_Supp_fig6.tiff",p,width = 15,height = 25,compression="lzw")
p
```

4. Combined above three figures

```{r}
fig_width=24
fig_height=36#
labels=letters[1:16]
Methods=c("Raw count HVGs","Raw count All","CarDEC (latent)","CarDEC (denoised HVGs)","CarDEC (denoised All)","scVI (latent)","scVI (denoised HVGs)","scVI (denoised All)","DCA (latent)","DCA (denoised HVGs)","DCA (denoised All)","MNN (denoised HVGs)","MNN (denoised All)","Scanorama (latent)","Scanorama (denoised HVGs)","Scanorama (denoised All)")
use_id=c(1,3:5,14:16,6:13)
get_draw_plot=function(plot_id=1,plist0){
  x=ifelse(plot_id%%2==1,0,0.5075)
  y=1-floor((plot_id+1)/2)/8
  width=1/2-0.015
  height=1/8-0.01 # total number of figures is 12
  pp=draw_plot(egg::ggarrange(plots=plist0,ncol = 2,draw = F,newpage = F),x = x,y = y,width = width,height = height)
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(pp)
}

get_label_pos=function(plot_id=1){
  x=ifelse(plot_id%%2==1,0,0.5075)
  y=1-floor((plot_id+1)/2)/8
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/8-1/100))
}

get_title_pos=function(plot_id=1){
  x=ifelse(plot_id%%2==1,0,0.5075)+0.5/2
  y=1-floor((plot_id+1)/2)/8
  #draw_label(labels[plot_id],x=x,y=y+1/6,hjust =1,vjust = 0.5,size = 35)
  return(c(x,y+1/8-1/300))
}

p=ggdraw()
for(i in 1:length(use_id)){
  p=p+get_draw_plot(i,get_plot4_sep(df_pseudotime_list[[use_id[i]]])[c(1,2)])
}

for (i in 1:length(use_id)){
  p=p+draw_label(labels[i],x=get_label_pos(i)[1],y=get_label_pos(i)[2],size=18,color="black",hjust = 0,vjust = 1)+
    draw_label(paste0(Methods[use_id[i]],collapse = ""),x=get_title_pos(i)[1],y=get_title_pos(i)[2],size=25,color = "black",vjust = 1)
}

```

```{r fig.width=24,fig.height=36}
p
```

```{r}
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Supp_fig456_combined.pdf",width = 24,height = 36)
ggsave(p,filename = "./revised_figures/CarDEC_monocyte_Supp_fig456_combined.tiff",width = 24,height = 36,compression="lzw")
```


```{r}
#save object 
rm(mtx)
rm(output)
rm(adata)
rm(cds)
rm(obj0)
rm(raw.data)
gc()
save.image(file="carDEC_monocyte_final_revised.RData")#include scanorama,carDEC_monocyte_final.RData not include scanorama
#load("carDEC_monocyte_final_revised.RData")
```

```{r}
#load("carDEC_monocyte_final_revised.RData")
```

```{r}
sessionInfo()
```
