该文详细介绍了通过单细胞测序数据对肿瘤微环境中的细胞亚群进行分析的过程，包括细胞类型划分、基因表达特征识别、细胞亚群注释、功能基因集表达分析和细胞亚群浸润分析。主要利用R包Seurat和其他工具，揭示了不同细胞亚群的特性和功能差异。
由于单纯对细胞类型的刻画不足以深入了解肿瘤微环境，采用细胞亚群之间的Marker基因进行重注释，发现独特基因表达的细胞亚群，并进行表型特征刻画。
第一步：导入R包
library(scibet) library(Seurat) library(scater) library(scran) library(dplyr) library(Matrix) library(cowplot) library(ggplot2) library(harmony)
第二步：读入数据（步骤略，这里就是读入（一）中处理后的数据）
第三步：细致划分细胞亚群
将降维聚类写成subType()方法
subType <- function(x){ set.seed(123) x <- NormalizeData(x, normalization.method = "LogNormalize") x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 3000) x <- ScaleData(x, vars.to.regress = c('nCount_RNA')) x <- RunPCA(x) x <- RunHarmony(object = x,group.by.vars = c('Dataset')) x <- RunUMAP(x,reduction = "harmony",dims = 1:30,seed.use = 12345) x <- FindNeighbors(x,reduction = 'harmony', dims = 1:30, verbose = FALSE) x <- FindClusters(x,resolution = 0.5, verbose = FALSE,random.seed=20220727) return(x) } BRCA_SingleCell.Fcell <- subType(BRCA_SingleCell.Fcell)
第四步：筛选细胞亚群差异基因，并清洗细胞亚群（如MKI67为细胞周期的细胞将其剔除）
将差异分析写成findMarker()方法
findMarker <- function(x){ Idents(x) <- 'seurat_clusters' cellType_markerGene <- FindAllMarkers(x,logfc.threshold = 0.5,only.pos = T,test.use = 'roc',min.pct = 0.2) cellType_markerGene <- subset(cellType_markerGene,!grepl(pattern = 'RP[LS]',gene)) cellType_markerGene <- subset(cellType_markerGene,!grepl(pattern = 'MT-',gene)) top5 <- cellType_markerGene %>% group_by(cluster) %>% top_n(n = 10, wt = myAUC) print(top5) return(cellType_markerGene) } cellType_markerGene.Fcell <- findMarker(BRCA_SingleCell.Fcell)展示特征基因
genes_to_check = c('PTPRC', 'CD3D', 'CD3E','BTLA', 'CD4','CD8A','CXCR6','CD27','CD69','ITGAE','CXCL13','LAG3','GZMA', 'GZMK', 'CD79A', 'MS4A1' ,'PDCD1','TIGIT','CTLA4','IL2RA','FOXP3','IL7R','CCR7', 'IGHG1', 'MZB1', 'SDC1', 'CD68', 'CD163', 'CD14', 'JCHAIN','TCF4','IRF4', 'TPSAB1' , 'TPSB2', # mast cells, 'RCVRN','FPR1' , 'VIM' , 'C1QA', 'C1QB', 'S100A9', 'S100A8','SPP1', 'MMP1','ATF4','IL6','THBS2','CYR61','CXCL12', 'LAMP3', 'IDO1','IDO2', 'CD1E','CD1C', 'KRT86','GNLY', 'FGF7','MME', 'ACTA2','MYH11','TAGLN', 'DCN', 'LUM', 'GSN' , 'FAP','FN1','THY1','COL1A1','COL3A1', 'PECAM1', 'VWF', 'EPCAM' , 'CD24', 'KRT19', 'KRT18','MKI67' ) options(repr.plot.height = 11, repr.plot.width = 8) DotPlot(BRCA_SingleCell.Fcell,group.by = 'seurat_clusters', features = unique(genes_to_check),cluster.idents = T) + coord_flip()结果如下：
通过细胞类型的经典Marker清洗混合的细胞
B_P_cell=c("CD79A","MS4A1","IGHG1","JCHAIN") T_cell=c("CD3D","CD3E") Epithelial=c("EPCAM","CD24") Endothelial =c("PECAM1","VWF") Myeloid=c("CD68","CD14") Mast_cell=c("TPSAB1","TPSB2") Fibro_cell=c("LUM","DCN") CellCycle=c("MKI67") top <- cellType_markerGene.Fcell %>% group_by(cluster) %>% top_n(n =10, wt = myAUC) clusters <- top$cluster %>% levels() for(i in clusters){ gene <- top[top$cluster == i,8] %>% unlist(.) if(length(intersect(gene, B_P_cell)) > 0){ print(paste(i, "B cell", sep = " likes ") %>% paste(., intersect(gene, B_P_cell), sep = ":")) } else if(length(intersect(gene, Myeloid)) > 0){ print(paste(i, "Myeloid cell", sep = " likes ") %>% paste(., intersect(gene, Myeloid), sep = ":")) } else if(length(intersect(gene, CellCycle)) > 0){ print(paste(i, "Cell in cycle", sep = " likes ") %>% paste(., intersect(gene, CellCycle), sep = ":")) } else if(length(intersect(gene, Epithelial)) > 0){ print(paste(i, "Epithelial cell", sep = " likes ") %>% paste(., intersect(gene, Epithelial), sep = ":")) } else if(length(intersect(gene, T_cell)) > 0){ print(paste(i, "T cell", sep = " likes ") %>% paste(., intersect(gene, T_cell), sep = ":")) } else if(length(intersect(gene, Mast_cell)) > 0){ print(paste(i, "Mast cell", sep = " likes ") %>% paste(., intersect(gene, Mast_cell), sep = ":")) } else if(length(intersect(gene, Endothelial)) > 0){ print(paste(i, "Endothelial cell", sep = " likes ") %>% paste(., intersect(gene, Endothelial), sep = ":")) } } BRCA_SingleCell.Fcell <- subset(BRCA_SingleCell.Fcell, seurat_clusters %in% c(0,1,2,3,4,7,9,11,12,14,15)) BRCA_SingleCell.Fcell <- subType(BRCA_SingleCell.Fcell)
第五步：根据细胞簇的Marker基因暂时注释细胞亚群
cellType_markerGene.Fcell <- findMarker(BRCA_SingleCell.Fcell) top <- cellType_markerGene.Fcell %>% group_by(cluster) %>% top_n(n =10, wt = myAUC) ####查看Marker基因#### top$cluster %>% levels(.) #查看细胞簇个数 subset(top, cluster == 13, select = gene) #查看Marker基因 ####注释细胞簇#### cluster.1=c(0) #TPM4+F cluster.2=c(1) #COL1A1+F cluster.3=c(2) #MGP+F cluster.4=c(3) #BST2+F cluster.5=c(4) #IGFBP7+F cluster.6=c(5) #LGALS1+F cluster.7=c(6) #ADIRF+F cluster.8=c(7) #FOS+F cluster.9=c(8) #APOD+F cluster.10=c(9) #B2M+F cluster.11=c(10) #GAPDH+F cluster.12=c(11) #HSPA1A+F cluster.13=c(12) #HLA-B+F current.cluster.ids <- c(cluster.1,cluster.2,cluster.3,cluster.4,cluster.5,cluster.6,cluster.7,cluster.8,cluster.9,cluster.10, cluster.11,cluster.12,cluster.13) new.cluster.ids <- c(rep("TPM4+F",length(cluster.1)), rep("COL1A1+F",length(cluster.2)), rep("MGP+F",length(cluster.3)), rep("BST2+F",length(cluster.4)), rep("IGFBP7+F",length(cluster.5)), rep("LGALS1+F",length(cluster.6)), rep("ADIRF+F",length(cluster.7)), rep("FOS+F",length(cluster.8)), rep("APOD+F",length(cluster.9)), rep("B2M+F",length(cluster.10)), rep("GAPDH+F",length(cluster.11)), rep("HSPA1A+F",length(cluster.12)), rep("HLA-B+F",length(cluster.13)) ) BRCA_SingleCell.Fcell@meta.data$cell_Subtype <- plyr::mapvalues(x = BRCA_SingleCell.Fcell$seurat_clusters, from = current.cluster.ids, to = new.cluster.ids)可通过DimPlot()可视化注释后的细胞簇
options(repr.plot.height = 6.5, repr.plot.width = 8) DimPlot(object = BRCA_SingleCell.Fcell,reduction = 'umap',group.by = c('cell_Subtype'), label = TRUE)可视化结果如下：
第六步：根据细胞亚群的经典Marker定义细胞簇
feature.genes <- c("FAP","PDPN","MMP2","PDGFRA","PDGFRL") options(repr.plot.height = 8, repr.plot.width = 10) FeaturePlot(BRCA_SingleCell.Fcell,reduction = "umap",pt.size = .5, features = feature.genes, max.cutoff = 3,min.cutoff = 0) options(repr.plot.height = 8, repr.plot.width = 10) VlnPlot(BRCA_SingleCell.Fcell, features = feature.genes)可视化基因表达：
最后根据表达确认细胞亚群
cluster.1=c(6) #MYH11+MCAM+ myCAFs cluster.2=c(4) #RGS5+MCAM+ myCAFs cluster.3=c(10) #VEGFA+ CAFs cluster.4=c(8) #CFD+ CAFs cluster.5=c(2,7) #PLA2G2A+CFD+ CAFs cluster.6=c(0,1) #PDGFC+ CAFs cluster.7=c(3) #CXCL11+ CAFs cluster.8=c(5,9,11,12) #MMP11+ CAFs current.cluster.ids <- c(cluster.1,cluster.2,cluster.3,cluster.4,cluster.5,cluster.6,cluster.7,cluster.8) new.cluster.ids <- c(rep("MYH11+MCAM+ myCAFs",length(cluster.1)), rep("RGS5+MCAM+ myCAFs",length(cluster.2)), rep("VEGFA+ CAFs",length(cluster.3)), rep("CFD+ CAFs",length(cluster.4)), rep("PLA2G2A+CFD+ CAFs",length(cluster.5)), rep("PDGFC+ CAFs",length(cluster.6)), rep("CXCL11+ CAFs",length(cluster.7)), rep("MMP11+ CAFs",length(cluster.8)) ) BRCA_SingleCell.Fcell@meta.data$cell_Subtype <- plyr::mapvalues(x = BRCA_SingleCell.Fcell$seurat_clusters, from = current.cluster.ids, to = new.cluster.ids) BRCA_SingleCell.Fcell@meta.data$cell_Subtype <- factor(BRCA_SingleCell.Fcell@meta.data$cell_Subtype, levels = c("MYH11+MCAM+ myCAFs","RGS5+MCAM+ myCAFs","VEGFA+ CAFs","CFD+ CAFs", "PLA2G2A+CFD+ CAFs","PDGFC+ CAFs","CXCL11+ CAFs","MMP11+ CAFs")) BRCA_SingleCell.Fcell$cell_Subtype <- factor(BRCA_SingleCell.Fcell$cell_Subtype, levels = c("MYH11+MCAM+ myCAFs","RGS5+MCAM+ myCAFs","VEGF+ CAFs","CFD+ CAFs", "PLA2G2A+CFD+ CAFs","PDGFC+ CAFs","CXCL11+ CAFs","MMP11+ CAFs"), labels = c("MYH11+MCAM+ Myofibro","RGS5+MCAM+ Myofibro","VEGFA+ Fibro","CFD+ Fibro", "PLA2G2A+CFD+ Fibro","PDGFC+ Fibro","CXCL11+ Fibro","MMP11+ Fibro")) options(repr.plot.height = 8.5, repr.plot.width = 10.5) DimPlot(object = BRCA_SingleCell.Fcell,reduction = 'umap',group.by = c('cell_Subtype'), label = TRUE)
第七步：根据功能基因集表达刻画细胞亚群功能
library(aplot) Collagens <- c('COL1A1','COL1A2','COL3A1','COL4A1','COL5A1','COL5A2','COL6A1','COL7A1','COL10A1','COL11A1','COL12A1') ECM <- c('BGN','DCN','LUM','TAGLN','ELN','FN1','FAP','POSTN') MMPs <- c('MMP1','MMP2','MMP3','MMP9','MMP10','MMP11','MMP14','MMP19') TGFb<-c('SERPINE1','CTHRC1','THBS2','THBS1','SULF1') angiogenesis <- c('EGFL6','PDGFC','VEGFA','LOX','TGFBI','ANGPTL4','CA9') Contractile <- c('PDGFA','ACTA2','MYL6','MYH9','MYH11','MCAM','PLN') RAS <- c('RASL12','RASGRP2') proinflammatory <- c('CFD','CFI','C3','C7','CCL21','CXCL14','CXCL12','IL33','IL6','IL7','CXCL3','CXCL2','CCL2') options(repr.plot.height = 10, repr.plot.width = 6) test <- DotPlot(object = BRCA_SingleCell.Fcell, features = c(Collagens,ECM,MMPs,TGFb,angiogenesis,Contractile,RAS,proinflammatory),cols = c('blue','red'),group.by = 'cell_Subtype')+ theme_bw()+ theme(axis.text.x = element_text(angle = -90,hjust = 0,vjust = 0)) + coord_flip() test CAFmarker_data <- test$data CAFmarker_data$Zscore <- CAFmarker_data$avg.exp.scaled CAFmarker_data$cluster <- CAFmarker_data$id CAFmarker_data$geneClass <- NA CAFmarker_data$geneClass <- ifelse(CAFmarker_data$features.plot %in% Collagens,'Collagens',CAFmarker_data$geneClass) CAFmarker_data$geneClass <- ifelse(CAFmarker_data$features.plot %in% ECM,'ECM',CAFmarker_data$geneClass) CAFmarker_data$geneClass <- ifelse(CAFmarker_data$features.plot %in% MMPs,'MMPs',CAFmarker_data$geneClass) CAFmarker_data$geneClass <- ifelse(CAFmarker_data$features.plot %in% TGFb,'TGFb',CAFmarker_data$geneClass) CAFmarker_data$geneClass <- ifelse(CAFmarker_data$features.plot %in% angiogenesis,'angiogenesis',CAFmarker_data$geneClass) CAFmarker_data$geneClass <- ifelse(CAFmarker_data$features.plot %in% Contractile,'Contractile',CAFmarker_data$geneClass) CAFmarker_data$geneClass <- ifelse(CAFmarker_data$features.plot %in% RAS,'RAS',CAFmarker_data$geneClass) CAFmarker_data$geneClass <- ifelse(CAFmarker_data$features.plot %in% proinflammatory,'proinflammatory',CAFmarker_data$geneClass) options(repr.plot.height = 11, repr.plot.width = 4) p1 <- ggplot(subset(CAFmarker_data,geneClass=='proinflammatory'), aes(y=features.plot, x=cluster, fill=Zscore))+ geom_raster()+scale_fill_gradient2(low="#003366", high="#990033", mid="white")+ xlab(NULL) + ylab(NULL)+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),legend.position="none",axis.ticks = element_blank(),axis.text.x = element_blank())+ geom_vline(xintercept=c(-0.5,2.5,3.5,5.5),size=.8) p3 <- ggplot(subset(CAFmarker_data,geneClass=='RAS'), aes(y=features.plot, x=cluster, fill=Zscore))+ geom_raster()+scale_fill_gradient2(low="#003366", high="#990033", mid="white")+ xlab(NULL) + ylab(NULL)+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),legend.position="none",axis.ticks = element_blank(),axis.text.x = element_blank())+ geom_vline(xintercept=c(-0.5,2.5,3.5,5.5),size=.8) p4 <- ggplot(subset(CAFmarker_data,geneClass=='Contractile'), aes(y=features.plot, x=cluster, fill=Zscore))+ geom_raster()+scale_fill_gradient2(low="#003366", high="#990033", mid="white")+ xlab(NULL) + ylab(NULL)+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),legend.position="none",axis.ticks = element_blank(),axis.text.x = element_blank())+ geom_vline(xintercept=c(-0.5,2.5,3.5,5.5),size=.8) p5 <- ggplot(subset(CAFmarker_data,geneClass=='angiogenesis'), aes(y=features.plot, x=cluster, fill=Zscore))+ geom_raster()+scale_fill_gradient2(low="#003366", high="#990033", mid="white")+ xlab(NULL) + ylab(NULL)+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),legend.position="none",axis.ticks = element_blank(),axis.text.x = element_blank())+ geom_vline(xintercept=c(-0.5,2.5,3.5,5.5),size=.8) p6 <- ggplot(subset(CAFmarker_data,geneClass=='TGFb'), aes(y=features.plot, x=cluster, fill=Zscore))+ geom_raster()+scale_fill_gradient2(low="#003366", high="#990033", mid="white")+ xlab(NULL) + ylab(NULL)+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),legend.position="none",axis.ticks = element_blank(),axis.text.x = element_blank())+ geom_vline(xintercept=c(-0.5,2.5,3.5,5.5),size=.8) p7 <- ggplot(subset(CAFmarker_data,geneClass=='MMPs'), aes(y=features.plot, x=cluster, fill=Zscore))+ geom_raster()+scale_fill_gradient2(low="#003366", high="#990033", mid="white")+ xlab(NULL) + ylab(NULL)+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),legend.position="none",axis.ticks = element_blank(),axis.text.x = element_blank())+ geom_vline(xintercept=c(-0.5,2.5,3.5,5.5),size=.8) p8 <- ggplot(subset(CAFmarker_data,geneClass=='ECM'), aes(y=features.plot, x=cluster, fill=Zscore))+ geom_raster()+scale_fill_gradient2(low="#003366", high="#990033", mid="white")+ xlab(NULL) + ylab(NULL)+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),legend.position="none",axis.ticks = element_blank(),axis.text.x = element_blank())+ geom_vline(xintercept=c(-0.5,2.5,3.5,5.5),size=.8) p9 <- ggplot(subset(CAFmarker_data,geneClass=='Collagens'), aes(y=features.plot, x=cluster, fill=Zscore))+ geom_raster()+scale_fill_gradient2(low="#003366", high="#990033", mid="white")+ theme(axis.text.x = element_text(angle = -90,hjust = 0,vjust = 0.5),axis.ticks.y = element_blank())+ xlab(NULL) + ylab(NULL)+ theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))+ geom_vline(xintercept=c(-0.5,2.5,3.5,5.5),size=.8) p10<- p1 %>% insert_bottom(p3,height = length(RAS)/length(proinflammatory)) %>% insert_bottom(p4,height = length(Contractile)/length(proinflammatory)) %>% insert_bottom(p7,height = length(MMPs)/length(proinflammatory)) %>% insert_bottom(p5,height = length(angiogenesis)/length(proinflammatory)) %>% insert_bottom(p6,height = length(TGFb)/length(proinflammatory)) %>% insert_bottom(p8,height = length(ECM)/length(proinflammatory)) %>% insert_bottom(p9,height = length(Collagens)/length(proinflammatory)) p10结果如下：
第八步：细胞亚群浸润分析
library(ggpubr) library(ggsignif) BRCA_Unpaired <- c(BRCA_TNBC,BRCA_HER.2,BRCA_PR,BRCA_ER,BRCA_TNBC_BRCA,BRCA_BRCA,BRCA_NM) BRCA_Paired <- c(BRCA_TNBC_Primary,BRCA_TNBC_Lymph,BRCA_Luminal.B_Primary,BRCA_Luminal.B_Lymph,BRCA_HER.2_Primary,BRCA_HER.2_Lymph, BRCA_ER_Primary,BRCA_ER_Lymph) BRCA_SingleCell.Fcell@meta.data$Paired <- "Unknown" BRCA_SingleCell.Fcell@meta.data$Paired[BRCA_SingleCell.Fcell@meta.data$orig.ident %in% BRCA_Unpaired] <- "Unpaired" BRCA_SingleCell.Fcell@meta.data$Paired[BRCA_SingleCell.Fcell@meta.data$orig.ident %in% BRCA_Paired] <- "Paired" cell.Paired <- aggregate(BRCA_SingleCell.Fcell$TumorType, list(BRCA_SingleCell.Fcell$cell_Subtype, BRCA_SingleCell.Fcell$TumorType, BRCA_SingleCell.Fcell$Patient.id, BRCA_SingleCell.Fcell$Paired), length) patient <- unique(cell.Paired$Group.3) for(i in patient){ #Lymph if(i == patient[1]){ cell.pro.Paired.L <- subset(cell.Paired, (Group.3 == i)&(Group.2 == "Lymph")) cell.pro.Paired.L$proportion <- cell.pro.Paired.L$x/sum(cell.pro.Paired.L$x) } else { single.Paired.L <- subset(cell.Paired, (Group.3 == i)&(Group.2 == "Lymph")) single.Paired.L$proportion <- single.Paired.L$x/sum(single.Paired.L$x) cell.pro.Paired.L <- rbind(cell.pro.Paired.L, single.Paired.L) } #Primary if(i == patient[1]){ cell.pro.Paired.P <- subset(cell.Paired, (Group.3 == i)&(Group.2 == "Primary")) cell.pro.Paired.P$proportion <- cell.pro.Paired.P$x/sum(cell.pro.Paired.P$x) } else { single.Paired.P <- subset(cell.Paired, (Group.3 == i)&(Group.2 == "Primary")) single.Paired.P$proportion <- single.Paired.P$x/sum(single.Paired.P$x) cell.pro.Paired.P <- rbind(cell.pro.Paired.P, single.Paired.P) } } cell.pro.Paired <- rbind(cell.pro.Paired.L, cell.pro.Paired.P) head(cell.Paired) head(cell.pro.Paired) cell.Paired <- aggregate(BRCA_SingleCell.Fcell$TumorType, list(BRCA_SingleCell.Fcell$cell_Subtype, BRCA_SingleCell.Fcell$TumorType, BRCA_SingleCell.Fcell$Patient.id, BRCA_SingleCell.Fcell$Paired), length) patient <- unique(cell.Paired$Group.3) for(i in patient){ #Lymph if(i == patient[1]){ cell.pro.Paired.L <- subset(cell.Paired, (Group.3 == i)&(Group.2 == "Lymph")) cell.pro.Paired.L$proportion <- cell.pro.Paired.L$x/sum(cell.pro.Paired.L$x) } else { single.Paired.L <- subset(cell.Paired, (Group.3 == i)&(Group.2 == "Lymph")) single.Paired.L$proportion <- single.Paired.L$x/sum(single.Paired.L$x) cell.pro.Paired.L <- rbind(cell.pro.Paired.L, single.Paired.L) } #Primary if(i == patient[1]){ cell.pro.Paired.P <- subset(cell.Paired, (Group.3 == i)&(Group.2 == "Primary")) cell.pro.Paired.P$proportion <- cell.pro.Paired.P$x/sum(cell.pro.Paired.P$x) } else { single.Paired.P <- subset(cell.Paired, (Group.3 == i)&(Group.2 == "Primary")) single.Paired.P$proportion <- single.Paired.P$x/sum(single.Paired.P$x) cell.pro.Paired.P <- rbind(cell.pro.Paired.P, single.Paired.P) } } cell.pro.Paired <- rbind(cell.pro.Paired.L, cell.pro.Paired.P) head(cell.Paired) head(cell.pro.Paired) options(repr.plot.height = 6, repr.plot.width = 7) ggboxplot(cell.pro.Paired, x = "Group.1", y = "proportion", color = "Group.2", palette = "jco", add = "jitter")+# palette可以按照期刊选择相应的配色，如"npg"等 stat_compare_means(aes(group = Group.2), label = "p.format")+ theme(panel.grid.major = element_line(colour = NA), panel.grid.minor = element_blank(), text=element_text(size = 12, family = "serif"), axis.text.x = element_text(angle = 90))结果如下：
第九步：HallMarker功能差异分析
library(msigdbr) h.human <- msigdbr(species = "Homo sapiens", category = "H") h.names <- unique(h.human$gs_name) h.sets <- vector("list", length=length(h.names)) names(h.sets) <- h.names for(i in names(h.sets)){ h.sets[[i]] <- subset(h.human, gs_name == i, "gene_symbol") %>% unlist(.) } map_names = function(seur=NULL, names=NULL){ # Map "names" to genes or features in a Seurat object # --------------------------------------------------- # seur = seurat object # names = list of names (genes or features) # returns list(genes=c(GENES), feats=(FEATS)) # Initialize variables names = as.list(names) genes = c() feats = c() # Get data and metadata data = seur@assays[['RNA']]@data meta = seur@meta.data # Map names if(!is.null(names)){ genes = sapply(names, function(a){intersect(a, rownames(data))}, simplify=F) feats = sapply(names, function(a){intersect(a, colnames(meta))}, simplify=F) } # Filter genes and feats genes = genes[lengths(genes) > 0] feats = feats[lengths(feats) > 0] # Fix gene names if(length(genes) > 0){ if(is.null(names(genes))){names(genes) = sapply(genes, paste, collapse='.')} } # Fix feat names if(length(feats) > 0){ if(is.null(names(feats))){names(feats) = sapply(feats, paste, collapse='.')} } return(list(genes=genes, feats=feats)) } score_cells = function(seur=NULL, names=NULL, combine_genes='mean', groups=NULL, group_stat='mean', cells.use=NULL){ # Score genes and features across cells and optionally aggregate # The steps are: # - calculate mean expression across genes (combine_genes = 'sum', 'mean', 'scale', 'scale2') # - calculate mean expression within each cell type (groups, group_stat = 'mean', 'alpha', 'mu') # Fix input arguments and get data for name mapping data = seur@assays[['RNA']]@data meta = seur@meta.data if(!is.null(groups)){groups = setNames(groups, colnames(seur@assays[['RNA']]@data))} scores = NULL # Map genes and feats res = map_names(seur=seur, names=names) genes = res$genes feats = res$feats genes.use = unique(do.call(c, genes)) # Subset cells if(!is.null(cells.use)){ data = data[,cells.use,drop=F] meta = meta[cells.use,,drop=F] if(!is.null(groups)){groups = groups[cells.use]} } group_genes = function(x, method){ # combine expression data across genes within a signature # x = [genes x cells] matrix # method = 'sum', 'mean', 'scale' # returns [genes x cells] or [1 x cells] matrix if(nrow(x) == 1){return(x[1,,drop=F])} if(method == 'sum'){ t(colSums(x)) } else if(method == 'mean'){ t(colMeans(x, na.rm=T)) } else if(method == 'scale'){ x = t(scale(t(x))) t(colMeans(x, na.rm=T)) } else if(method == 'scale2'){ x = t(scale(t(x), center=F)) t(colMeans(x, na.rm=T)) } else if(method == 'none'){ x } else { stop('Error: invalid combine_genes method') } } group_cells = function(x, groups, method){ # combine expression data across cells # x = [genes x cells] matrix # group_stat = 'alpha', 'mu', or 'mean' # returns [genes x groups] matrix if(is.null(groups)){return(x)} if(method %in% c('n', 'sum')){ if(method == 'n'){x = x > 0} x = t(data.frame(aggregate(t(x), list(groups), sum, na.rm=T), row.names=1)) } else { if(method == 'alpha'){x = x > 0} if(method == 'mu'){x[x == 0] = NA} x = t(data.frame(aggregate(t(x), list(groups), mean, na.rm=T), row.names=1)) } x[is.na(x)] = 0 x } # Calculate scores names.use = unique(c(names(genes), names(feats))) # Speed improvements (fast indexing for flat structures) name_map = sapply(names.use, function(a) c(genes[[a]], feats[[a]]), simplify=F) do.flat = all(lengths(name_map) == 1) if(do.flat == TRUE){ genes[['flat']] = do.call(c, genes) feats[['flat']] = do.call(c, feats) names.iter = 'flat' combine_genes = 'none' } else { names.iter = names.use } backup = scores scores = lapply(names.iter, function(name){ # Combine data and metadata if(name %in% names(genes)){ si = data[genes[[name]],,drop=F] } else { si = c() } if(name %in% names(feats)){ if(is.null(si)){ si = t(meta[,feats[[name]],drop=F]) } else { si = rBind(si, t(meta[,feats[[name]],drop=F])) } } si = as.matrix(as.data.frame(si)) si = group_genes(si, method=combine_genes) si = group_cells(si, groups=groups, method=group_stat) si = data.frame(t(si)) }) # Collapse scores if(do.flat == TRUE){ scores = scores[[1]][,make.names(name_map[names.use]),drop=F] } else { do.collapse = all(lapply(scores, ncol) == 1) if(do.collapse == TRUE){ scores = as.data.frame(do.call(cbind, scores)) } } # Fix names names.use = make.names(names.use) names(scores) = names.use # Combine data if(!is.null(backup)){ if(is.data.frame(scores)){ scores = cbind.data.frame(scores, backup) } else { scores = c(scores, backup) } } if(nrow(scores) == 0){scores = NULL} return(scores) } map_names = function(seur=NULL, names=NULL){ # Map "names" to genes or features in a Seurat object # --------------------------------------------------- # seur = seurat object # names = list of names (genes or features) # returns list(genes=c(GENES), feats=(FEATS)) # Initialize variables names = as.list(names) genes = c() feats = c() # Get data and metadata data = seur@assays[['RNA']]@data meta = seur@meta.data # Map names if(!is.null(names)){ genes = sapply(names, function(a){intersect(a, rownames(data))}, simplify=F) feats = sapply(names, function(a){intersect(a, colnames(meta))}, simplify=F) } # Filter genes and feats genes = genes[lengths(genes) > 0] feats = feats[lengths(feats) > 0] # Fix gene names if(length(genes) > 0){ if(is.null(names(genes))){names(genes) = sapply(genes, paste, collapse='.')} } # Fix feat names if(length(feats) > 0){ if(is.null(names(feats))){names(feats) = sapply(feats, paste, collapse='.')} } return(list(genes=genes, feats=feats)) } score_cells = function(seur=NULL, names=NULL, combine_genes='mean', groups=NULL, group_stat='mean', cells.use=NULL){ # Score genes and features across cells and optionally aggregate # The steps are: # - calculate mean expression across genes (combine_genes = 'sum', 'mean', 'scale', 'scale2') # - calculate mean expression within each cell type (groups, group_stat = 'mean', 'alpha', 'mu') # Fix input arguments and get data for name mapping data = seur@assays[['RNA']]@data meta = seur@meta.data if(!is.null(groups)){groups = setNames(groups, colnames(seur@assays[['RNA']]@data))} scores = NULL # Map genes and feats res = map_names(seur=seur, names=names) genes = res$genes feats = res$feats genes.use = unique(do.call(c, genes)) # Subset cells if(!is.null(cells.use)){ data = data[,cells.use,drop=F] meta = meta[cells.use,,drop=F] if(!is.null(groups)){groups = groups[cells.use]} } group_genes = function(x, method){ # combine expression data across genes within a signature # x = [genes x cells] matrix # method = 'sum', 'mean', 'scale' # returns [genes x cells] or [1 x cells] matrix if(nrow(x) == 1){return(x[1,,drop=F])} if(method == 'sum'){ t(colSums(x)) } else if(method == 'mean'){ t(colMeans(x, na.rm=T)) } else if(method == 'scale'){ x = t(scale(t(x))) t(colMeans(x, na.rm=T)) } else if(method == 'scale2'){ x = t(scale(t(x), center=F)) t(colMeans(x, na.rm=T)) } else if(method == 'none'){ x } else { stop('Error: invalid combine_genes method') } } group_cells = function(x, groups, method){ # combine expression data across cells # x = [genes x cells] matrix # group_stat = 'alpha', 'mu', or 'mean' # returns [genes x groups] matrix if(is.null(groups)){return(x)} if(method %in% c('n', 'sum')){ if(method == 'n'){x = x > 0} x = t(data.frame(aggregate(t(x), list(groups), sum, na.rm=T), row.names=1)) } else { if(method == 'alpha'){x = x > 0} if(method == 'mu'){x[x == 0] = NA} x = t(data.frame(aggregate(t(x), list(groups), mean, na.rm=T), row.names=1)) } x[is.na(x)] = 0 x } # Calculate scores names.use = unique(c(names(genes), names(feats))) # Speed improvements (fast indexing for flat structures) name_map = sapply(names.use, function(a) c(genes[[a]], feats[[a]]), simplify=F) do.flat = all(lengths(name_map) == 1) if(do.flat == TRUE){ genes[['flat']] = do.call(c, genes) feats[['flat']] = do.call(c, feats) names.iter = 'flat' combine_genes = 'none' } else { names.iter = names.use } backup = scores scores = lapply(names.iter, function(name){ # Combine data and metadata if(name %in% names(genes)){ si = data[genes[[name]],,drop=F] } else { si = c() } if(name %in% names(feats)){ if(is.null(si)){ si = t(meta[,feats[[name]],drop=F]) } else { si = rBind(si, t(meta[,feats[[name]],drop=F])) } } si = as.matrix(as.data.frame(si)) si = group_genes(si, method=combine_genes) si = group_cells(si, groups=groups, method=group_stat) si = data.frame(t(si)) }) # Collapse scores if(do.flat == TRUE){ scores = scores[[1]][,make.names(name_map[names.use]),drop=F] } else { do.collapse = all(lapply(scores, ncol) == 1) if(do.collapse == TRUE){ scores = as.data.frame(do.call(cbind, scores)) } } # Fix names names.use = make.names(names.use) names(scores) = names.use # Combine data if(!is.null(backup)){ if(is.data.frame(scores)){ scores = cbind.data.frame(scores, backup) } else { scores = c(scores, backup) } } if(nrow(scores) == 0){scores = NULL} return(scores) } row_AUC <- rownames(AUC) AUC <- cbind(AUC, cells = row_AUC) cells_meta <- data.frame(cells = rownames(BRCA_SingleCell.Fcell@meta.data), cell_Subtype = BRCA_SingleCell.Fcell@meta.data$cell_Subtype, TumorType = BRCA_SingleCell.Fcell@meta.data$TumorType) cells_meta %>% head() AUC <- merge(cells_meta, AUC, by = "cells") AUC <- AUC[,-1] row_AUC <- rownames(AUC) AUC <- cbind(AUC, cells = row_AUC) cells_meta <- data.frame(cells = rownames(BRCA_SingleCell.Fcell@meta.data), cell_Subtype = BRCA_SingleCell.Fcell@meta.data$cell_Subtype, TumorType = BRCA_SingleCell.Fcell@meta.data$TumorType) cells_meta %>% head() AUC <- merge(cells_meta, AUC, by = "cells") AUC <- AUC[,-1] pathway.name <- AUC_Cells_diff$pathway pathway.name <- strsplit(pathway.name, split = "HALLMARK_") pathway.name <- unlist(pathway.name) pathway.name <- pathway.name[seq(2,length(pathway.name), by = 2)] AUC_Cells_diff$pathway <- pathway.name AUC_Cells_diff <- AUC_Cells_diff[AUC_Cells_diff$pathway %in% c("ALLOGRAFT_REJECTION","COAGULATION","COMPLEMENT", "IL6_JAK_STAT3_SIGNALING","INFLAMMATORY_RESPONSE", "INTERFERON_ALPHA_RESPONSE","INTERFERON_GAMMA_RESPONSE", #Immune "BILE_ACID_METABOLISM","CHOLESTEROL_HOMEOSTASIS","FATTY_ACID_METABOLISM", "GLYCOLYSIS","HEME_METABOLISM","OXIDATIVE_PHOSPHORYLATION", "XENOBIOTIC_METABOLISM", #Metabolism "E2F_TARGETS","G2M_CHECKPOINT","MITOTIC_SPINDLE","MYC_TARGETS_V1", "MYC_TARGETS_V2", #Proliferation "ANGIOGENESIS","EPITHELIAL_MESENCHYMAL_TRANSITION", #Metastasis "ANDROGEN_RESPONSE","APOPTOSIS","ESTROGEN_RESPONSE_EARLY", "ESTROGEN_RESPONSE_LATE","HEDGEHOG_SIGNALING","HYPOXIA", "IL2_STAT5_SIGNALING","KRAS_SIGNALING_DN","KRAS_SIGNALING_UP", "MTORC1_SIGNALING","NOTCH_SIGNALING","PI3K_AKT_MTOR_SIGNALING", "PROTEIN_SECRETION","REACTIVE_OXYGEN_SPECIES_PATHWAY","TGF_BETA_SIGNALING", "TNFA_SIGNALING_VIA_NFKB","UNFOLDED_PROTEIN_RESPONSE", "WNT_BETA_CATENIN_SIGNALING"),] #Signaling AUC_Cells_diff$pathway <- factor(AUC_Cells_diff$pathway, levels = c("ALLOGRAFT_REJECTION","COAGULATION","COMPLEMENT", "IL6_JAK_STAT3_SIGNALING","INFLAMMATORY_RESPONSE", "INTERFERON_ALPHA_RESPONSE","INTERFERON_GAMMA_RESPONSE", #Immune "BILE_ACID_METABOLISM","CHOLESTEROL_HOMEOSTASIS","FATTY_ACID_METABOLISM", "GLYCOLYSIS","HEME_METABOLISM","OXIDATIVE_PHOSPHORYLATION", "XENOBIOTIC_METABOLISM", #Metabolism "E2F_TARGETS","G2M_CHECKPOINT","MITOTIC_SPINDLE","MYC_TARGETS_V1", "MYC_TARGETS_V2", #Proliferation "ANGIOGENESIS","EPITHELIAL_MESENCHYMAL_TRANSITION", #Metastasis "ANDROGEN_RESPONSE","APOPTOSIS","ESTROGEN_RESPONSE_EARLY", "ESTROGEN_RESPONSE_LATE","HEDGEHOG_SIGNALING","HYPOXIA", "IL2_STAT5_SIGNALING","KRAS_SIGNALING_DN","KRAS_SIGNALING_UP", "MTORC1_SIGNALING","NOTCH_SIGNALING","PI3K_AKT_MTOR_SIGNALING", "PROTEIN_SECRETION","REACTIVE_OXYGEN_SPECIES_PATHWAY","TGF_BETA_SIGNALING", "TNFA_SIGNALING_VIA_NFKB","UNFOLDED_PROTEIN_RESPONSE", "WNT_BETA_CATENIN_SIGNALING"),#Signaling labels = c("Allograft rejection","Coagulation","Complement", "IL6 JAK STAT3 signaling","Inflammatory response", "Interferon alpha response","Interferon gamma response", "Bile acid metabolism","Cholesterol homeostasis","Fatty acid matebolism", "Glycolysis","Heme Metabolism","Oxidative phosphorylation", "Xenobiotic metabolism", "E2f tagets","G2M checkpoint","Mitotic spindle","MYC targets v1", "MYC targets v2", "Angiogenesis","Epithelial mesenchymal transition", "Androgen response","Apoptosis","Estrogen response early", "Estrogen response late","Hedgehog signaling","Hypoxia", "IL2 STAT5 signaling","KRAS signaling dn","KRAS signaling up", "mTORC1 signaling","Notch signaling","PI3K AKT MTOR signaling", "Protein secretion","Reactive oxygen species pathway","TGF beta signaling", "TNFA signaling via NFKB","Unfolded protein response", "WNT beta catenin signaling")) FDR <- AUC_Cells_diff[is.finite(AUC_Cells_diff$FDR),1] LOG2FC <- AUC_Cells_diff[is.finite(AUC_Cells_diff$LOG2FC),2] AUC_Cells_diff[is.infinite(AUC_Cells_diff$LOG2FC),2] <- ifelse(AUC_Cells_diff[is.infinite(AUC_Cells_diff$LOG2FC),2] < 0, -max(abs(LOG2FC)), max(abs(LOG2FC))) options(repr.plot.height = 7, repr.plot.width = 15) ggplot(AUC_Cells_diff, aes(pathway, forcats::fct_rev(pre_cellType), fill = LOG2FC, size = FDR)) + geom_point(shape = 21, stroke = 0.1) + geom_hline(yintercept = seq(1.5, 9.5, 1), size = 0.5, col = "black") + geom_vline(xintercept = seq(1.5, 38.5, 1), size = 0.5, col = "black") + scale_x_discrete(position = "top") + scale_radius(range = c(2, 6)) + scale_fill_gradient2(low = "blue", mid = "white",high = "red", breaks = c(-round(max(abs(AUC_Cells_diff$LOG2FC)),1), 0, round(max(abs(AUC_Cells_diff$LOG2FC)),1)), limits = c(-round(max(abs(AUC_Cells_diff$LOG2FC)),1), round(max(abs(AUC_Cells_diff$LOG2FC)),1))) + #scale_fill_gradient(low = "blue" , high = "red") + #scale_fill_distiller(palette = "Spectral") + #theme_minimal() + theme(legend.position = "bottom", panel.grid.major = element_blank(), axis.ticks.y = element_blank(), axis.ticks.x = element_blank(), legend.text = element_text(size = 8), legend.title = element_text(size = 8), axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5, size = 10), axis.text.y = element_text(size = 12)) + guides(size = guide_legend(override.aes = list(fill = NA, color = "black", stroke = .25), label.position = "bottom", title.position = "right", order = 1), fill = guide_colorbar(ticks.colour = NA, title.position = "top", order = 2)) + labs(size = "FDR", fill = "logFC", x = NULL, y = NULL)结果类似如下：（随便找了个图，反正可以调整）
至此，本节内容结束，主要是通过对细胞进行重注释，分出细胞亚群，从而进行功能相关的分析，以刻画细胞亚群的部分特征。部分代码可能需要修改才能使用哦！主要是观其教程，通其意！
扫一扫
2540
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
