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

Annotation

library(DoubletFinder)

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

patient #5 ~ #16 are used

# 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"

#the overall height is set if each sample generate
# a square graph and graphs are layout as ncol=2
# and each graph has a height of 5
overall.height <- 5.0*((length(samples)+1) %/% 2 )
print(overall.height)

[1] 60

read in 10x raw matrix and draw VlnPlots

mtx.path <- '/home/public/ghx/wuLab/Bladder/MIBC'
mtx.list <- list()
for(sample in samples){
    temp.path <- file.path(mtx.path, sample)
    temp.mtx <- Read10X(temp.path)
    # there will be some infinity and a lot of sequencing error on barcodes
    temp.mtx <- temp.mtx[, colSums(temp.mtx) >= 10 &
                           colSums(temp.mtx) <= 500000]
    # add prefix to cell names
    temp.cells <- colnames(temp.mtx)
    temp.cells <- paste0(sample,'.',temp.cells)
    colnames(temp.mtx) <- temp.cells
    mtx.list[[sample]] <- temp.mtx
}

plots <- list()
for(sample in samples){

    UMI.counts <- data.frame(UMI.count = colSums(mtx.list[[sample]]),
                             cluster = 'whole')
    p <- ggplot(UMI.counts, aes(x = cluster, y = UMI.count, fill = cluster)) +
          geom_violin(trim = TRUE) + scale_y_log10() +
          geom_hline(aes(yintercept=1000)) + geom_hline(aes(yintercept=70000))+
          geom_jitter(size=0.01) +
          labs(title=sample) + theme(legend.position = 'none')
    plots[[sample]] <- p

}
p <- arrangeGrob(grobs = plots, ncol = 2,
                 top = 'nUMI')
ggsave(p, filename = '/home/public/ghx/wuLab/Bladder/1217_all/nUMI.jpg',
       limitsize = FALSE,
       width = 10, height = overall.height)

setting cutoff of cells and backgrounds for each sample

# cut.off.min is the minumum gene count for a barcode to be recognized as a cell
# cut.off.max is the max term of above
# cut.off.max.bk is the maximum for background
cut.off.max.bk <- c('MIBC5P' = 1000,'MIBC5T' = 1000,
                    'MIBC6P' = 500,'MIBC6T' = 800,
                    'MIBC7P' = 300,'MIBC7T' = 300,
                    'MIBC8P' = 500,'MIBC8T' = 500,
                    'MIBC9P' = 1000,'MIBC9T' = 1000,
                    'MIBC10P' = 500,'MIBC10T' = 500,
                    'MIBC11P' = 300,'MIBC11T' = 1100,
                    'MIBC12P' = 1000,'MIBC12T' = 500,
                    'MIBC13P' = 100,'MIBC13T' = 300,
                    'MIBC14P' = 300,'MIBC14T' = 1000,
                    'MIBC15P' = 300,'MIBC15T' = 300,
                    'MIBC16P' = 800,'MIBC16T' = 800)

cut.off.min <- c('MIBC5P' = 2000,'MIBC5T' = 1200,
                'MIBC6P' = 800,'MIBC6T' = 800,
                'MIBC7P' = 300,'MIBC7T' = 800,
                'MIBC8P' = 800,'MIBC8T' = 800,
                'MIBC9P' = 1000,'MIBC9T' = 1000,
                'MIBC10P' = 700,'MIBC10T' = 700,
                'MIBC11P' = 500,'MIBC11T' = 1200,
                'MIBC12P' = 1000,'MIBC12T' = 700,
                'MIBC13P' = 300,'MIBC13T' = 1000,
                'MIBC14P' = 1000,'MIBC14T' = 1000,
                'MIBC15P' = 1000,'MIBC15T' = 1000,
                'MIBC16P' = 1000,'MIBC16T' = 1000)

cut.off.max <- c('MIBC5P' = 100000,'MIBC5T' = 100000,
                'MIBC6P' = 30000,'MIBC6T' = 70000,
                'MIBC7P' = 30000,'MIBC7T' = 100000,
                'MIBC8P' = 30000,'MIBC8T' = 30000,
                'MIBC9P' = 70000,'MIBC9T' = 50000,
                'MIBC10P' = 70000,'MIBC10T' = 40000,
                'MIBC11P' = 70000,'MIBC11T' = 50000,
                'MIBC12P' = 100000,'MIBC12T' = 70000,
                'MIBC13P' = 30000,'MIBC13T' = 70000,
                'MIBC14P' = 70000,'MIBC14T' = 70000,
                'MIBC15P' = 70000,'MIBC15T' = 70000,
                'MIBC16P' = 100000,'MIBC16T' = 100000)

using soupX to remove background

library(SoupX)

soup.mtx <- list()
for(sample in samples){
    dat <- mtx.list[[sample]]
    dat <- dat[, colSums(dat)<=cut.off.max[sample]]
    datCells <- dat[, colSums(dat)>cut.off.min[sample]]
    sc <- SoupChannel(tod = dat, toc = datCells, keepDroplets = TRUE,
                      channelName = sample,
                      dataType = "10X", isV3 = TRUE)
    sc <- estimateSoup(sc, keepDroplets = TRUE,
                       soupRange = c(30,cut.off.max.bk[sample]))
    sc$soupProfile <- sc$soupProfile[ !is.nan(sc$soupProfile$est), ]
    useToEst <- estimateNonExpressingCells(sc,
                                          nonExpressedGeneList = cell.type.markers)
    sc <- calculateContaminationFraction(sc,
                                         cell.type.markers,
                                         useToEst = useToEst)
    soup.mtx[[sample]] <- adjustCounts(sc)
}

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 22.58%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 11.60%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 20.57%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 8.38%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 34.01%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 23.92%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 14.32%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 8.25%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 49.16%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 13.87%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 22.77%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 33.39%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 59.83%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 67.77%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 47.06%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 39.52%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 44.28%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 23.49%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 23.71%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 75.89%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 26.18%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 30.04%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 41.38%

No clusters found or supplied, using every cell as its own cluster.

Estimated global contamination fraction of 74.85%

matrixPath <- '/home/public/ghx/wuLab/Bladder/1217_all/matrices'
for(sample in samples){
    saveRDS(mtx.list[[sample]],
            file = file.path(matrixPath, paste0(sample, '_raw_matrix.rds')))
    saveRDS(soup.mtx[[sample]],
            file = file.path(matrixPath, paste0(sample, '_soup_matrix.rds')))
}

create seurat object and run normal pipeline, DoubletFinder, PhaseDetection, etc.

seu.list <- list()
raw.normalized.data <- list()
for(sample in samples){
    temp.seu <- CreateSeuratObject(counts = soup.mtx[[sample]],
                          project = sample,
                          min.features = 0)
    temp.seu <- Nml_Scl_Fvg_PCA_Fnb(seu = temp.seu)
    temp.seu <- FindClusters(object = temp.seu,
                resolution = 0.25, verbose = F)
    temp.seu <- RunUMAP(object = temp.seu, dims = 1:25, verbose = F,
               umap.method = 'umap-learn', metric = 'correlation')


    # compare background
    raw.seu <- CreateSeuratObject(counts = mtx.list[[sample]][, Cells(temp.seu)],
                          project = 'tool',
                          min.features = 0)
    raw.seu <- NormalizeData(object = raw.seu, verbose = F)
    raw.mtx <- raw.seu@assays$RNA@data[unique(unlist(cell.type.markers)),]
    row.names(raw.mtx) <- paste0('raw_', row.names(raw.mtx))
    raw.mtx <- t(raw.mtx)
    raw.normalized.data[[sample]] <- raw.mtx

    # put raw.mtx into seu as meta.data

    temp.seu@meta.data <- cbind(temp.seu@meta.data,as.data.frame(raw.mtx))
    seu.list[[sample]] <- temp.seu
}

rm(list('raw.seu', 'mtx.list', 'soup.mtx'))

gc()

A matrix: 2 × 6 of type dbl

	used	(Mb)	gc trigger	(Mb)	max used	(Mb)
Ncells	5311972	283.7	15174464	810.5	38464072	2054.3
Vcells	13180174847	100556.8	20737561180	158215.1	17210272365	131304.0

# read former annotations
all.anno <- readRDS(file = paste0('/home/public/ghx/wuLab/Bladder/1101_1st_clean_data/',
                   'All_annotation.rds'))
all.anno$cell.type <- as.character(all.anno$cell.type.1)
anno.cells <- row.names(all.anno)
anno.cells <- paste0(all.anno$orig.ident, '.', anno.cells)
row.names(all.anno) <- anno.cells
# load doublet cells' annotation
db.metadata <- readRDS('/home/public/ghx/wuLab/Bladder/test/doublet_metadata.rds')
db.cells <- row.names(db.metadata)
db.cells <- paste0(db.metadata$orig.ident, '.', db.cells)
row.names(db.metadata) <- db.cells

for(sample in samples){
    intersect.cells <- intersect(anno.cells, Cells(seu.list[[sample]]))
    if(length(intersect.cells)>0){
        cell.type <- rep(NA, length(Cells(seu.list[[sample]])))
        names(cell.type) <- Cells(seu.list[[sample]])
        cell.type[intersect.cells] <- all.anno[intersect.cells,'cell.type']
        cell.type[intersect(Cells(seu.list[[sample]]), db.cells)] <- 'doublet'
        cell.type <- factor(cell.type, levels = c(levels(all.anno$cell.type.1), 'doublet'))
    }else{
        cell.type <- 'undefined'
    }
    seu.list[[sample]]@meta.data$cell.type <- cell.type
}

for(sample in samples){
    saveRDS(seu.list[[sample]],
            file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                             paste0(sample,'_seu.rds')))
}

# compare some gene distribution before and after SoupX removes background
resultPath <- '/home/public/ghx/wuLab/Bladder/1217_all/graph/compareSoup'
for(sample in samples){
    temp.seu <- seu.list[[sample]]
    for(ct in names(cell.type.markers)){
        ct.markers <- cell.type.markers[[ct]]
        ct.markers <- ct.markers[1:min(10, length(ct.markers))]
        raw.ct.markers <- paste0('raw_', ct.markers)
        features.plot <- c()
        for(i in 1:length(ct.markers)){
            features.plot <- c(features.plot,
                               c(ct.markers[i], raw.ct.markers[i]))
        }
        p <- FeaturePlot(object = temp.seu, ncol = 2, pt.size = 0.01,
                         features = features.plot) +
             labs(title = paste0(sample, '-', ct))
        ggsave(p, limitsize = FALSE,
               width = 10, height = 5 * ((length(features.plot)+1) %/% 2),
               file = file.path(resultPath,
                                paste0(sample,'_',ct,'_compare.jpg')))
    }
}

Warning message in FeaturePlot(object = temp.seu, ncol = 2, pt.size = 0.01, features = features.plot):
“All cells have the same value (0) of DES.”
Warning message in FeaturePlot(object = temp.seu, ncol = 2, pt.size = 0.01, features = features.plot):
“All cells have the same value (0) of raw_DES.”
Warning message in FeaturePlot(object = temp.seu, ncol = 2, pt.size = 0.01, features = features.plot):
“All cells have the same value (0) of KRT13.”
Warning message in FeaturePlot(object = temp.seu, ncol = 2, pt.size = 0.01, features = features.plot):
“All cells have the same value (0) of raw_KRT13.”
Warning message in FeaturePlot(object = temp.seu, ncol = 2, pt.size = 0.01, features = features.plot):
“All cells have the same value (0) of AGR2.”
Warning message in FeaturePlot(object = temp.seu, ncol = 2, pt.size = 0.01, features = features.plot):
“All cells have the same value (0) of raw_AGR2.”

# draw FeaturePlot for non-immune and immune cell types
resultPath <- '/home/public/ghx/wuLab/Bladder/1217_all/graph/FeaturePlots'
setwd(resultPath)

#a toy graph shows up if no marker is fed
data(diamonds)
set.seed(42)
small <- diamonds[sample(nrow(diamonds), 1000), ]
toy.p <- ggplot(small)+geom_bar(aes(x=clarity, fill=cut))+coord_polar()

for(sample in samples){
    temp.seu <- seu.list[[sample]]
    genes <- row.names(temp.seu)
    dir.create(sample)
    # overall view of temp sample
    p1 <- DimPlot(object = temp.seu, reduction = 'umap',
        label = T, pt.size = 0.1) +
        labs(title='cluster')
    p2 <- DimPlot(object = temp.seu, reduction = 'umap', group.by = 'cell.type',
        label = T, pt.size = 0.1) +
        labs(title='cell.type/doublets')
    layout <- matrix(c(1,2), ncol = 2, byrow = TRUE)
    p <- arrangeGrob(grobs = list(p1,p2), ncol = 1, layout_matrix = layout,
            top = sample)
    ggsave(plot = p,
           filename = file.path(sample,paste0(sample,'_0.clusters.jpg')),
           width = 10, height = 5)
    # FeaturePlot: non-immune cell markers
    for(ct in names(ni.marker.list)){
        features.plot <- ni.marker.list[[ct]][ni.marker.list[[ct]] %in% genes]
        if(length(features.plot) > 0){
            if(length(features.plot)==1){
                temp.height <- 10
            }else{
                temp.height <- 5 * ((length(features.plot)+1) %/% 2)
            }
            p <- FeaturePlot(object = temp.seu, features = features.plot, pt.size = 0.1, ncol = 2)
            ggsave(plot = p, limitsize = FALSE,
                   filename = file.path(sample,paste0(sample,'_',ct,'.jpg')),
                   width = 10, height = temp.height)
        }
    }
    # FeaturePlot: immune cell markers
    for(ct in names(i.marker.list)){
        features.plot.p <- i.marker.list[[ct]][['p']][i.marker.list[[ct]][['p']] %in% genes]
        features.plot.n <- i.marker.list[[ct]][['n']][i.marker.list[[ct]][['n']] %in% genes]
        if(length(features.plot.p) > 0){
            if(length(features.plot.p)==1){
                temp.p.height <- 2
            }else{
                temp.p.height <- ((length(features.plot.p)+1) %/% 2)
            }
            p1 <- FeaturePlot(object = temp.seu, features = features.plot.p, cols = c('grey','blue'),
                              pt.size = 0.1, ncol = 2)
        }else{
            p1 <- toy.p
            temp.p.height <- 2
        }
        if(length(features.plot.n) > 0){
            if(length(features.plot.n)==1){
                temp.n.height <- 2
            }else{
                temp.n.height <- ((length(features.plot.n)+1) %/% 2)
            }
            p2 <- FeaturePlot(object = temp.seu, features = features.plot.n, cols = c('yellow','red'),
                              pt.size = 0.1, ncol = 2)
        }else{
            p2 <- toy.p
            temp.n.height <- 2
        }
        layout <- matrix(c(rep(c(1,1),temp.p.height),
                           rep(c(2,2), temp.n.height)), ncol = 2, byrow = TRUE)
        p <- arrangeGrob(grobs = list(p1,p2), ncol = 1, layout_matrix = layout,
            top = paste0(sample,'_',ct))
        ggsave(plot = p, limitsize = FALSE,
               filename = file.path(sample,paste0(sample,'_',ct,'.jpg')),
               width = 10, height = 5*(temp.p.height+temp.n.height))
    }
}

Warning message in FeaturePlot(object = temp.seu, features = features.plot, pt.size = 0.1, :
“All cells have the same value (0) of DES.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of CEACAM8.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of EPO.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of EPO.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot, pt.size = 0.1, :
“All cells have the same value (0) of ALAS2.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of CEACAM8.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of EPO.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of LILRA4.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of EPO.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of CEACAM8.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot, pt.size = 0.1, :
“All cells have the same value (0) of KRT13.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot, pt.size = 0.1, :
“All cells have the same value (0) of UPK1B.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of CEACAM8.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of EPO.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of CEACAM8.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of EPO.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of CEACAM8.”
Warning message in FeaturePlot(object = temp.seu, features = features.plot.p, cols = c("grey", :
“All cells have the same value (0) of EPO.”

seu.list

$MIBC5P
An object of class Seurat
33538 features across 8468 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC5T
An object of class Seurat
33538 features across 10714 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC6P
An object of class Seurat
33538 features across 18020 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC6T
An object of class Seurat
33538 features across 20137 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC7P
An object of class Seurat
33538 features across 14824 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC7T
An object of class Seurat
33538 features across 15522 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC8P
An object of class Seurat
33538 features across 22121 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC8T
An object of class Seurat
33538 features across 5612 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC9P
An object of class Seurat
33538 features across 9632 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC9T
An object of class Seurat
33538 features across 9044 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC10P
An object of class Seurat
33538 features across 9662 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC10T
An object of class Seurat
33538 features across 8582 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC11P
An object of class Seurat
33538 features across 18981 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC11T
An object of class Seurat
33538 features across 20796 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC12P
An object of class Seurat
33538 features across 17995 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC12T
An object of class Seurat
33538 features across 15912 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC13P
An object of class Seurat
33538 features across 7082 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC13T
An object of class Seurat
33538 features across 9196 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC14P
An object of class Seurat
33538 features across 6416 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC14T
An object of class Seurat
33538 features across 10024 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC15P
An object of class Seurat
33538 features across 4889 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC15T
An object of class Seurat
33538 features across 10878 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC16P
An object of class Seurat
33538 features across 8895 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

$MIBC16T
An object of class Seurat
33538 features across 12154 samples within 1 assay
Active assay: RNA (33538 features)
 2 dimensional reductions calculated: pca, umap

run DoubletFinder

pk.list <- c('MIBC4P' = 0.01,'MIBC4T' = 0.005,
            'MIBC5P' = 0.005,'MIBC5T' = 0.005,
            'MIBC6P' = 0.005,'MIBC6T' = 0.01,
            'MIBC7P' = 0.005,'MIBC7T' = 0.005,
            'MIBC8P' = 0.005,'MIBC8T' = 0.005,
            'MIBC9P' = 0.005,'MIBC9T' = 0.005,
            'MIBC10P' = 0.005,'MIBC10T' = 0.005,
            'MIBC11P' = 0.005,'MIBC11T' = 0.001,
            'MIBC12P' = 0.001,'MIBC12T' = 0.001,
            'MIBC13P' = 0.001,'MIBC13T' = 0.001,
            'MIBC14P' = 0.001,'MIBC14T' = 0.001,
            'MIBC15P' = 0.001,'MIBC15T' = 0.001,
            'MIBC16P' = 0.001,'MIBC16T' = 0.001
            )

for(sample in samples[15:24]){
    seu.list[[sample]] <- doubletFinder_v3(seu.list[[sample]], PCs = 1:10, pN = 0.25,
                                           pK = pk.list[sample],
                        nExp = 5, reuse.pANN = F, sct = FALSE)
    pann.name <- names(seu.list[[sample]]@meta.data)[ grep(x = names(seu.list[[sample]]@meta.data),
                                        pattern = 'pann', ignore.case = T)]
    df.name <- names(seu.list[[sample]]@meta.data)[ grep(x = names(seu.list[[sample]]@meta.data),
                                            pattern = 'df.class', ignore.case = T)]
    seu.list[[sample]]@meta.data$db.score <- seu.list[[sample]][[pann.name]][,1]
    seu.list[[sample]][[pann.name]] <- NULL
    seu.list[[sample]][[df.name]] <- NULL
}

[1] "Creating 5998 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."
[1] "Creating 5304 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."
[1] "Creating 2361 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."
[1] "Creating 3065 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."
[1] "Creating 2139 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."
[1] "Creating 3341 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."
[1] "Creating 1630 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."
[1] "Creating 3626 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."
[1] "Creating 2965 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."
[1] "Creating 4051 artificial doublets..."
[1] "Creating Seurat object..."
[1] "Normalizing Seurat object..."
[1] "Finding variable genes..."
[1] "Scaling data..."

Centering and scaling data matrix

[1] "Running PCA..."
[1] "Calculating PC distance matrix..."
[1] "Computing pANN..."
[1] "Classifying doublets.."

# plots <- list()
for(sample in samples){
    p <- FeaturePlot(seu.list[[sample]], features = 'db.score',
                     cols = c('yellow','red'), pt.size = 0.1) +
        labs(title = sample)
    plots[[sample]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                 top = 'MIBC5-16:pANN')
ggsave(p, filename = '/home/public/ghx/wuLab/Bladder/1217_all/pANN.jpg',
       limitsize = FALSE,
       width = 10, height = overall.height)

calculate MT-Genes and find dead cells

genes <- row.names(seu.list[['MIBC10P']])
mt.genes <- genes[grep(x = genes, pattern = '^MT-', ignore.case = T)]
for(sample in samples){
    temp.mt.matrix <- seu.list[[sample]]@assays$RNA@counts[mt.genes, ]
    temp.mt.percentage <- colSums(temp.mt.matrix) / colSums(seu.list[[sample]]@assays$RNA@counts)
    seu.list[[sample]]@meta.data$mt.percentage <- temp.mt.percentage
    seu.list[[sample]]@meta.data$dead.cells <- temp.mt.percentage>0.25
}

plots <- list()
for(sample in samples){
    p1 <- FeaturePlot(seu.list[[sample]], features = 'mt.percentage',
                      cols = c('grey','red'), pt.size = 0.01) +
            labs(title = paste0(sample,'-mt.percentage'))
    p2 <- DimPlot(seu.list[[sample]], group.by = 'dead.cells',
                  cols = c('grey','red'), pt.size = 0.01) +
            labs(title = paste0(sample,'-dead.cells.',as.character(sum(seu.list[[sample]]@meta.data$dead.cells))))
    plots[[paste0(sample,'-mt.percentage')]] <- p1
    plots[[paste0(sample,'-dead.cells')]] <- p2
}

p1 <- arrangeGrob(grobs = plots[1:(length(plots)/2)], ncol = 2,
                 top = 'MIBC5-10:dead.cells')
ggsave(p1, filename = '/home/public/ghx/wuLab/Bladder/1217_all/graph/overview/deadCells5-10.jpg',
       limitsize = FALSE,
       width = 10, height = overall.height)
p2 <- arrangeGrob(grobs = plots[(length(plots)/2+1):length(plots)], ncol = 2,
                 top = 'MIBC11-16:dead.cells')
ggsave(p2, filename = '/home/public/ghx/wuLab/Bladder/1217_all/graph/overview/deadCells11-16.jpg',
       limitsize = FALSE,
       width = 10, height = overall.height)

find marker genes for each sample

library(future)

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

structure(function (expr, envir = parent.frame(), substitute = TRUE,
    lazy = FALSE, seed = NULL, globals = TRUE, workers = 8, 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 = 8))

marker.list <- list()
for(sample in samples){
    temp.markers <- FindAllMarkers(object = seu.list[[sample]], verbose = F)
    marker.list[[sample]] <- temp.markers
}

"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC5P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC5T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC6P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC6T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC7P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC7T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC8P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC8T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC9P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC9T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC10P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC10T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC11P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC11T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC12P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC12T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC13P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC13T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC14P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC14T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC15P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC15T_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC16P_cluster_markers.csv"
"","x"
"1","/home/public/ghx/wuLab/Bladder/1217_all/markers/MIBC16T_cluster_markers.csv"

plan("sequential")

for(sample in samples){
    temp.markers <- marker.list[[sample]]
    write.csv(temp.markers, file.path('/home/public/ghx/wuLab/Bladder/1217_all/markers',
                        paste0(sample,'_cluster_markers.csv')))
    tops <- extractTopNMarkers(df = temp.markers, top_n = 10)
    p <- DoHeatmap(object = seu.list[[sample]], features = tops$gene, draw.lines = FALSE,
                   size = 5, angle = 0, disp.min = -2, disp.max = 2) +
            theme(text = element_text(size = 6)) +
            labs(title = sample)
    ggsave(p, width = 10, height = 10,limitsize = FALSE,
           filename = file.path('/home/public/ghx/wuLab/Bladder/1217_all/graph/heatmaps',
                                 paste0(sample, '_heatmap.jpg')))
}

calculate cell cycle scores

plots <- list()
for(sample in samples){
    seu.list[[sample]] <- CellCycleScoring(object = seu.list[[sample]], set.ident = FALSE,
                 s.features = cc.genes$s.genes,
                 g2m.features = cc.genes$g2m.genes)
    p <- DimPlot(seu.list[[sample]], pt.size = 0.01, group.by = 'Phase') +
            labs(title = sample)
    plots[[sample]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                 top = 'MIBC5-16:phase')
ggsave(p, filename = '/home/public/ghx/wuLab/Bladder/1217_all/cellPhase.jpg',
       limitsize = FALSE,
       width = 10, height = overall.height)

Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: MLF1IP, not searching for symbol synonyms”
Warning message:
“The following features are not present in the object: FAM64A, HN1, not searching for symbol synonyms”

annotate clusters of each sample based on FeaturePlot

ct.order <- c('T', 'NK', 'B', 'Pla', 'Mye',
              'Mast','Fib','SMC','Epi','Umb','Endo','UK')

#MIBC5P checked
sample <- 'MIBC5P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,1)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(8,11)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(4)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(2,3)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(6,12,13)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC5T checking - need to increase resolution to separate cluster 7
sample <- 'MIBC5T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,8)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c(2,3,9)] <- 'Pla'

cell.type[seu.list[[sample]]@active.ident %in% c(13)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(11,14)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(1,5)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(12)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(4)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC6P checking - need to increase resolution to separate cluster 2
sample <- 'MIBC6P'
cell.type <- as.character(seu.list[[sample]]@meta.data$cell.type)
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(14)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(5,10)] <- 'SMC'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC6T checked
sample <- 'MIBC6T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,1,3,7)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(2)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c(4,9,13)] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(15)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(14)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC7P checked
sample <- 'MIBC7P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,1)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(2,4)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(12)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC7T checked
sample <- 'MIBC7T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,1)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(2,4)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(12)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC7T checked
sample <- 'MIBC7T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,2,5,7)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(13)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(6,8)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(14)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(1,4)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(15)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC8P checked
sample <- 'MIBC8P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,2,9,12)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(1)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(4)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC8T checked
sample <- 'MIBC8T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,3,4,7)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(13)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c(2)] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(12)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC9P checked
sample <- 'MIBC9P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(1,2)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(4)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(0,5)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC9T checked
#There was a mistake, that I mislabeled cluster13 as plasma cells,
#however it should be Myeloid cells
sample <- 'MIBC9T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,1,2,4)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(15)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c(11,13)] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(16)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(6,12,14)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(8,9)] <- 'UK'


seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC10P checked
sample <- 'MIBC10P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,6)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(2)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(1)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(3,5)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(4)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(12,13,14)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC10T checked
sample <- 'MIBC10T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,2)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(3,4,8)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(1)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC11P checked
sample <- 'MIBC11P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(2)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(1,9)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(4)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(0,7)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(12)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC11T checking - need to increase resolution to separate cluter 5
sample <- 'MIBC11T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(3,7,10,12)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(13)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(17)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(6,8,16,18)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(15)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(14)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(0,1,2,4,5,19)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC12P checked: 7 are dead cells
sample <- 'MIBC12P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(1)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(2)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(0)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(4,10)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC12T checked
sample <- 'MIBC12T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,2,11)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(1,6,8)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(4)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(5,12,13)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC13P checked: 1 are dead cells, 12 are erythrocytes,
sample <- 'MIBC13P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,2,3,4)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(5,7)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(2)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(10,13)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC13T checked: 6 are dead cells, 12&13 are erythrocytes(but in G2S phase)
sample <- 'MIBC13T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,1)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(4)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(3,11)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC14P checked: 8 have G2S phase
sample <- 'MIBC14P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,6,9)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(1)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(4,7)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(2)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(3,11)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(13,14)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(10)] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(12)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC14T checked: 0,1 are background, 4 are dead cells
sample <- 'MIBC14T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(2,5)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(8)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(7)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC15P checed: 4 are dead cells
sample <- 'MIBC15P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c(6)] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c(5)] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c(1)] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c(3)] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(2)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(9,10)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(7,8)] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC15T checking: 3&2(part) are dead cells, needs higher resolution to separate 2
sample <- 'MIBC15T'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c(0,8)] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(9)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

#MIBC16P checing
sample <- 'MIBC16P'
cell.type <- rep(NA,length(Cells(sample)))
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'T'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'NK'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'B'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Pla'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Mye'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Mast'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Fib'
cell.type[seu.list[[sample]]@active.ident %in% c(2)] <- 'Epi'
cell.type[seu.list[[sample]]@active.ident %in% c(9,10)] <- 'Endo'
cell.type[seu.list[[sample]]@active.ident %in% c(7,8)] <- 'SMC'
cell.type[seu.list[[sample]]@active.ident %in% c()] <- 'Umb'
cell.type[seu.list[[sample]]@active.ident %in% c(11)] <- 'UK'

seu.list[[sample]]@meta.data$cell.type <- factor(cell.type,
                levels = ct.order)

plots <- list()
for(sample in samples){
    p <- DimPlot(seu.list[[sample]], label = T, pt.size = 0.1, group.by = 'cell.type') +
         labs(title = sample) #+ theme(legend.position = 'none')
    plots[[sample]] <- p
}
p <- arrangeGrob(grobs = plots, ncol = 2,
                 top = 'MIBC5-16:celltype')
ggsave(p, filename = '/home/public/ghx/wuLab/Bladder/1217_all/graph/overview/cellType.jpg',
       limitsize = FALSE,
       width = 10, height = overall.height)

for(sample in samples){
    temp.df <- seu.list[[sample]]@meta.data
    raw.genes <- row.names(temp.df)[grep(x = row.names(temp.df),
                                         pattern = '^raw_')]
    for(g in raw.genes){
        temp.df[[g]] <- NULL
    }
    saveRDS(temp.df,
            file = file.path('/home/public/ghx/wuLab/Bladder/1217_all/seu',
                            paste0(sample,'_annotation.rds')))
}