
require(gdata)
library(readxl)
library(dplyr)
library(patchwork)
library(presto)
library(uwot)
library(harmony)
library(magrittr)
library(gridExtra)
library(reticulate)
library(ggbeeswarm)
library(ggrepel)
library(plyr)
library(ggplot2)
library(gridExtra)
library(RColorBrewer)
library(MASS)
library(matrixStats)
library(viridis)
library(grid)
library(singlecellmethods)
library(umap)
library(tidyverse)
library(ggpubr)
library(rstatix)
library(Seurat)
library(CCA)
library(tidyverse)
library(ggthemes)
library(symphony)
library(pheatmap)
library(ggrastr)


meta_colors <- list(
    
    "cell_type" = c(
        "B cell/plasma cell" = "#FDBF6F",
        "B cell" = "#FDBF6F",
        "DC" = "#CAB2D6",
        "Fibroblast" = "#08519C",
        "Stromal cell" = "#08519C",
        "Endothelial cell" = "#A6CEE3",
        "Macrophage" = "#6A3D9A",
        "Monocyte" = "#6A3D9A",
        "Myeloid cell" = "#6A3D9A",
        "T cell" = "#B2DF8A",
        "NK" = "#33A02C"
    ),

    "class" = c(
        "M" = "#9E0142",
        "T + M" = "#F46D43",
        "T + B" = "#FEE08B",
        "T + F" = "#E6F598",
        "F" = "#ABDDA4",
        "E + F + M" = "#66C2A5",
        "OA" = "grey"
    )
)


#' Basic function to plot cells, colored and faceted by metadata variables
#' 
#' @param metadata metadata, with UMAP labels in UMAP1 and UMAP2 slots
#' @param title Plot title
#' @param color.by metadata column name for phenotype labels
#' @param facet.by metadata column name for faceting
#' @param color.mapping custom color mapping
#' @param show.legend Show cell type legend

plotBasic = function(umap_labels,                # metadata, with UMAP labels in UMAP1 and UMAP2 slots
                        title = 'Query',         # Plot title
                        color.by = 'cell_type',  # metadata column name for coloring
                        facet.by = NULL,         # (optional) metadata column name for faceting
                        color.mapping = NULL,    # custom color mapping
                        legend.position = 'right') {  # Show cell type legend
    
    p = umap_labels %>%
            dplyr::sample_frac(1L) %>% # permute rows randomly
            ggplot(aes(x = UMAP1, y = UMAP2)) + 
            geom_point_rast(aes(col = get(color.by)), size = 0.3, stroke = 0.2, shape = 16)
        if (!is.null(color.mapping)) { p = p + scale_color_manual(values = color.mapping) }
    
    # Default formatting
    p = p + theme_bw() +
            labs(title = title, color = color.by) + 
            theme(plot.title = element_text(hjust = 0.5)) +
            theme(legend.position=legend.position) +
            theme(legend.text = element_text(size=8), legend.title=element_text(size=12)) + 
            guides(colour = guide_legend(override.aes = list(size = 4))) + guides(alpha = FALSE)

    if(!is.null(facet.by)) {
        p = p + facet_wrap(~get(facet.by)) +
                theme(strip.text.x = element_text(size = 12)) }    
    return(p)
}


mRNA_exprs <- readRDS("/data/srlab/fzhang/amp/results/2020_01_22_AMP_cite_seq_QC/mRNA_exprs_2020-01-24.rds")
meta_all <- readRDS("/data/srlab1/public/srcollab/AMP_Phase_2/singlecell_result/fine_cluster_all_celltypes_82samples_2021-04-23.rds")
mRNA_exprs <- mRNA_exprs[, meta_all$cell]
dim(mRNA_exprs)

exprs_norm <- normalizeData(mRNA_exprs, method = "log")
dim(exprs_norm)
adt_exprs_norm <- readRDS("/data/srlab/fzhang/amp/results/2020_01_22_AMP_cite_seq_QC/adt_exprs_norm_filter_2020-02-18.rds") # adt_exprs_norm are CLR normalized and QCed
adt_exprs_norm <- adt_exprs_norm[, meta_all$cell]
dim(adt_exprs_norm)


meta_all <- meta_all[, c("cell", "sample", "cluster_number", "cluster_name", "cell_type")]


# Feature selection using vst on raw counts
meta_all$sample <- as.character(meta_all$sample)

# Select variable genes per sample and then take the union
var_genes <- vargenes_vst(mRNA_exprs, topn = 1000, groups = meta_all$sample)
length(var_genes)


# Check if VST exclude cell cycle and MT and ribosome genes
mt_ri <- grep("^MT-|^RPL|^RPS|MALAT1|MIR-", row.names(exprs_norm), value = TRUE)
cycle_prolif <- c(Seurat::cc.genes$s.genes, Seurat::cc.genes$g2m.genes)
genes_exclude <- union(mt_ri, cycle_prolif)
var_genes <- var_genes[-which(var_genes %in% genes_exclude)]
length(var_genes)
rm(mRNA_exprs) # Don't need raw count data anymore


# Normalized data with highly variable genes
ref_exp = exprs_norm[var_genes, ]
dim(ref_exp)
rm(exprs_norm)


# Calculate and save the mean and standard deviations for each gene
vargenes_means_sds = tibble(symbol = var_genes, mean = Matrix::rowMeans(ref_exp))
vargenes_means_sds$stddev = singlecellmethods::rowSDs(ref_exp, vargenes_means_sds$mean)
head(vargenes_means_sds)
dim(vargenes_means_sds)
saveRDS(vargenes_means_sds, "all_cells_vargenes_means_sds_2021-05-20.rds")


# Use the most cell type-specific proteins based on KL divergence for each protein
adt_exprs_norm <- adt_exprs_norm[c("CD11b_prot", "CD64_prot", "CD34_prot", "Podoplanin_prot", "CD3_prot", "CD11c_prot", "CD31_prot", "CD206/MMR_prot",
                                   "FR-beta_prot", "CD278/ICOS_prot", "CD90/THY1_prot", "HLA-DR_prot", "CD45RO_prot", "EGFR_prot",
                                   "CD141/Thrombomodulin_prot", "CD27(LG.3A10)_prot", "CD112/Nectin-2_prot", "CD4_prot", "CD45RA_prot",
                                   "CD86_prot", "CD146_prot", "CD45(2D1)_prot", "CD14_prot", "CD127/IL-7R-alpha_prot", "CD19_prot", "CD279/PD-1_prot",
                                   "CD144/VE-Cadherin_prot", "CD107a/LAMP1_prot", "CD69_prot", "CD44_prot", "CD155/PVR_prot", "CD8a_prot", "CD20_prot",
                                   "CD21_prot", "CD68_prot", "CD163_prot", "CD192/CCR2_prot", "CD1c_prot", "CD56/NCAM_prot"),]
dim(adt_exprs_norm)

# adt_exprs_norm <- adt_exprs_norm[c("CD34_prot", "Podoplanin_prot", "CD3_prot", "CD11c_prot",
#                                    "FR-beta_prot", "CD278/ICOS_prot", "CD90/THY1_prot", "CD45RO_prot","CD27(LG.3A10)_prot", "CD4_prot", "CD45RA_prot",
#                                    "CD86_prot", "CD146_prot", "CD45(2D1)_prot", "CD14_prot", "CD127/IL-7R-alpha_prot", "CD19_prot", "CD8a_prot", "CD20_prot",
#                                    "CD21_prot", "CD68_prot", "CD163_prot", "CD192/CCR2_prot", "CD1c_prot", "CD56/NCAM_prot"),]
# dim(adt_exprs_norm)


mRNA_scale <- singlecellmethods::scaleData(ref_exp) 
# Exactly the same results with singlecellmethods::scaleDataWithStats(ref_exp, vargenes_means_sds$mean, vargenes_means_sds$stddev, 1)
protein_scale <- singlecellmethods::scaleData(adt_exprs_norm)

rm(ref_exp)
rm(adt_exprs_norm)
dim(mRNA_scale)
dim(protein_scale)


gc()


# This step can be PCA or CCA
# Run CCA
set.seed(0)
system.time({
    res_cca = cc(t(mRNA_scale), t(protein_scale))
    })

# Save the results
saveRDS(res_cca, "res_cca_all_328201cells_VST5751genes_39proteins_2021-05-20.rds")

rm(mRNA_scale)
rm(protein_scale)


## Extracting loadings
loadings <- res_cca$xcoef[,1:20]
cca_res <- res_cca$scores$xscores[,1:20]
dim(loadings)
dim(cca_res)


## Build reference
set.seed(0)
# We have tested parameter theta (0,1,2) and select theta = 1 
system.time({
    ref_harmObj <- HarmonyMatrix(cca_res, meta_all, c("sample"), 
                                 theta = c(1), plot_convergence = TRUE, 
                                 epsilon.cluster = -Inf, epsilon.harmony = -Inf, 
                                 max.iter.harmony = 10, max.iter.cluster = 30, 
                                 do_pca = F, verbose = T, return_object = TRUE
                                 )
})
rm(cca_res)


# vargenes_means_sds <- readRDS("all_cells_vargenes_means_sds_2021-05-13.rds")


# To run the next function buildReferenceFromHarmonyObj(), you need to input the saved gene loadings (loadings) and vargenes_means_sds.
reference = symphony::buildReferenceFromHarmonyObj(
                           ref_harmObj,            # output object from HarmonyMatrix()
                           meta_all,           # reference cell metadata
                           vargenes_means_sds,     # gene names, means, and std devs for scaling
                           loadings,               # genes x PCs matrix
                           verbose = TRUE,         # verbose output
                           do_umap = TRUE,         # Set to TRUE only when UMAP model was saved for reference
                           save_uwot_path = './all_cells_uwot_model_2021-05-20')


# Save all cell reference
saveRDS(reference, './all_cells_reference_2021-05-20.rds') 

# meta_data: metadata
# vargenes: variable genes, means, and standard deviations used for scaling
# loadings: gene loadings for projection into pre-Harmony PC space
# R: Soft cluster assignments
# Z_orig: Pre-Harmony PC embedding
# Z_corr: Harmonized PC embedding
# centroids: locations of final Harmony soft cluster centroids
# cache: pre-calculated reference-dependent portions of the mixture model
# umap: UMAP coordinates
# save_uwot_path: path to saved uwot model (for query UMAP projection into reference UMAP coordinates)


class(reference)


reference <- readRDS("all_cells_reference_2021-05-20.rds")


# res_cca_all_328201cells_VST5751genes_39proteins_2021-05-20.rds
umap_labels = cbind(reference$meta_data, reference$umap$embedding)

# fig.size(3, 5)
options(repr.plot.height = 4, repr.plot.width = 5)
ggplot() +
  geom_point(
    data = umap_labels[sample(nrow(umap_labels)),],
    mapping = aes_string(x = "UMAP1", y = "UMAP2", fill = "cell_type"),
    size = 0.1, stroke = 0.0001, shape = 21
  ) +
#   facet_wrap( ~ new_cluster_name) +
  scale_fill_manual(values = meta_colors$cell_type, name = "") + 
  theme_bw(base_size = 10) +
  theme(
    legend.position = "none",
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(color="black", size=10)
  ) 
ggsave(file = paste("all_celltype_umap", ".png", sep = ""), width = 5, height = 4, dpi = 300)


# adt_exprs <- readRDS("/data/srlab/fzhang/amp/results/2020_01_22_AMP_cite_seq_QC/adt_exprs_2020-01-24.rds")
# adt_exprs <- adt_exprs[, meta_all$cell]
# adt_exprs_norm <- adt_exprs %>% singlecellmethods::normalizeData(method = "cellCLR")

adt_exprs_norm <- readRDS("/data/srlab/fzhang/amp/results/2020_01_22_AMP_cite_seq_QC/adt_exprs_norm_filter_2020-02-18.rds") # adt_exprs_norm are CLR normalized and QCed
adt_exprs_norm <- adt_exprs_norm[, umap_labels$cell]
dim(adt_exprs_norm)


# rownames(adt_exprs_norm) <- substr(rownames(adt_exprs_norm), 1, nchar(rownames(adt_exprs_norm) )-5)
plot_protein <- c( "CD4_prot", "CD8a_prot", "CD127/IL-7R-alpha_prot", "CD56/NCAM_prot", 
                   "CD14_prot", "CD163_prot", "CD192/CCR2_prot", "CD1c_prot", 
                   "Podoplanin_prot", "CD90/THY1_prot", "CD34_prot", "CD31_prot", 
                   "CD146_prot", "CD19_prot", "CD20_prot"
                 )

rownames(adt_exprs_norm) <- substr(rownames(adt_exprs_norm), 1, nchar(rownames(adt_exprs_norm) )-5)
plot_protein <- substr(plot_protein, 1, nchar(plot_protein)-5)


myplots <- list()
for (i in 1:length(plot_protein)) {
    protein <- plot_protein[i]
    
    max.cutoff = quantile(adt_exprs_norm[protein,], .99)
    min.cutoff = quantile(adt_exprs_norm[protein,], .01)

    tmp <- sapply(X = adt_exprs_norm[protein,], FUN = function(x) {
        return(ifelse(test = x > max.cutoff, yes = max.cutoff,
            no = x))
    })
    tmp <- sapply(X = tmp, FUN = function(x) {
        return(ifelse(test = x < min.cutoff, yes = min.cutoff,
            no = x))
    })
  
    umap_labels$protein <- as.numeric(tmp)
    

    
   ind <- paste("p", i, sep = "")
   ind <- ggplot(
                 data = umap_labels[sample(nrow(umap_labels)),],
                 aes(x = UMAP1, y = UMAP2)) + 
      geom_point(mapping = aes(color = protein), size = 0.000001) + 
      scale_color_viridis(option = "plasma", end = .9) +
      labs(x="", y="")+
      theme_bw(base_size = 15)+
      theme(
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        panel.grid = element_blank(),
        plot.title = element_text(color="black", size=35), # face="bold.italic"
        legend.position = "none") +
      labs(title = protein)
   myplots[[i]] <- ind
    
    
}

options(repr.plot.height = 12, repr.plot.width = 22)
p <- do.call("grid.arrange", c(myplots, ncol = 5))
ggsave(file = paste("protein_umap", ".png", sep = ""), p, width = 22, height = 9, dpi = 300)


# Fine grained subpopulations are generally maintained in the global transcriptomic space.
# To more accurately define subpopulations, we recommend our cell type-specific references.

colors37 = c("#466791","#60bf37","#953ada","#4fbe6c","#ce49d3","#a7b43d","#5a51dc","#d49f36","#552095","#507f2d","#db37aa","#84b67c","#a06fda","#df462a","#5b83db","#c76c2d","#4f49a3","#82702d","#dd6bbb","#334c22","#d83979","#55baad","#dc4555","#62aad3","#8c3025","#417d61","#862977","#bba672","#403367","#da8a6d","#a79cd4","#71482c","#c689d0","#6b2940","#d593a7","#895c8b","#bd5975")
cluster_center <- umap_labels %>%
                  group_by(cluster_name) %>%
                  summarise_at(vars(UMAP1, UMAP2), funs(median(., na.rm=TRUE)))
cluster_center <- as.data.frame(cluster_center)
cluster_center$cluster_name <- as.character(cluster_center$cluster_name)

options(repr.plot.height = 8, repr.plot.width = 10)
ggplot(umap_labels[sample(nrow(umap_labels)),],
       aes(x = UMAP1, y = UMAP2, fill= cluster_name)
      ) +
  geom_point(size = 0.1, stroke = 0.0001, shape = 21, alpha = 0.9) +
  geom_label_repel(
    data = cluster_center,
    aes(label = cluster_name, fill = cluster_name),
#     fontface = 'bold', 
    size = 2, 
    box.padding = unit(0.3, "lines"),
    point.padding = unit(0.2, "lines"),
    segment.color = 'grey50'
  ) +
  scale_fill_manual(values = c(colors37, colors37, colors37), name = "") +
  theme_bw(base_size = 10) +
  theme(
    legend.position = "none",
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(color="black", size=10, face = "italic")
)

Warning message:
“ggrepel: 31 unlabeled data points (too many overlaps). Consider increasing max.overlaps”


## Observations:
# - Cells are clustered into major cell types based on transcriptomic and proteomic differences;
# - Protein surface markers are helpful - more clearly seperate NK cells from CD8 T cells; seperate DCs from macrophages, etc
# - Cell type-specific analysis are done on each cell type with fine-grained clustering. Cell type-speicfic references are provided as well.


# In order to map a new query dataset onto the reference, you will need a reference object saved from the steps above, as well as query cell expression and metadata.
# The query dataset is assumed to have been normalized in the same manner as the reference cells (here, default is log(CP10k+1) normalization).


query_exp <- readRDS("/data/srlab/fzhang/amp/data/amp_ra_phase1/celseq_synovium_log2_5265cells_paper.rds")
query_metadata <- readRDS("/data/srlab/fzhang/amp/data/amp_ra_phase1/celseq_synovium_meta_5265cells_paper.rds")
colnames(query_metadata)[2] <- "cell"
query_metadata <- query_metadata[, c("cell", "sample", "cell_type", "disease", "plate", "cluster", "Mahalanobis_20")]
dim(query_exp)
dim(query_metadata)
query_metadata[1:4,]


# Map query
query = mapQuery(query_exp,             # query gene expression (genes x cells)
                 query_metadata,        # query metadata (cells x attributes)
                 reference,             # Symphony reference object
                 vars = c('sample'),    # also integrate over query batches
                 do_normalize = FALSE,  # perform log(CP10k) normalization on query
                 do_umap = TRUE)        # project query cells into reference UMAP

Scaling and synchronizing query gene expression
Found 4676 reference variable genes in query dataset
Project query cells using reference gene loadings
Clustering query cells to reference centroids
Correcting query batch effects
UMAP
All done!


# Predict query cell types using k-NN
query = knnPredict(query, reference, reference$meta_data$cell_type, k = 5)


# head(query$meta_data)


# Visualize of mapping
# Sync the column names for both data frames
reference$meta_data$cell_type_pred_knn = NA
reference$meta_data$ref_query = 'reference'
query$meta_data$ref_query = 'query'

# Add the UMAP coordinates to the metadata
meta_data_combined = rbind(query$meta_data[, c("cell", "sample", "cell_type", "cell_type_pred_knn", "ref_query")], reference$meta_data[, c("cell", "sample", "cell_type", "cell_type_pred_knn", "ref_query")])
umap_combined = rbind(query$umap, reference$umap$embedding)
umap_combined_labels = cbind(meta_data_combined, umap_combined)


# Plot UMAP visualization of all cells
# plotBasic(umap_combined_labels, title = 'Reference and query cells', color.by = 'ref_query')

options(repr.plot.height = 4, repr.plot.width = 6)
ggplot() +
  geom_point(
    data = umap_combined_labels[rev(order(umap_combined_labels$ref_query)),],
    mapping = aes_string(x = "UMAP1", y = "UMAP2", fill = "ref_query"),
    size = 0.2, stroke = 0.0001, shape = 21
  ) +
#   facet_wrap( ~ new_cluster_name) +
#   scale_fill_manual(values = meta_colors$cell_type, name = "") + 
  theme_bw(base_size = 10) +
  theme(
#     legend.position = "none",
#     axis.text = element_blank(),
#     axis.ticks = element_blank(),
#     panel.grid = element_blank(),
    plot.title = element_text(color="black", size=10)
  )


options(repr.plot.height = 3, repr.plot.width = 7)
plotBasic(umap_combined_labels, title = 'Reference and query cells', 
          color.by = 'cell_type', facet.by = 'ref_query', color.mapping = meta_colors$cell_type)


## Cell type mappings heatmap: compare real cell type labels with predicted labels
res = symphony:::evaluate(query$meta_data$cell_type_pred_knn, query$meta_data$cell_type)

# # fig.size(4,5)
Conf = t(res$Conf) / rowSums(t(res$Conf))
Conf <- Conf[complete.cases(Conf), ] # Remove NA value rows becuase ther was no endothelial cells and no NK cells in the query

options(repr.plot.height = 3.5, repr.plot.width = 4)
pheatmap(Conf)


## Summary:
# - We are able to map the query cells into the right cell types in the reference panel.
# - No endothelial cells in the query, we the query cells are not mapped to the endothelial cells in the reference.
# - A small number of CD8 T cells in the query are mapped to the NK cells in the reference, this might be acceptable since these CD8 T cells are cytotoxic so share many similar genes with NK cells.
# - To more accurately map cells onto subpopulations, we recommend mapping query cells onto cell type-specific references.


sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server release 6.5 (Santiago)

Matrix products: default
BLAS/LAPACK: /PHShome/fz049/miniconda3/lib/libopenblasp-r0.3.9.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggrastr_0.2.1           symphony_1.0            ggthemes_4.2.0         
 [4] CCA_1.2                 fields_11.6             spam_2.5-1             
 [7] dotCall64_1.0-0         fda_5.1.5.1             Seurat_3.9.9.9027      
[10] rstatix_0.6.0           ggpubr_0.4.0            forcats_0.5.1          
[13] stringr_1.4.0           purrr_0.3.4             readr_1.4.0            
[16] tidyr_1.1.2             tibble_3.1.0            tidyverse_1.3.0        
[19] umap_0.2.6.0            singlecellmethods_0.1.0 viridis_0.5.1          
[22] viridisLite_0.3.0       matrixStats_0.58.0      MASS_7.3-53            
[25] RColorBrewer_1.1-2      plyr_1.8.6              ggrepel_0.9.1          
[28] ggbeeswarm_0.6.0        ggplot2_3.3.3           reticulate_1.18        
[31] gridExtra_2.3           magrittr_2.0.1          harmony_1.0            
[34] uwot_0.1.10             Matrix_1.3-2            presto_1.0.0           
[37] data.table_1.13.6       Rcpp_1.0.6              patchwork_1.1.1        
[40] dplyr_1.0.4             readxl_1.3.1            gdata_2.18.0           

loaded via a namespace (and not attached):
  [1] utf8_1.1.4            tidyselect_1.1.0      htmlwidgets_1.5.3    
  [4] Rtsne_0.15            munsell_0.5.0         codetools_0.2-16     
  [7] ica_1.0-2             pbdZMQ_0.3-3          future_1.21.0        
 [10] miniUI_0.1.1.1        withr_2.4.1           colorspace_2.0-0     
 [13] uuid_0.1-2            rstudioapi_0.13       ROCR_1.0-7           
 [16] ggsignif_0.6.0        tensor_1.5            listenv_0.8.0        
 [19] labeling_0.4.2        repr_0.19.2           polyclip_1.10-0      
 [22] farver_2.1.0          parallelly_1.23.0     vctrs_0.3.6          
 [25] generics_0.1.0        R6_2.5.0              rsvd_1.0.2           
 [28] bitops_1.0-6          spatstat.utils_1.20-2 assertthat_0.2.1     
 [31] promises_1.1.1        scales_1.1.1          beeswarm_0.2.3       
 [34] gtable_0.3.0          Cairo_1.5-12.2        globals_0.14.0       
 [37] goftest_1.2-2         rlang_0.4.10          splines_3.6.3        
 [40] lazyeval_0.2.2        broom_0.7.4.9000      reshape2_1.4.4       
 [43] abind_1.4-5           modelr_0.1.8          backports_1.2.1      
 [46] httpuv_1.5.5          tools_3.6.3           ellipsis_0.3.1       
 [49] gplots_3.0.1.1        ggridges_0.5.3        base64enc_0.1-3      
 [52] ps_1.5.0              rpart_4.1-15          openssl_1.3          
 [55] deldir_0.1-23         pbapply_1.4-2         cowplot_1.1.1        
 [58] zoo_1.8-6             haven_2.4.1           cluster_2.1.0        
 [61] fs_1.5.0              RSpectra_0.16-0       scattermore_0.7      
 [64] openxlsx_4.2.3        lmtest_0.9-36         reprex_0.3.0         
 [67] RANN_2.6.1            fitdistrplus_1.1-3    hms_1.0.0            
 [70] mime_0.9              evaluate_0.14         xtable_1.8-4         
 [73] rio_0.5.16            compiler_3.6.3        maps_3.3.0           
 [76] KernSmooth_2.23-15    crayon_1.4.1          htmltools_0.5.1.1    
 [79] mgcv_1.8-33           later_1.1.0.1         lubridate_1.7.4      
 [82] DBI_1.0.0             dbplyr_2.0.0          car_3.0-2            
 [85] cli_2.3.1             parallel_3.6.3        igraph_1.2.6         
 [88] pkgconfig_2.0.3       foreign_0.8-71        IRdisplay_0.7.0      
 [91] plotly_4.9.0          xml2_1.3.2            vipor_0.4.5          
 [94] rvest_0.3.6           digest_0.6.27         sctransform_0.3.2    
 [97] RcppAnnoy_0.0.18      spatstat.data_1.7-0   cellranger_1.1.0     
[100] leiden_0.3.6          curl_4.3              shiny_1.6.0          
[103] gtools_3.8.1          lifecycle_1.0.0       nlme_3.1-151         
[106] jsonlite_1.7.2        carData_3.0-4         askpass_1.0          
[109] fansi_0.4.2           pillar_1.5.1          lattice_0.20-41      
[112] fastmap_1.1.0         httr_1.4.2            survival_2.44-1.1    
[115] glue_1.4.2            zip_2.1.1             spatstat_1.64-1      
[118] png_0.1-7             class_7.3-18          stringi_1.5.3        
[121] caTools_1.17.1.2      IRkernel_0.8.15       irlba_2.3.3          
[124] future.apply_1.7.0

symbol	mean	stddev
IGHA1	0.3966540	0.7501617
IGLC2	0.8484138	1.1498599
IGHG2	0.0721903	0.3591387
IGLC3	0.4520035	0.8971236
JCHAIN	0.2230015	0.6807754
SPP1	0.3448693	0.9752211

	used	(Mb)	gc trigger	(Mb)	max used	(Mb)
Ncells	3191131	170.5	6467374	345.4	6467374	345.4
Vcells	1978966319	15098.4	4128779532	31500.1	9545045347	72823.0

	cell	sample	cell_type	disease	plate	cluster	Mahalanobis_20
4264	S006_L1Q1_A03	301-0159	T cell	OA	S006	SC-T1	OA
4265	S006_L1Q1_A05	301-0159	T cell	OA	S006	SC-T4	OA
4266	S006_L1Q1_A07	301-0159	T cell	OA	S006	SC-T1	OA
4267	S006_L1Q1_A09	301-0159	T cell	OA	S006	SC-T1	OA

Build a reference for AMP phase II RA CITE-seq, and map AMP phase I RA scRNA-seq to this reference¶

Fan Zhang¶

2021-05-14¶

Load QCed cells from mRNA and protein matrices¶

Build reference¶

Visualize results¶

Visualize cell type lineage protein markers¶

Protein expression of cell type lineage markers help idenitfy key subpopulations¶

Now let's project the fine grained clusters defined by cell type-specific analysis onto coarse anlaysis/global embeddings:¶

Map query¶

Map AMP phase I RA cells (query) generated from CEL-Seq2 to the built reference¶

Reference: AMP phase II; Query: AMP phase I¶