[1]:

library(Seurat)
library(gridExtra)
library(Matrix)
library(ggplot2)

Mye harmony#

[2]:

library(harmony)

Loading required package: Rcpp


[3]:

library(pheatmap)

[27]:

library(future)

[4]:

source('/home/public/ghx/software/my_seurat.r')

[5]:

RDSPath <- '/home/public/ghx/wuLab/Bladder/1217_all/seu'

[6]:

ct <- 'Mye'

[6]:

mtx <- readRDS(file = file.path(RDSPath,
                        paste0('V1_Celltype_',ct,'_matrix.rds')))
anno <- readRDS(file = file.path(RDSPath,
                        paste0('V1_Celltype_',ct,'_annotation.rds')))

[7]:

seu <- CreateSeuratObject(counts = mtx, meta.data = anno, project = ct)

[8]:

samples <- levels(seu@meta.data$orig.ident)

[9]:

table(data.frame(origin = seu@meta.data$orig.ident))


 MIBC5P  MIBC5T  MIBC6P  MIBC6T  MIBC7P  MIBC7T  MIBC8P  MIBC8T  MIBC9P  MIBC9T
    418     206     771     538    3297     803     795     165     698     489
MIBC10P MIBC10T MIBC11P MIBC11T MIBC12P MIBC12T MIBC13P MIBC13T MIBC14P MIBC14T
    612     288    3162    1492     909    2800     654     326     513     125
MIBC15P MIBC15T MIBC16P MIBC16T
    167    1553     423     315

First, integrate all T cell raw data without doing any batch effect removal#

[10]:

seu <- Nml_Scl_Fvg_PCA_Fnb(seu = seu)

[11]:

seu <- FindNeighbors(seu, reduction = 'pca', dims = 1:20, verbose = F)

[12]:

seu <- FindClusters(object = seu,
            resolution = 0.5, verbose = F)

[13]:

seu <- RunUMAP(object = seu, dims = 1:25, verbose = F,
               umap.method = 'umap-learn', metric = 'correlation')

[93]:

table(seu@active.ident)


   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
5214 3082 2363 1747 1690 1408  992  878  849  742  583  545  452  372  354  155
  16
  93

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
p1 <- DimPlot(seu, label = T, pt.size = 0.1) +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'orig.ident') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'status') +
     labs(title = 'status') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3), ncol = 2,
                top = paste0('MIBC5-16: raw merge'))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                            'Merge_v1_Mye_raw_overview.jpg'))
grid.arrange(p)

17.png 18.png

[100]:

seu@meta.data$Mye.ct <- 'Mac'
seu@meta.data$Mye.ct[seu@active.ident %in% c(14)] <- 'cDC'
seu@meta.data$Mye.ct[seu@active.ident %in% c(15)] <- 'pDC'
seu@meta.data$Mye.ct[seu@active.ident %in% c(0,2,12)] <- 'Neu'

[101]:

saveRDS(seu@meta.data,
        file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                         'V1_Celltype_Mac_annotation.rds'))

[ ]:

options(repr.plot.width=10, repr.plot.height=20, repr.plot.res = 200)
FeaturePlot(object = seu, features = c('CD68', 'HLA-DRA',
                     'CD14', 'FCGR3A','CLEC9A', 'XCR1','LILRA4','TCF4'),
            reduction = 'umap', ncol = 2,
            pt.size = 0.01)

19.png 20.png

[16]:

options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
FeaturePlot(object = seu, features = c('CSF3R', 'SIGLEC5','LILRA5', 'S100A8'),
            reduction = 'umap', ncol = 2,
            pt.size = 0.01)
../../_images/notebooks_Cell_type_analysis_Mye_harmony_24_0.png

Harmony: remove batch effect between patients#

[31]:

options(repr.plot.width=5, repr.plot.height=5, repr.plot.res = 100)
seu <- RunHarmony(seu, theta = 10, verbose = F,
                  group.by.vars = 'patient', plot_convergence = T)
../../_images/notebooks_Cell_type_analysis_Mye_harmony_26_0.png 
[26]:

seu <- RunUMAP(object = seu, reduction.name = 'umap_harmony',
               umap.method = 'umap-learn', metric = 'correlation',
               reduction = 'harmony', dims = 1:20)

Warning message:
“Cannot add objects with duplicate keys (offending key: UMAP_), setting key to 'umap_harmony_'”

[27]:

seu <- FindNeighbors(seu, reduction = 'harmony', dims = 1:20, verbose = F)

[28]:

seu <- FindClusters(seu, verbose = F,resolution = 0.5)

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
p1 <- DimPlot(seu, label = T, pt.size = 0.1,
              reduction = 'umap_harmony') +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'orig.ident',
              reduction = 'umap_harmony') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'status',
              reduction = 'umap_harmony') +
     labs(title = 'status') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'patient',
              reduction = 'umap_harmony') +
     labs(title = 'patient') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16: raw merge'))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                            'Merge_v1_Mye_harmony_patient_overview.jpg'))
grid.arrange(p)

21.png 22.png

[ ]:

options(repr.plot.width=10, repr.plot.height=20, repr.plot.res = 200)
FeaturePlot(object = seu, features = c('CD68', 'HLA-DRA',
                     'CD14', 'FCGR3A','CLEC9A', 'XCR1','LILRA4','TCF4'),
            reduction = 'umap_harmony', ncol = 2,
            pt.size = 0.01)

image.png image.png

use batch-effect-removal embedding(harmony) to reconstruct scaled expression data#

[9]:

# X=U*Sigma*V, U and V are Unitary Matrix, e.g. V^-1=t(V)
# U=X*t(V)*(1/Sigma)
# X(nGene*nCell)=U(nGene*nPC)*Sigma(nPC*nPC)*V(nPC*nCell)
#X <- seu@assays$RNA@scale.data
Sigma <- seu@reductions$pca@stdev
V <- seu@reductions$pca@cell.embeddings[, 1:100]
V <- apply(V,MARGIN = 2,FUN = function(x){
    return(x/sqrt(sum(x^2)))
})
# here V(nCell*nPC) is actually t(V)
U <- seu@assays$RNA@scale.data %*% V %*% diag(1/Sigma)

[10]:

# after batch effect removal, X.new=U*Sigma*V.new
V.new <- seu@reductions$harmony@cell.embeddings[, 1:100]
V.new <- apply(V.new,MARGIN = 2,FUN = function(x){
    return(x/sqrt(sum(x^2)))
})
#clean.data.full <- U %*% diag(Sigma) %*% t(V.new)

[11]:

seu <- SetAssayData(seu, slot = 'scale.data', new.data = U %*% diag(Sigma) %*% t(V.new))

[ ]:

options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
p <- DoHeatmap(object = seu, cells = sub.cells,draw.lines = FALSE,
          features = top10$gene, size = 5, angle = 0, disp.min = -2, disp.max = 2) +
theme(text = element_text(size = 4))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/T/Merge',
                            'Merge_v1_T_harmony_patient_recons_heatmap.jpg'))
p

23.png

[ ]:

library(future)
plan("multiprocess", workers = 8)
options(future.globals.maxSize = 10000 * 1024^2)
markers.new <- FindAllMarkers(seu, verbose = F, slot = 'scale.data')
plan("sequential")

[ ]:

top10 <- extractTopNMarkers(df = markers.new, top_n = 20)

[ ]:

options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
p <- DoHeatmap(object = seu, cells = sub.cells,draw.lines = FALSE,
          features = top10$gene, size = 5, angle = 0, disp.min = -2, disp.max = 2) +
theme(text = element_text(size = 4))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/T/Merge',
                            'Merge_v1_T_harmony_patient_recons_heatmap.jpg'))
p

[ ]:

write.csv(markers.new,
          file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/T/Merge',
                    'Merge_v1_T_harmony_patient_reconstuct_cluster_markers.csv'))

only take Macrophages out to do harmony#

[15]:

Mac.clusters <- c(1,3,4,5,6,7,8,9,10,11,13,16)
temp.cells <- Cells(seu)[ seu@active.ident %in% Mac.clusters]
temp.mtx <- mtx[, temp.cells]
temp.anno <- anno[temp.cells, ]
seu.sub <- CreateSeuratObject(counts = temp.mtx, meta.data = temp.anno, project = 'Mac')

[16]:

seu.sub <- Nml_Scl_Fvg_PCA_Fnb(seu = seu.sub)

[17]:

seu.sub <- FindNeighbors(seu.sub, reduction = 'pca', dims = 1:20, verbose = F)

[18]:

seu.sub <- FindClusters(object = seu.sub,
            resolution = 0.5, verbose = F)

[19]:

seu.sub <- RunUMAP(object = seu.sub, dims = 1:25, verbose = F,
               umap.method = 'umap-learn', metric = 'correlation')

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
p1 <- DimPlot(seu.sub, label = T, pt.size = 0.3) +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu.sub, label = T, pt.size = 0.3, group.by = 'orig.ident') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu.sub, label = T, pt.size = 0.3, group.by = 'status') +
     labs(title = 'status') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3), ncol = 2,
                top = paste0('MIBC5-16: raw merge'))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                            'Merge_v1_Mac_raw_overview.jpg'))
grid.arrange(p)

24.png 25.png

[23]:

options(repr.plot.width=5, repr.plot.height=5, repr.plot.res = 100)
seu.sub <- RunHarmony(seu.sub, theta = 10, verbose = F,
                  group.by.vars = 'patient', plot_convergence = T)
../../_images/notebooks_Cell_type_analysis_Mye_harmony_52_0.png 
[24]:

seu.sub <- RunUMAP(object = seu.sub, reduction.name = 'umap_harmony',
               umap.method = 'umap-learn', metric = 'correlation',
               reduction = 'harmony', dims = 1:20)

Warning message:
“Cannot add objects with duplicate keys (offending key: UMAP_), setting key to 'umap_harmony_'”

[25]:

seu.sub <- FindNeighbors(seu.sub, reduction = 'harmony', dims = 1:20, verbose = F)

[26]:

seu.sub <- FindClusters(seu.sub, verbose = F,resolution = 1.0)

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
p1 <- DimPlot(seu.sub, label = T, pt.size = 0.5,
              reduction = 'umap_harmony') +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu.sub, label = T, pt.size = 0.5, group.by = 'orig.ident',
              reduction = 'umap_harmony') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu.sub, label = T, pt.size = 0.5, group.by = 'status',
              reduction = 'umap_harmony') +
     labs(title = 'status') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu.sub, label = T, pt.size = 0.5, group.by = 'patient',
              reduction = 'umap_harmony') +
     labs(title = 'patient') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16: raw merge'))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                            'Merge_v1_Mac_harmony_overview.jpg'))
grid.arrange(p)

26.png 27.png

[ ]:

options(repr.plot.width=10, repr.plot.height=20, repr.plot.res = 200)
FeaturePlot(object = seu.sub, features = c('FCER1A', 'PTGS2',
                     'CD14', 'FCGR3A','CCL5', 'MRC1','MRC1','TCF4'),
            reduction = 'umap_harmony', ncol = 2,
            pt.size = 0.01)

29.png 30.png

[49]:

marker.list <- list()
marker.list[['M1']] <- c('CD163','CD80', 'CD86', 'CD36', 'CD68',
                       'FCGR2A', 'FCGR3A', 'IFNGR1','HLA-DRA', 'PTGS2',
                       'NOS2', 'IRF5', 'STAT1', 'CCL2', 'CCL3',
                       'CCL4', 'CCL5', 'CCL8', 'CCL11', 'CCL15',
                       'CCL19', 'CCL20', 'CXCL1', 'CXCL2', 'CXCL3',
                       'CXCL5', 'CXCL8', 'CXCL9', 'CXCL10', 'CXCL11',
                       'CXCL13', 'CXCL1', 'IL1B', 'IL6', 'IL12A',
                       'IL15', 'IL17A', 'IL18', 'IL23A', 'IFNG', 'TNF')
marker.list[['M2a']] <- c('CD163', 'CD200R1', 'CLEC10A', 'CXCR1', 'CXCR2',
                        'CD209', 'CLEC7A', 'FCER1A', 'IL1R2', 'IL4R',
                        'MRC1', 'HLA-DRA', 'IRF4', 'PPARG', 'STAT6',
                        'CCL1', 'CCL2', 'CCL14', 'CCL17', 'CCL18',
                        'CCL22', 'CCL23', 'CCL24', 'CCL26', 'IL1R1',
                        'IL10', 'IL12A', 'TGFB1','CD86', 'IL4R')
marker.list[['M2b']] <- c('PTGS2', 'IRF4', 'SOCS3', 'SPHK1', 'SPHK2', 'CCL1',
                        'CCL20', 'CXCL1', 'CXCL2', 'CXCL3', 'CSF3',
                        'CSF2', 'IL1B', 'IL6', 'IL10', 'TNF')
marker.list[['M2c']] <- c('CCR2', 'SLAMF1', 'CD163', 'IL4R', 'MRC1',
                        'MSR1', 'SCARB1', 'TLR1', 'ARG1', 'IRF4', 'SOCS3', 'TLR8')
marker.list[['M2d']] <- c('CCL5', 'CXCL10', 'CXCL16', 'IL10', 'IL12A',
                        'TNF', 'VEGFA')
marker.list[['HCM']] <- c('FCGR3A', 'CX3CR1', 'ITGAM', 'CD163', 'HLA-DRA',
                        'CSF1R', 'SELL', 'CD14', 'CCR2')
marker.list[['HIM']] <- c('CCR2', 'ITGAM', 'FCGR3A', 'CD163', 'CCR5', 'CSF1R',
                        'SELL', 'CD14', 'CX3CL1', 'HLA-DRA')
marker.list[['HNCM']] <- c('CCR2', 'CD163', 'SELL', 'ITGAM', 'CD14', 'CSF1R',
                         'FCGR3A', 'CX3CL1', 'HLA-DRA')
marker.list[['GMDSC']] <- c('HLA-DRA', 'FCGR3A', 'IL10', 'MMP9', 'TGFB1', 'VEGFA')
marker.list[['MMDSC']] <- c('FUT4', 'CEACAM8','HLA-DRA','ITGAM', 'CD14',
                          'CD33', 'IL4R','IL10', 'MMP9', 'TGFB1', 'VEGFA')

[62]:

resultPath <- '/home/public/ghx/wuLab/Bladder/0119_discussion/FeaturePlot'
for(ct in names(marker.list)){
    temp.markers <- marker.list[[ct]]
    gene.nb <- length(temp.markers)
    for(i in 0:(gene.nb %/% 8)){
        features.plot <- temp.markers[(i*8+1):min(gene.nb,(i+1)*8)]
        if(length(features.plot)==0) {next}
        if(length(features.plot)==1){
            temp.height <- 10
        }else{
            temp.height <- 5 * ((length(features.plot)+1) %/% 2)
        }
        p <- FeaturePlot(object = seu.sub, features = features.plot,
                         pt.size = 0.1, ncol = 2)
        ggsave(plot = p, limitsize = FALSE,
               filename = file.path(resultPath,'Mac',
                                    paste0('Mac_',ct,'_',i,'_FeaturePlot.jpg')),
               width = 10, height = temp.height)
        p <- VlnPlot(object = seu.sub, features = features.plot,
                         pt.size = 0.1, ncol = 2)
        ggsave(plot = p, limitsize = FALSE,
               filename = file.path(resultPath,'Mac',
                                    paste0('Mac_',ct,'_',i,'_VlnPlot.jpg')),
               width = 10, height = temp.height)
    }
}

Warning message in FetchData(object = object, vars = c(dims, "ident", features), :
“The following requested variables were not found: NA”
Warning message in FetchData(object = object, vars = features, slot = slot):
“The following requested variables were not found: NA”

[29]:

library(future)
plan("multiprocess", workers = 12)
options(future.globals.maxSize = 10000 * 1024^2)
markers <- FindAllMarkers(seu.sub, verbose = F)
plan("sequential")

[ ]:

options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
top10 <- extractTopNMarkers(df = markers, top_n = 20)
p <- DoHeatmap(object = seu.sub, draw.lines = FALSE,
          features = top10$gene, size = 5, angle = 0,
               disp.min = -2, disp.max = 2) +
theme(text = element_text(size = 4))
p

28.png

[31]:

ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                            'Merge_v1_Mac_heatmap_cluster.jpg'))

[32]:

write.csv(markers, file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                            'Mac_markers_cluster.csv'))

[33]:

stats <- table(seu.sub@active.ident, seu.sub@meta.data$orig.ident)
stats <- apply(stats, MARGIN = 2,
               FUN = function(x) {return(x/sum(x))})
stats <- t(apply(stats, MARGIN = 1,
               FUN = function(x) {return(x/max(x))}))
samples <- colnames(stats)
patient <- substr(samples, start = 5,
                  stop = sapply(samples,FUN=nchar)-1)
patient <- factor(patient, levels=5:16)
status <- substr(samples, start = sapply(samples,FUN=nchar),
                  stop = sapply(samples,FUN=nchar))
status <- factor(status, levels=c('P','T'))
row.anno <- data.frame(row.names = colnames(stats),
                       patient = patient, status = status)

[34]:

p <- pheatmap(t(stats), cutree_rows = 6, cutree_cols = 5,
              filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                                paste0('Merge_v1_Mac_heatmap_sample_cluster.jpg')),
         annotation_row = row.anno)
options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
p

[35]:

p
../../_images/notebooks_Cell_type_analysis_Mye_harmony_69_0.png 
[36]:

# X=U*Sigma*V, U and V are Unitary Matrix, e.g. V^-1=t(V)
# U=X*t(V)*(1/Sigma)
# X(nGene*nCell)=U(nGene*nPC)*Sigma(nPC*nPC)*V(nPC*nCell)
#X <- seu@assays$RNA@scale.data
Sigma <- seu.sub@reductions$pca@stdev
V <- seu.sub@reductions$pca@cell.embeddings[, 1:100]
V <- apply(V,MARGIN = 2,FUN = function(x){
    return(x/sqrt(sum(x^2)))
})
# here V(nCell*nPC) is actually t(V)
U <- seu.sub@assays$RNA@scale.data %*% V %*% diag(1/Sigma)
# after batch effect removal, X.new=U*Sigma*V.new
V.new <- seu.sub@reductions$harmony@cell.embeddings[, 1:100]
V.new <- apply(V.new,MARGIN = 2,FUN = function(x){
    return(x/sqrt(sum(x^2)))
})
seu.sub <- SetAssayData(seu.sub, slot = 'scale.data',
                    new.data = U %*% diag(Sigma) %*% t(V.new))

[ ]:

top10 <- extractTopNMarkers(df = markers, top_n = 20)
p <- DoHeatmap(object = seu.sub, draw.lines = FALSE,
          features = top10$gene, size = 5, angle = 0,
               disp.min = -2, disp.max = 2) +
theme(text = element_text(size = 4))
options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
p

31.png

[72]:

cell.type <- rep(NA,length(Cells(seu.sub)))
names(cell.type) <- Cells(seu.sub)
cell.type[seu.sub@active.ident == 0] <- 'CD16-CD14-(S100A)'
cell.type[seu.sub@active.ident == 1] <- 'CD16++CD14+(apoliprotein)'
cell.type[seu.sub@active.ident == 2] <- 'CD16-CD14-(S100A)'
cell.type[seu.sub@active.ident == 3] <- 'CD1c+CD1e+DC(FCER1A)'
cell.type[seu.sub@active.ident == 4] <- 'CD16-CD14+(CCL20)'
cell.type[seu.sub@active.ident == 5] <- 'CD16-CD14-(HSP)'
cell.type[seu.sub@active.ident == 6] <- 'CD16-CD14-(S100B)'
cell.type[seu.sub@active.ident == 7] <- 'CD16++CD14-(CDKN1C)'
cell.type[seu.sub@active.ident == 8] <- 'CD16+CD14+(CCL4)'
cell.type[seu.sub@active.ident == 9] <- 'CD16++CD14+(IFNresponding)'
cell.type[seu.sub@active.ident == 10] <- 'CD16+CD14++(SPP1)'
cell.type[seu.sub@active.ident == 11] <- 'CD16+CD14++(RNASE1)'
cell.type[seu.sub@active.ident == 12] <- 'DOUBLETS(T&MAC)'
cell.type[seu.sub@active.ident == 13] <- 'CD16++CD14+(Complement)'
cell.type[seu.sub@active.ident == 14] <- 'CD1e+DC(TXN)'
seu.sub@meta.data$anno.v0 <- cell.type

[75]:

seu.sub <- subset(seu.sub, cells = Cells(seu.sub)[ seu.sub@active.ident!=12 ])

Warning message:
“All object keys must be alphanumeric characters, followed by an underscore ('_'), setting key to 'umapharmony_'”

[76]:

saveRDS(seu.sub
        file = file.path(RDSPath,
                         'V1_Celltype_Mac_seu.rds'))

[3]:

seu.sub <- readRDS(file.path(RDSPath,
                    'V1_Celltype_Mac_seu.rds'))

[5]:

saveRDS(seu.sub@meta.data, file = file.path(RDSPath,
                                  'V1_Celltype_Mac_annotation_withAnno.rds'))

[ ]:

DimPlot(seu.sub, label = T, pt.size = 1, reduction = 'umap_harmony',
              group.by = 'anno.v0') +
     labs(title = 'cell.type.v1') + theme(legend.position = 'none')

32.png

[104]:

ct.table <- table(seu.sub@meta.data$orig.ident, seu.sub@active.ident)
write.csv(ct.table, file.path('/home/public/ghx/wuLab/Bladder/0129_Mac',
                              'Mac_cluster_by_Patient.csv'))

only take cDC out to do harmony#

[79]:

cDC.clusters <- c(14)
temp.cells <- Cells(seu)[ seu@active.ident %in% cDC.clusters]
temp.mtx <- mtx[, temp.cells]
temp.anno <- anno[temp.cells, ]
seu.sub <- CreateSeuratObject(counts = temp.mtx,
                              meta.data = temp.anno,
                              project = 'cDC')

[80]:

seu.sub <- Nml_Scl_Fvg_PCA_Fnb(seu = seu.sub)

[81]:

options(repr.plot.width=5, repr.plot.height=5, repr.plot.res = 100)
seu.sub <- RunHarmony(seu.sub, theta = 10, verbose = F,
                  group.by.vars = 'patient', plot_convergence = T)
../../_images/notebooks_Cell_type_analysis_Mye_harmony_84_0.png 
[82]:

seu.sub <- RunUMAP(object = seu.sub, reduction.name = 'umap_harmony',
               umap.method = 'umap-learn', metric = 'correlation',
               reduction = 'harmony', dims = 1:20)

[83]:

seu.sub <- FindNeighbors(seu.sub, reduction = 'harmony', dims = 1:20, verbose = F)

[86]:

seu.sub <- FindClusters(seu.sub, verbose = F,resolution = 1.2)

[90]:

options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
p1 <- DimPlot(seu.sub, label = T, pt.size = 1,
              reduction = 'umap_harmony') +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu.sub, label = T, pt.size = 1, group.by = 'orig.ident',
              reduction = 'umap_harmony') +
     labs(title = 'sample.origin') + theme(legend.position = 'none')
p3 <- DimPlot(seu.sub, label = T, pt.size = 1, group.by = 'status',
              reduction = 'umap_harmony') +
     labs(title = 'status') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu.sub, label = T, pt.size = 1, group.by = 'patient',
              reduction = 'umap_harmony') +
     labs(title = 'patient') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16: cDC merge'))
ggsave(plot = p, limitsize = FALSE, width = 10, height = 10,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                            'Merge_v1_cDC_harmony_overview.jpg'))
grid.arrange(p)
../../_images/notebooks_Cell_type_analysis_Mye_harmony_88_0.png 
[88]:

markers <- FindAllMarkers(seu.sub, verbose = F)

[ ]:

top10 <- extractTopNMarkers(df = markers, top_n = 25)
p <- DoHeatmap(object = seu.sub, draw.lines = FALSE,
          features = top10$gene, size = 5, angle = 0,
               disp.min = -2, disp.max = 2) +
        theme(text = element_text(size = 4))
ggsave(plot = p, limitsize = FALSE, width = 10, height = 10,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                                'Merge_v1_cDC_harmony_heatmap.jpg'))
options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
p

image.png

[7]:

write.csv(markers,file.path('/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge',
                            'cDC_markers_cluster.csv'))

Error in is.data.frame(x): object 'markers' not found
Traceback:

1. write.csv(markers, file.path("/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge",
 .     "cDC_markers_cluster.csv"))
2. eval.parent(Call)
3. eval(expr, p)
4. eval(expr, p)
5. write.table(markers, file.path("/home/public/ghx/wuLab/Bladder/1217_all/CellTypes/Mye/Merge",
 .     "cDC_markers_cluster.csv"), col.names = NA, sep = ",", dec = ".",
 .     qmethod = "double")
6. is.data.frame(x)

calculate batch effect removed data and combine with neu&DC#

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
p1 <- DimPlot(seu.sub, label = T, pt.size = 0.1,
              reduction = 'umap_harmony') +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu.sub, label = T, pt.size = 0.1, group.by = 'orig.ident',
              reduction = 'umap_harmony') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu.sub, label = T, pt.size = 0.1, group.by = 'status',
              reduction = 'umap_harmony') +
     labs(title = 'status') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu.sub, label = T, pt.size = 0.1, group.by = 'patient',
              reduction = 'umap_harmony') +
     labs(title = 'patient') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16: raw merge'))
grid.arrange(p)

33.png 34.png

[ ]:



[262]:

stats <- table(seu.sub@active.ident, seu.sub@meta.data$orig.ident)
stats <- apply(stats, MARGIN = 2,
               FUN = function(x) {return(x/sum(x))})
stats <- t(apply(stats, MARGIN = 1,
               FUN = function(x) {return(x/max(x))}))
samples <- colnames(stats)
patient <- substr(samples, start = 5,
                  stop = sapply(samples,FUN=nchar)-1)
patient <- factor(patient, levels=5:16)
status <- substr(samples, start = sapply(samples,FUN=nchar),
                  stop = sapply(samples,FUN=nchar))
status <- factor(status, levels=c('P','T'))
row.anno <- data.frame(row.names = colnames(stats),
                       patient = patient, status = status)

[263]:

p <- pheatmap(t(stats), cutree_rows = 6, cutree_cols = 5,
         annotation_row = row.anno)
options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
p
../../_images/notebooks_Cell_type_analysis_Mye_harmony_98_0.png 
[170]:

# X=U*Sigma*V, U and V are Unitary Matrix, e.g. V^-1=t(V)
# U=X*t(V)*(1/Sigma)
# X(nGene*nCell)=U(nGene*nPC)*Sigma(nPC*nPC)*V(nPC*nCell)
#X <- seu@assays$RNA@scale.data
Sigma <- seu@reductions$pca@stdev
V <- seu@reductions$pca@cell.embeddings[, 1:100]
V <- apply(V,MARGIN = 2,FUN = function(x){
    return(x/sqrt(sum(x^2)))
})
# here V(nCell*nPC) is actually t(V)
U <- seu@assays$RNA@scale.data %*% V %*% diag(1/Sigma)
seu <- SetAssayData(seu, slot = 'scale.data',
                    new.data = U %*% diag(Sigma) %*% t(V))

[172]:

bind.scaleData <- cbind(seu.sub@assays$RNA@scale.data,
                        seu@assays$RNA@scale.data[,setdiff(Cells(seu),Cells(seu.sub))])
bind.scaleData <- bind.scaleData[, Cells(seu)]
seu <- SetAssayData(seu, slot = 'scale.data',
                    new.data = bind.scaleData)

[208]:

seu <- RunPCA(seu, npcs = 100, verbose = F)
seu <- FindNeighbors(seu,dims = 1:20,verbose = F)

[212]:

seu <- RunUMAP(object = seu,
               reduction.name = 'umap_harmony_cData',
               umap.method = 'umap-learn', metric = 'correlation',
               dims = 1:20)

[213]:

seu <- FindClusters(seu, verbose = F,resolution = 0.3)

[232]:

my.anno <- as.character(seu@active.ident)
my.anno[ my.anno==6 & Cells(seu) %in% Cells(seu.sub)] <- 11
seu@meta.data$my.anno <- factor(my.anno, levels = 0:11)

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
p1 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'my.anno',
              reduction = 'umap_harmony_cData') +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'orig.ident',
              reduction = 'umap_harmony_cData') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'status',
              reduction = 'umap_harmony_cData') +
     labs(title = 'status') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'patient',
              reduction = 'umap_harmony_cData') +
     labs(title = 'patient') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16: harmony merge'))
grid.arrange(p)

35.png 36.png

[ ]:

options(repr.plot.width=10, repr.plot.height=20, repr.plot.res = 200)
FeaturePlot(object = seu, features = c('CD68', 'HLA-DRA',
                     'CD14', 'FCGR3A','CLEC9A', 'XCR1','LILRA4','TCF4'),
            reduction = 'umap_harmony_cData', ncol = 2,
            pt.size = 0.01)

37.png 38.png

[235]:

library(future)
plan("multiprocess", workers = 12)
options(future.globals.maxSize = 10000 * 1024^2)
markers <- FindMyMarkers(seu = seu, group.by = 'my.anno')
plan("sequential")

[ ]:

options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
top10 <- extractTopNMarkers(df = markers, top_n = 20, group.order = 0:11)
p <- DoHeatmap(object = seu, draw.lines = FALSE, group.by = 'my.anno',
          features = top10$gene, size = 5, angle = 0,
               disp.min = -2, disp.max = 2) +
theme(text = element_text(size = 4))
p

39.png

[240]:

write.csv(markers, file.path('/home/public/ghx/wuLab/Bladder/0109_Mye_annotation',
                             'markers_Mye_cluster_cData.csv'))
saveRDS(seu, file = file.path('/home/public/ghx/wuLab/Bladder/0109_Mye_annotation',
                              'Mye_cData_seu.rds'))

[ ]:

Bind harmony corrected PCs of Macs with non-corrected PCs of Neus and DCs#

[162]:

bind.embedding <- rbind(seu.sub@reductions$harmony@cell.embeddings,
                        seu@reductions$pca@cell.embeddings[setdiff(Cells(seu),Cells(seu.sub)), ])
bind.embedding <- bind.embedding[Cells(seu), ]
bind.reduc <- CreateDimReducObject(embeddings = bind.embedding,
                                   key = 'bind_pca',
                                   assay = 'RNA')
seu[['bind_pca']] <- bind.reduc

Warning message:
“All keys should be one or more alphanumeric characters followed by an underscore '_', setting key to bind_”

[166]:

seu <- RunUMAP(object = seu, reduction = 'bind_pca',
               reduction.name = 'umap_harmony',
               umap.method = 'umap-learn', metric = 'correlation', dims = 1:20)

[167]:

seu <- FindNeighbors(seu, reduction = 'bind_pca', dims = 1:20, verbose = F)

[168]:

seu <- FindClusters(seu, verbose = F,resolution = 0.3)

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
p1 <- DimPlot(seu, label = T, pt.size = 0.1,
              reduction = 'umap_harmony') +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'orig.ident',
              reduction = 'umap_harmony') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'status',
              reduction = 'umap_harmony') +
     labs(title = 'status') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'patient',
              reduction = 'umap_harmony') +
     labs(title = 'patient') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16: harmony merge'))
grid.arrange(p)

40.png 41.png

put all Mye cells together, annotated with Mac, Neu and DC#

[48]:

seu <- FindClusters(seu, resolution = 1.0, verbose = F)

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
p1 <- DimPlot(seu, label = T, pt.size = 0.1,
              reduction = 'umap_harmony') +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'orig.ident',
              reduction = 'umap_harmony') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'status',
              reduction = 'umap_harmony') +
     labs(title = 'status') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'patient',
              reduction = 'umap_harmony') +
     labs(title = 'patient') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('Mye'))
grid.arrange(p)

42.png 43.png

[ ]:

options(repr.plot.width=20, repr.plot.height=10, repr.plot.res = 200)
p1 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'Mye.ct',
              reduction = 'umap_harmony') +
     labs(title = 'cell type') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'sub.ct',
              reduction = 'umap_harmony') +
     labs(title = 'sub cell type') + theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2), ncol = 2,
                top = paste0('Mye'))
grid.arrange(p)

44.png

[50]:

plan("multiprocess", workers = 8)
options(future.globals.maxSize = 10000 * 1024^2)
markers <- FindAllMarkers(seu, logfc.threshold = 0.1, verbose = F)
plan("sequential")

[51]:

write.csv(markers, file.path('/home/public/ghx/wuLab/Bladder/0508_Mye',
                             'markers_Mye_cluster_res1.0.csv'))

[67]:

ct.annotation.txt <- 'Cluster 0: neu – CMTM2+SELL+
Cluster 1: mac/mono – EREG+inflammasome+
Cluster 2: neu – IL1RN+
Cluster 3: mac/mono – APOE+SPP1+
Cluster 4: DC – CD1C+CLEC10A+
Cluster 5: mac/mono – IL10+HS_ signaling+
Cluster 6: neu – IFN_signaling+
Cluster 7: mac/mono – CXCL10+TIMP1+
Cluster 8: neu – SLPI+MMP9+
Cluster 9: mac/mono – CDKN1C+TCF7L2+
Cluster 10: DC – CD207+CD1E+
Cluster 12: mac/mono – RNASE1+LYVE1+
Cluster 13: mac/mono – CCL4+CLL3+CD14+
Cluster 14: mac/mono – CCL4+CLL3+CD14-
Cluster 15: DC – CLEC9A+WDFY4+
Cluster 16: pDC – TCF4+IL3RA+
Cluster 17: DC – “mregDC”
Cluster 11: remove
Cluster 18: remove
Cluster 19: remove'
ct.annotation.line <- strsplit(ct.annotation.txt, split = '\n')[[1]]
ct.annotation.line <- strsplit(ct.annotation.line, split = ': ', fixed = T)
ct.mapping <- sapply(ct.annotation.line, FUN = function(x){
    return(x[2])
})
c.names <- sapply(ct.annotation.line, FUN = function(x){
    return(x[1])
})
names(ct.mapping) <- substr(c.names, start = 9, stop = nchar(c.names))

[68]:

ct.mapping
0
'neu – CMTM2+SELL+'
1
'mac/mono – EREG+inflammasome+'
2
'neu – IL1RN+'
3
'mac/mono – APOE+SPP1+'
4
'DC – CD1C+CLEC10A+'
5
'mac/mono – IL10+HS_ signaling+'
6
'neu – IFN_signaling+'
7
'mac/mono – CXCL10+TIMP1+'
8
'neu – SLPI+MMP9+'
9
'mac/mono – CDKN1C+TCF7L2+'
10
'DC – CD207+CD1E+'
12
'mac/mono – RNASE1+LYVE1+'
13
'mac/mono – CCL4+CLL3+CD14+'
14
'mac/mono – CCL4+CLL3+CD14-'
15
'DC – CLEC9A+WDFY4+'
16
'pDC – TCF4+IL3RA+'
17
'DC – “mregDC”'
11
'remove'
18
'remove'
19
'remove'
[ ]:

big.cell.type <- as.character(seu@active.ident)
sub.cell.type <- as.character(seu@active.ident)
for(cluster in names(ct.mapping)){
    big.cell.type[big.cell.type %in% cluster] <-
}
cell.type[cell.type %in% c()]

[ ]:



[ ]:



[ ]:



[ ]: