[7]:

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

Integrate#

[8]:

library(permute)

[9]:

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

[10]:

options(warn=-1)

[11]:

# idxs indicates that which patients will be used
idxs <- 5:16
samples <- paste0('MIBC',
                  rep(idxs, each = 2),
                  rep(c('P','T'), times = length(idxs)))
print(samples)

 [1] "MIBC5P"  "MIBC5T"  "MIBC6P"  "MIBC6T"  "MIBC7P"  "MIBC7T"  "MIBC8P"
 [8] "MIBC8T"  "MIBC9P"  "MIBC9T"  "MIBC10P" "MIBC10T" "MIBC11P" "MIBC11T"
[15] "MIBC12P" "MIBC12T" "MIBC13P" "MIBC13T" "MIBC14P" "MIBC14T" "MIBC15P"
[22] "MIBC15T" "MIBC16P" "MIBC16T"

[12]:

overall.height <- 5.0*((length(samples)+1) %/% 2 )
print(overall.height)

[1] 60

load SoupX-corrected raw matrix and annotations generate by “annotation_5-16.ipynb”#

[5]:

matrixPath <- '/home/public/ghx/wuLab/Bladder/1217_all/matrices'
annotationPath <- '/home/public/ghx/wuLab/Bladder/1217_all/seu'
seu.list <- list()
for(sample in samples){
    mtx <- readRDS(file = file.path(matrixPath,
                                    paste0(sample, '_soup_matrix.rds')))
    anno <- readRDS(file.path(annotationPath,
                            paste0(sample,'_annotation.rds')))
    seu.list[[sample]] <- CreateSeuratObject(counts = mtx,
                                             min.cells = 0, min.features = 0,
                                             meta.data = anno, project = sample)
}

remove unknown cells, dead cells, erythrocyte and proliferating cells#

[14]:

for(sample in samples){
    removed.cells <- Cells(seu.list[[sample]])[is.na(seu.list[[sample]]@meta.data$cell.type) |
                                               seu.list[[sample]]@meta.data$dead.cells |
                                               seu.list[[sample]]@meta.data$cell.type %in% c('UK')]
    seu.list[[sample]] <- SubsetData(object = seu.list[[sample]],
                           cells = setdiff(Cells(seu.list[[sample]]),
                                           removed.cells))
}

merge cell matrix and annotations#

[31]:

anno.used <- c('orig.ident','nCount_RNA','nFeature_RNA','seurat_clusters',
'cell.type','db.score','mt.percentage','S.Score',
'G2M.Score','Phase')

[32]:

mtx.all <- as.matrix(matrix(rep(0,33538),ncol = 1), "dgCMatrix")
anno.all <- data.frame()
for(sample in samples){
    mtx.all <- cbind(mtx.all, seu.list[[sample]]@assays$RNA@counts)
    temp.anno <- seu.list[[sample]]@meta.data
    temp.anno <- temp.anno[anno.used]
    temp.anno$orig.ident <- sample
    temp.anno$sample.cluster <- paste0(sample, '-', temp.anno$seurat_clusters)
    anno.all <- rbind(anno.all, temp.anno)
}

[33]:

mtx.all <- mtx.all[,-1]

[35]:

dim(mtx.all)
dim(anno.all)
  1. 33538
  2. 246754
  1. 246754
  2. 11
[36]:

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

[37]:

seu.all

An object of class Seurat
33538 features across 246754 samples within 1 assay
Active assay: RNA (33538 features)

[38]:

seu.all@meta.data$orig.ident <- factor(seu.all@meta.data$orig.ident,
                                       levels = samples)

[45]:

rm(seu.list)
gc()
A matrix: 2 × 6 of type dbl
used(Mb)gc trigger(Mb)max used(Mb)
Ncells 2408413 128.7 4431859 236.7 4431859 236.7
Vcells7090047805409.3433530036633075.8541912545741344.7
 
[46]:

rm(mtx.all)
gc()
A matrix: 2 × 6 of type dbl
used(Mb)gc trigger(Mb)max used(Mb)
Ncells 2408942 128.7 4431859 236.7 4431859 236.7
Vcells7090058855409.3346824029326460.6541912545741344.7
 
[48]:

table(seu.all@meta.data$cell.type)


     T     NK      B    Pla    Mye   Mast    Fib    SMC    Epi    Umb   Endo
120535  13106   3468   6555  26329  22099  10018   3000  35507   2832   3305
    UK
     0

[1]:

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

Error in Nml_Scl_Fvg_PCA_Fnb(seu = seu.all): could not find function "Nml_Scl_Fvg_PCA_Fnb"
Traceback:


[ ]:

seu.all@active.ident <- factor(seu.all@active.ident,
                          levels = 0:(length(unique(seu.all@active.ident))-1))

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
p1 <- DimPlot(seu.all, label = T, pt.size = 0.1) +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu.all, label = F, pt.size = 0.1, group.by = 'orig.ident') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu.all, label = T, pt.size = 0.1, group.by = 'cell.type') +
     labs(title = 'cell.type') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu.all, label = T, pt.size = 0.1, group.by = 'status') +
     labs(title = 'status') #+ 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/graph/overview',
                            'Merge_raw_overview.jpg'))
grid.arrange(p)

1.png

[20]:

ggsave(plot = p1, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_raw_cluster.jpg'))
ggsave(plot = p2, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_raw_origin.jpg'))
ggsave(plot = p3, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_raw_cellType.jpg'))
ggsave(plot = p4, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_raw_status.jpg'))

[21]:

plots <- list()
for(sample in samples){
    p <- DimPlot(seu.all, label = T, pt.size = 0.01, sizes.highlight = 0.2,
                 cells.highlight = Cells(seu.all)[seu.all@meta.data$orig.ident == sample]) +
         labs(title = sample) #+ theme(legend.position = 'none')
    plots[[sample]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                top = paste0('MIBC5-16: sampleOrigin'))
ggsave(plot = p, limitsize = FALSE,
       width = 10, height = overall.height,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_raw_highlight_sampleOrigin.jpg'))

[22]:

plots <- list()
for(ct in levels(seu.all@meta.data$cell.type)){
    temp.cells <- Cells(seu.all)[seu.all@meta.data$cell.type == ct]
    p <- DimPlot(seu.all, label = T, pt.size = 0.1, sizes.highlight = 0.2,
                 cells.highlight = temp.cells) +
         labs(title = ct) #+ theme(legend.position = 'none')
    plots[[ct]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                top = paste0('MIBC5-16: cell type'))
ggsave(plot = p, limitsize = FALSE,
       width = 10, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_raw_highlight_cellType.jpg'))

[39]:

plots <- list()
for(s in c('P','T')){
    temp.cells <- Cells(seu.all)[seu.all@meta.data$status == s]
    p <- DimPlot(seu.all, label = T, pt.size = 0.01, sizes.highlight = 0.2,
                 cells.highlight = temp.cells) +
         labs(title = s) #+ theme(legend.position = 'none')
    plots[[s]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                top = paste0('MIBC5-16: status'))
ggsave(plot = p, limitsize = FALSE,
       width = 10, height = 5,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_raw_highlight_status.jpg'))

[63]:

saveRDS(seu.all, file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                  'Merge_raw_seu.rds'))

[80]:

saveRDS(seu.all@assays$RNA@counts, file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                  'Merge_raw_mtx.rds'))
saveRDS(seu.all@meta.data, file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                  'Merge_raw_annotation.rds'))

[82]:

for(ct in levels(seu.all@meta.data$cell.type)){
    temp.cells <- Cells(seu.all)[seu.all@meta.data$cell.type == ct]
    temp.mtx <- seu.all@assays$RNA@counts[, temp.cells]
    temp.anno <- seu.all@meta.data[temp.cells, ]
    saveRDS(temp.mtx,
            file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                            paste0('Celltype_',ct,'_matrix.rds')))
    saveRDS(temp.anno,
            file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                            paste0('Celltype_',ct,'_annotation.rds')))
}

remove doublets by cell types#

find cell type markers genes#

[71]:

cell.subset <- c()
for(ct in levels(seu.all@meta.data$cell.type)){
    temp.cells <- Cells(seu.all)[seu.all@meta.data$cell.type == ct]
    temp.sub <- sample(temp.cells,
                       size = min(2000, length(temp.cells)),
                       replace = F)
    cell.subset <- c(cell.subset,temp.sub)
}

[72]:

seu.sub <- CreateSeuratObject(counts = seu.all@assays$RNA@counts[, cell.subset],
                              project = 'sub')
seu.sub <- NormalizeData(seu.sub, verbose = F)

[75]:

seu.sub@meta.data$cell.type <- seu.all@meta.data[ cell.subset, ]$cell.type

[73]:

library(future)

[74]:

plan("multiprocess", workers = 12)
plan()
options(future.globals.maxSize = 10000 * 1024^2)

structure(function (expr, envir = parent.frame(), substitute = TRUE,
    lazy = FALSE, seed = NULL, globals = TRUE, workers = 12, 
    gc = FALSE, earlySignal = FALSE, label = NULL, ...) 
{
    if (substitute) 
        expr <- substitute(expr)
    fun <- if (supportsMulticore(warn = TRUE)) 
        multicore
    else multisession
    fun(expr = expr, envir = envir, substitute = FALSE, lazy = lazy, 
        seed = seed, globals = globals, workers = workers, gc = gc, 
        earlySignal = earlySignal, label = label, ...)
}, class = c("FutureStrategy", "tweaked", "multiprocess", "future",
"function"), init = FALSE, call = plan("multiprocess", workers = 12))
[76]:

ct.markers <- FindMyMarkers(seu = seu.sub,
                            group.by = 'cell.type')

[77]:

plan("sequential")

[ ]:

options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
tops <- extractTopNMarkers(df = ct.markers, top_n = 20)
p <- DoHeatmap(object = seu.all,
               group.by = 'cell.type', cells = cell.subset,
               features = tops$gene, draw.lines = FALSE,
               size = 5, angle = 0, disp.min = -2, disp.max = 2) +
        theme(text = element_text(size = 6)) +
        labs(title = 'Merge_raw')
ggsave(p, width = 10, height = 10,limitsize = FALSE,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                             'Merge_raw_heatmap.jpg'))
p

2.png

[79]:

write.csv(ct.markers,
          file.path('/home/public/ghx/wuLab/Bladder/1217_all/markers',
                    'Merge_raw_ct_markers.csv'))

find doublets pipeline is in ‘preprocessing/doublets_by_celltype/db_by_ct_1217_all.ipynb’#

[6]:

seu.all <- readRDS(file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                    'Merge_raw_seu.rds'))
seu.all@meta.data$status <- sapply(as.character(seu.all@meta.data$orig.ident), USE.NAMES = F, FUN = function(x){
    a <- substring(x, first = nchar(x), last = nchar(x))
    return(a)
})

[7]:

new.anno <- readRDS(file.path('/home/public/ghx/wuLab/Bladder/1217_all/doublet',
                              'new_annotation.rds'))
seu.all@meta.data$cell.type.v1 <- new.anno[Cells(seu.all)]

[8]:

table(new.anno)

new.anno
      T      NK       B     Pla     Mye    Mast     Fib     SMC     Epi     Umb
 100492    9855    2520    5521   21365   19345    7119    3263   27704    2014
   Endo      UK doublet
   2833    3018   41705

[9]:

plots <- list()
for(ct in levels(seu.all@meta.data$cell.type.v1)){
    # p1: doublets or cell cycle cells in ct cells
    temp.cells.1 <- Cells(seu.all)[seu.all@meta.data$cell.type == ct &
                                   seu.all@meta.data$cell.type.v1 %in% c('UK','doublet')]
    p1 <- DimPlot(seu.all, label = T, pt.size = 0.1, sizes.highlight = 0.2,
                 cells.highlight = temp.cells.1) +
         labs(title = paste0(ct,'_doublet&UK')) #+ theme(legend.position = 'none')
    # p2: clean cells in ct cells
    temp.cells.2 <- Cells(seu.all)[seu.all@meta.data$cell.type.v1 == ct]
    p2 <- DimPlot(seu.all, label = T, pt.size = 0.1, sizes.highlight = 0.2,
                 cells.highlight = temp.cells.2) +
         labs(title = paste0(ct,'_clean')) #+ theme(legend.position = 'none')
    plots[[paste(ct,'1')]] <- p1
    plots[[paste(ct,'2')]] <- p2
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                top = paste0('MIBC5-16: cell type'))
ggsave(plot = p, limitsize = FALSE,
       width = 10, height = 60,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_highlight_doublets.jpg'))

remove doublets and cell cycle cells, then do re-preprocessing allover again#

[13]:

new.anno <- readRDS(file.path('/home/public/ghx/wuLab/Bladder/1217_all/doublet',
                              'new_annotation.rds'))

[14]:

cells.keep <- names(new.anno)[ !new.anno %in% c('UK','doublet')]

[15]:

anno <- readRDS(file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                 'Merge_raw_annotation.rds'))
mtx <- readRDS(file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                'Merge_raw_mtx.rds'))
anno <- anno[cells.keep, ]
mtx <- mtx[,cells.keep]

[16]:

seu.all <- CreateSeuratObject(counts = mtx,
                              meta.data = anno,
                              project = 'Merge_v1')

[23]:

seu.all@meta.data$cell.type.v1 <- new.anno[Cells(seu.all)]

[17]:

seu.all@meta.data$status <- sapply(as.character(seu.all@meta.data$orig.ident),
                                   USE.NAMES = F, FUN = function(x){
    a <- substring(x, first = nchar(x), last = nchar(x))
    return(a)
})
seu.all@meta.data$patient <- sapply(X = as.character(seu.all@meta.data$orig.ident),
                  USE.NAMES = F, FUN = function(x){
                      a <- substring(text = x, first = 1, last = nchar(x)-1)
                  })

[18]:

table(seu.all@meta.data$orig.ident)


 MIBC5P  MIBC5T  MIBC6P  MIBC6T  MIBC7P  MIBC7T  MIBC8P  MIBC8T  MIBC9P  MIBC9T
   6148    6266   13536   14178   10890   10800   16693    3460    7026    7130
MIBC10P MIBC10T MIBC11P MIBC11T MIBC12P MIBC12T MIBC13P MIBC13T MIBC14P MIBC14T
   7566    7267   14161   13672   12947   12575    3842    6078    4938    1346
MIBC15P MIBC15T MIBC16P MIBC16T
   3966    7629    6756    3161

[19]:

seu.all <- Nml_Scl_Fvg_PCA_Fnb(seu = seu.all)
seu.all <- RunUMAP(object = seu.all, min.dist = 0.7,
                   dims = 1:25, verbose = F,
               umap.method = 'umap-learn', metric = 'correlation')
seu.all <- FindClusters(object = seu.all, resolution = 0.25, verbose = F)

[25]:

stats <- table(seu.all@meta.data$orig.ident, seu.all@meta.data$cell.type.v1)
stats <- as.data.frame(stats[, -ncol(stats)])
names(stats) <- c('origin', 'cell.type', 'cell.number')
p5 <- ggplot(stats, aes(fill=cell.type, y=cell.number, x=origin)) +
        geom_bar( stat="identity", position="fill", color = 'black', ) +
        theme(axis.text.x = element_text(size=10, angle = 90))

[26]:

ggsave(plot = p1, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_cluster.jpg'))
ggsave(plot = p2, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_origin.jpg'))
ggsave(plot = p3, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_cellType.jpg'))
ggsave(plot = p4, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_status.jpg'))
ggsave(plot = p5, limitsize = FALSE, width = 10, height = 10,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_composition.jpg'))

[31]:

plots <- list()
for(sample in samples){
    p <- DimPlot(seu.all, label = T, pt.size = 0.01, sizes.highlight = 0.1,
                 cells.highlight = Cells(seu.all)[seu.all@meta.data$orig.ident == sample]) +
         labs(title = sample) + theme(legend.position = 'none')
    plots[[sample]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                top = paste0('MIBC5-16: sampleOrigin'))
ggsave(plot = p, limitsize = FALSE,
       width = 10, height = overall.height,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_highlight_sampleOrigin.jpg'))

[32]:

plots <- list()
for(ct in levels(seu.all@meta.data$cell.type.v1)){
    if(ct %in% c('UK','doublet')) {next}
    temp.cells <- Cells(seu.all)[seu.all@meta.data$cell.type.v1 == ct]
    p <- DimPlot(seu.all, label = T, pt.size = 0.1, sizes.highlight = 0.1,
                 cells.highlight = temp.cells) +
         labs(title = ct) + theme(legend.position = 'none')
    plots[[ct]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                top = paste0('MIBC5-16: cell type'))
ggsave(plot = p, limitsize = FALSE,
       width = 20, height = 60,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_highlight_cellType.jpg'))

[33]:

plots <- list()
for(s in c('P','T')){
    temp.cells <- Cells(seu.all)[seu.all@meta.data$status == s]
    p <- DimPlot(seu.all, label = T, pt.size = 0.01, sizes.highlight = 0.1,
                 cells.highlight = temp.cells) +
         labs(title = s) + theme(legend.position = 'none')
    plots[[s]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                top = paste0('MIBC5-16: status'))
ggsave(plot = p, limitsize = FALSE,
       width = 30, height = 15,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_highlight_status.jpg'))

[29]:

plots <- list()
for(sample in samples){
    temp.cells <- Cells(seu.all)[seu.all@meta.data$orig.ident == sample]
    other.cells <- setdiff(Cells(seu.all),temp.cells)
    seu.all@meta.data$temp.cluster <- seu.all@meta.data$cell.type.v1
    seu.all@meta.data[other.cells,]$temp.cluster <- NA
    p <- DimPlot(seu.all, group.by = 'temp.cluster',
                 label = T, pt.size = 0.01,
                 cells = c(other.cells,temp.cells)) +
         labs(title = sample) + theme(legend.position = 'none')
    plots[[sample]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                top = paste0('MIBC5-16: sampleOrigin'))
ggsave(plot = p, limitsize = FALSE,
       width = 10, height = overall.height,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_split_sampleOrigin.jpg'))

[30]:

plots <- list()
for(s in c('P','T')){
    temp.cells <- Cells(seu.all)[seu.all@meta.data$status == s]
    other.cells <- setdiff(Cells(seu.all),temp.cells)
    seu.all@meta.data$temp.cluster <- seu.all@meta.data$cell.type.v1
    seu.all@meta.data[other.cells,]$temp.cluster <- NA
    p <- DimPlot(seu.all, group.by = 'temp.cluster',
                 label = T, pt.size = 0.01,
                 cells = c(other.cells,temp.cells)) +
         labs(title = s) + theme(legend.position = 'none')
    plots[[s]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                top = paste0('MIBC5-16: status'))
ggsave(plot = p, limitsize = FALSE,
       width = 30, height = 15,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_split_status.jpg'))

batch effect removal: only remove difference between patients, but not Tumor/Paratumor#

[27]:

library(harmony)

Loading required package: Rcpp


[32]:

seu.all <- RunHarmony(seu.all, 'patient', theta = 10, max.iter.harmony = 20,
                      plot_convergence = F,verbose = F)

[33]:

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

[34]:

p1 <- DimPlot(seu.all, label = T, pt.size = 0.1,
              reduction = 'umap_harmony') +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu.all, label = F, pt.size = 0.1, group.by = 'orig.ident',
              reduction = 'umap_harmony') +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu.all, label = T, pt.size = 0.1, group.by = 'cell.type.v1',
              reduction = 'umap_harmony') +
     labs(title = 'cell.type') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu.all, label = T, pt.size = 0.1, group.by = 'status',
              reduction = 'umap_harmony') +
     labs(title = 'status') #+ theme(legend.position = 'none')

p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16:merge_v1 harmony'))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Merge_v1_overview_harmony.jpg'))

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
grid.arrange(p)

3.png 4.png

[34]:

saveRDS(seu.all, file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                  'Merge_v1_seu.rds'))

[35]:

saveRDS(seu.all@assays$RNA@counts, file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                  'Merge_v1_mtx.rds'))
saveRDS(seu.all@meta.data, file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                  'Merge_v1_annotation.rds'))

[22]:

for(ct in levels(seu.all@meta.data$cell.type)){
    if(ct == 'T') next
    temp.cells <- Cells(seu.all)[seu.all@meta.data$cell.type.v1 == ct]
    temp.mtx <- seu.all@assays$RNA@counts[, temp.cells]
    temp.anno <- seu.all@meta.data[temp.cells, ]
    saveRDS(temp.mtx,
            file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                            paste0('V1_Celltype_',ct,'_matrix.rds')))
    saveRDS(temp.anno,
            file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                            paste0('V1_Celltype_',ct,'_annotation.rds')))
}

test sample difference by celltype#

[74]:

plots <- list()
var.genes <- seu.all@assays$RNA@var.features
for(ct in levels(seu.all@meta.data$cell.type.v1)){
    mean.expression <- data.frame()
    sample.used <- c()
    for(sample in samples){
        temp.cells <- Cells(seu.all)[seu.all@meta.data$orig.ident == sample &
                                     seu.all@meta.data$cell.type.v1 == ct]
        if(length(temp.cells)<=1){
            next
        }
        temp.mean <- rowMeans(seu.all@assays$RNA@data[var.genes, temp.cells])
        if(!is.na(temp.mean[1])){
            mean.expression <- rbind(mean.expression,temp.mean)
            sample.used <- c(sample.used, sample)
            row.names(mean.expression) <- sample.used
        }

    }
    colnames(mean.expression) <- var.genes
    mean.expression <- as.matrix(mean.expression)
    mean.pca <- prcomp(x = mean.expression)

    plot.data <- as.data.frame(mean.pca$x)
    plot.data$sample <- row.names(plot.data)
    plot.data$condition <- substr(row.names(plot.data),
                                  start = nchar(row.names(plot.data)),
                                  stop = nchar(row.names(plot.data)))

    p <- ggplot(data = plot.data, aes(x=PC1, y=PC2, colour=condition)) +
    geom_point(size=5,shape=16) +
    geom_text(aes(label = sample), nudge_y = -0.25) +
    labs(title = ct)

    plots[[ct]] <- p
}

Error in names(x) <- value: 'names' attribute [1000] must be the same length as the vector [0]
Traceback:

1. `colnames<-`(`*tmp*`, value = var.genes)

[ ]:

options(repr.plot.width=10, repr.plot.height=30, repr.plot.res = 200)
grid.arrange(grobs = plots, ncol = 2, width = c(1,1), heigths = c(0.5,0.5),
            top = 'pca by celltype')

5.png 6.png i7mage.png 8.png 9.png 10.png

[1]:

anno <- readRDS(file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                  'Merge_v1_annotation.rds'))

[2]:

library(pheatmap)

[ ]:

stats <- table(anno$cell.type.v1,
               anno$orig.ident)
stats <- stats[rowSums(stats)>0,colSums(stats)>0]
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=unique(patient)[order(as.integer(unique(patient)))])
status <- substr(samples, start = sapply(samples,FUN=nchar),
                  stop = sapply(samples,FUN=nchar))
status <- factor(status)
row.anno <- data.frame(row.names = colnames(stats),
                       patient = patient, status = status)
pheatmap(t(stats),
         cutree_rows = min(6,nrow(stats)),
         cutree_cols = min(5,ncol(stats)),
         annotation_row = row.anno)

11.png

[ ]:

CD45+ cells#

[7]:

new.anno <- readRDS(file.path('/home/public/ghx/wuLab/Bladder/1217_all/doublet',
                              'new_annotation.rds'))

[8]:

cells.keep <- names(new.anno)[ !new.anno %in% c('UK','doublet')]

[9]:

anno <- readRDS(file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                 'Merge_raw_annotation.rds'))
mtx <- readRDS(file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                'Merge_raw_mtx.rds'))
anno <- anno[cells.keep, ]
mtx <- mtx[,cells.keep]

[11]:

anno$cell.type.v1 <- new.anno[row.names(anno)]
cells.use <- row.names(anno)[ anno$cell.type.v1 %in% c('T','NK','B','Pla','Mye','Mast')]

[13]:

anno <- anno[cells.use,]
mtx <- mtx[, cells.use]

[14]:

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

[15]:

seu@meta.data$status <- sapply(as.character(seu@meta.data$orig.ident),
                                   USE.NAMES = F, FUN = function(x){
    a <- substring(x, first = nchar(x), last = nchar(x))
    return(a)
})
seu@meta.data$patient <- sapply(X = as.character(seu@meta.data$orig.ident),
                  USE.NAMES = F, FUN = function(x){
                      a <- substring(text = x, first = 1, last = nchar(x)-1)
                  })

[16]:

table(seu@meta.data$orig.ident)


 MIBC5P  MIBC5T  MIBC6P  MIBC6T  MIBC7P  MIBC7T  MIBC8P  MIBC8T  MIBC9P  MIBC9T
   4557    5281   12142   13683    9841    8131   15625    2987    3958    6131
MIBC10P MIBC10T MIBC11P MIBC11T MIBC12P MIBC12T MIBC13P MIBC13T MIBC14P MIBC14T
   6196    6871    7048    4486    8701   12015    3762    5909    4196    1073
MIBC15P MIBC15T MIBC16P MIBC16T
   3120    5939    4330    3116

[17]:

seu <- Nml_Scl_Fvg_PCA_Fnb(seu = seu)
seu <- RunUMAP(object = seu, min.dist = 0.7,
                   dims = 1:25, verbose = F,
               umap.method = 'umap-learn', metric = 'correlation')
seu <- FindClusters(object = seu, resolution = 0.25, verbose = F)

[19]:

seu@active.ident <- factor(seu@active.ident,
                          levels = 0:(length(unique(seu@active.ident))-1))
seu@meta.data$cell.type.v1 <- factor(seu@meta.data$cell.type.v1,
                                     levels = c('T','NK','B','Pla','Mye','Mast'))

[ ]:



[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
cells.reorder <- Cells(seu)[ shuffle(length(Cells(seu)))]
p1 <- DimPlot(seu, label = T, pt.size = 0.1,
              cells = cells.reorder) +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = F, pt.size = 0.1,
              group.by = 'orig.ident', cells = cells.reorder) +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu, label = T, pt.size = 0.1,
              group.by = 'cell.type.v1', cells = cells.reorder) +
     labs(title = 'cell.type') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu, label = T, pt.size = 0.1,
              group.by = 'status', cells = cells.reorder) +
     labs(title = 'status') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16:immune_v1'))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Immune_v1_overview.jpg'))
grid.arrange(p)

12.png 13.png

[27]:

library(harmony)

Loading required package: Rcpp


[ ]:

options(repr.plot.width=10, repr.plot.height=5, repr.plot.res = 200)
seu <- RunHarmony(seu, theta = 10, verbose = F,
                  group.by.vars = 'patient', plot_convergence = T)

14.png

[51]:

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

[30]:

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

[31]:

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

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
cells.reorder <- Cells(seu)[ !seu@active.ident %in% c(12,14,15,16,20,21)]
cells.reorder <- cells.reorder[ shuffle(length(cells.reorder))]
p1 <- DimPlot(seu, label = T, pt.size = 0.1,
              reduction = 'umap_harmony', cells = cells.reorder) +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'orig.ident',
              reduction = 'umap_harmony', cells = cells.reorder) +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'cell.type.v1',
              reduction = 'umap_harmony', cells = cells.reorder) +
     labs(title = 'status') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'status',
              reduction = 'umap_harmony', cells = cells.reorder) +
     labs(title = 'patient') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16: immune harmony'))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Immune_v1_harmony_overview.jpg'))
grid.arrange(p)

15.png 16.png

[47]:

stats <- table(seu@meta.data$cell.type.v1, seu@meta.data$orig.ident)
nb.cluster <- nrow(stats)
stats <- as.data.frame(stats)
names(stats) <- c('cluster', 'origin', 'cell.number')
color <- rep(c(rep('white',nb.cluster),
               rep('black',nb.cluster)), 12)
p5 <- ggplot(stats, aes(fill=cluster, y=cell.number, x=origin)) +
        geom_bar( stat="identity", position="fill", color = color) +
        theme(axis.text.x = element_text(size=10, angle = 90))

[53]:

ggsave(plot = p1, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Immune_v1_cluster.jpg'))
ggsave(plot = p2, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Immune_v1_origin.jpg'))
ggsave(plot = p3, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Immune_v1_cellType.jpg'))
ggsave(plot = p4, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Immune_v1_status.jpg'))
ggsave(plot = p5, limitsize = FALSE, width = 10, height = 10,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'Immune_v1_composition.jpg'))

[54]:

saveRDS(seu, file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                       'Immune_v1_seu.rds'))

CD45- cells#

[55]:

new.anno <- readRDS(file.path('/home/public/ghx/wuLab/Bladder/1217_all/doublet',
                              'new_annotation.rds'))

[56]:

cells.keep <- names(new.anno)[ !new.anno %in% c('UK','doublet')]

[57]:

anno <- readRDS(file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                 'Merge_raw_annotation.rds'))
mtx <- readRDS(file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                                'Merge_raw_mtx.rds'))
anno <- anno[cells.keep, ]
mtx <- mtx[,cells.keep]

[58]:

anno$cell.type.v1 <- new.anno[row.names(anno)]
cells.use <- row.names(anno)[ anno$cell.type.v1 %in% c('SMC','Fib','Epi','Umb','Endo')]

[59]:

anno <- anno[cells.use,]
mtx <- mtx[, cells.use]

[60]:

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

[61]:

seu@meta.data$status <- sapply(as.character(seu@meta.data$orig.ident),
                                   USE.NAMES = F, FUN = function(x){
    a <- substring(x, first = nchar(x), last = nchar(x))
    return(a)
})
seu@meta.data$patient <- sapply(X = as.character(seu@meta.data$orig.ident),
                  USE.NAMES = F, FUN = function(x){
                      a <- substring(text = x, first = 1, last = nchar(x)-1)
                  })

[ ]:

seu <- Nml_Scl_Fvg_PCA_Fnb(seu = seu)
seu <- RunUMAP(object = seu, min.dist = 0.7,
                   dims = 1:25, verbose = F,
               umap.method = 'umap-learn', metric = 'correlation')
seu <- FindClusters(object = seu, resolution = 0.25, verbose = F)

[ ]:

seu@active.ident <- factor(seu@active.ident,
                          levels = 0:(length(unique(seu@active.ident))-1))
seu@meta.data$cell.type.v1 <- factor(seu@meta.data$cell.type.v1,
                                     levels = c('SMC','Fib','Epi','Umb','Endo'))

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
cells.reorder <- Cells(seu)[ shuffle(length(Cells(seu)))]
p1 <- DimPlot(seu, label = T, pt.size = 0.1,
              cells = cells.reorder) +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = F, pt.size = 0.1,
              group.by = 'orig.ident', cells = cells.reorder) +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu, label = T, pt.size = 0.1,
              group.by = 'cell.type.v1', cells = cells.reorder) +
     labs(title = 'cell.type') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu, label = T, pt.size = 0.1,
              group.by = 'status', cells = cells.reorder) +
     labs(title = 'status') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16:nonimmune_v1'))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'nonImmune_v1_overview.jpg'))
grid.arrange(p)

[ ]:

options(repr.plot.width=5, repr.plot.height=10, repr.plot.res = 200)
seu <- RunHarmony(seu, theta = 10, verbose = F,
                  group.by.vars = 'patient', plot_convergence = T)

[ ]:

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

[ ]:

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

[ ]:

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

[ ]:

options(repr.plot.width=20, repr.plot.height=20, repr.plot.res = 200)
cells.reorder <- Cells(seu)[ shuffle(length(Cells(seu)))]
p1 <- DimPlot(seu, label = T, pt.size = 0.1,
              reduction = 'umap_harmony', cells = cells.reorder) +
     labs(title = 'cluster') #+ theme(legend.position = 'none')
p2 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'orig.ident',
              reduction = 'umap_harmony', cells = cells.reorder) +
     labs(title = 'sample.origin') #+ theme(legend.position = 'none')
p3 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'cell.type.v1',
              reduction = 'umap_harmony', cells = cells.reorder) +
     labs(title = 'status') #+ theme(legend.position = 'none')
p4 <- DimPlot(seu, label = T, pt.size = 0.1, group.by = 'status',
              reduction = 'umap_harmony', cells = cells.reorder) +
     labs(title = 'patient') #+ theme(legend.position = 'none')
p <- arrangeGrob(grobs = list(p1,p2,p3,p4), ncol = 2,
                top = paste0('MIBC5-16: nonImmune harmony'))
ggsave(plot = p, limitsize = FALSE, width = 20, height = 20,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'nonImmune_v1_harmony_overview.jpg'))
grid.arrange(p)

[ ]:

stats <- table(seu@meta.data$cell.type.v1, seu@meta.data$orig.ident)
nb.cluster <- nrow(stats)
stats <- as.data.frame(stats)
names(stats) <- c('cluster', 'origin', 'cell.number')
color <- rep(c(rep('white',nb.cluster),
               rep('black',nb.cluster)), 12)
p5 <- ggplot(stats, aes(fill=cluster, y=cell.number, x=origin)) +
        geom_bar( stat="identity", position="fill", color = color) +
        theme(axis.text.x = element_text(size=10, angle = 90))

[ ]:

ggsave(plot = p1, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'nonImmune_v1_cluster.jpg'))
ggsave(plot = p2, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'nonImmune_v1_origin.jpg'))
ggsave(plot = p3, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'nonImmune_v1_cellType.jpg'))
ggsave(plot = p4, limitsize = FALSE, width = 30, height = 30,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'nonImmune_v1_status.jpg'))
ggsave(plot = p5, limitsize = FALSE, width = 10, height = 10,
       filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/overview',
                            'nonImmune_v1_composition.jpg'))

[ ]:

saveRDS(seu, file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                       'nonImmune_v1_seu.rds'))

[ ]: