（三）单细胞数据分析——细胞亚群的表型特征刻画
2024-05-14 03:14:37
由于单纯对细胞类型的刻画不足以深入了解肿瘤微环境，采用细胞亚群之间的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+ CAFscurrent.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){#Lymphif(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)}#Primaryif(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){#Lymphif(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)}#Primaryif(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 variablesnames = as.list(names)genes = c()feats = c()# Get data and metadatadata = seur@assays[['RNA']]@datameta = seur@meta.data# Map namesif(!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 featsgenes = genes[lengths(genes) > 0]feats = feats[lengths(feats) > 0]# Fix gene namesif(length(genes) > 0){if(is.null(names(genes))){names(genes) = sapply(genes, paste, collapse='.')}}# Fix feat namesif(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 mappingdata = seur@assays[['RNA']]@datameta = seur@meta.dataif(!is.null(groups)){groups = setNames(groups, colnames(seur@assays[['RNA']]@data))}scores = NULL# Map genes and featsres = map_names(seur=seur, names=names)genes = res$genesfeats = res$featsgenes.use = unique(do.call(c, genes))# Subset cellsif(!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] matrixif(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] matrixif(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)] = 0x}# Calculate scoresnames.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 = scoresscores = lapply(names.iter, function(name){# Combine data and metadataif(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 scoresif(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 namesnames.use = make.names(names.use)names(scores) = names.use# Combine dataif(!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 variablesnames = as.list(names)genes = c()feats = c()# Get data and metadatadata = seur@assays[['RNA']]@datameta = seur@meta.data# Map namesif(!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 featsgenes = genes[lengths(genes) > 0]feats = feats[lengths(feats) > 0]# Fix gene namesif(length(genes) > 0){if(is.null(names(genes))){names(genes) = sapply(genes, paste, collapse='.')}}# Fix feat namesif(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 mappingdata = seur@assays[['RNA']]@datameta = seur@meta.dataif(!is.null(groups)){groups = setNames(groups, colnames(seur@assays[['RNA']]@data))}scores = NULL# Map genes and featsres = map_names(seur=seur, names=names)genes = res$genesfeats = res$featsgenes.use = unique(do.call(c, genes))# Subset cellsif(!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] matrixif(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] matrixif(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)] = 0x}# Calculate scoresnames.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 = scoresscores = lapply(names.iter, function(name){# Combine data and metadataif(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 scoresif(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 namesnames.use = make.names(names.use)names(scores) = names.use# Combine dataif(!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"),#Signalinglabels = 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)结果类似如下：（随便找了个图，反正可以调整）
至此，本节内容结束，主要是通过对细胞进行重注释，分出细胞亚群，从而进行功能相关的分析，以刻画细胞亚群的部分特征。部分代码可能需要修改才能使用哦！主要是观其教程，通其意！
（三）单细胞数据分析——细胞亚群的表型特征刻画相关推荐
- 单细胞时代 || 细胞身份概念的演变
原创 周运来 单细胞天地 2021-02-12 17:30 收录于合集#单细胞时代7个 作者 | 周运来 男, 一个长大了才会遇到的帅哥, 稳健,潇洒,大方,靠谱. 一段生信缘,一棵技能树. 生信技 ...
- 单细胞数据分析工具scvi介绍
之前听师兄在组会上介绍了scvi,恰巧最近在用scvi做一些分析,感觉这个方法很精巧效果也不错,因此研究了一下原始论文,这里简要介绍一下这个单细胞数据分析工具的新星. scvi是Nature Meth ...
- seurat提取表达矩阵_单细胞数据分析神器——Seurat
近年来,单细胞技术日益火热,并且有着愈演愈烈的趋势.在2015年至2017年,甚至对某细胞群体或组织进行单细胞测序,解析其细胞成分就能发一篇CNS级别的文章.近两三年,单细胞技术从最开始的基因组,转录 ...
- 单细胞测序流程（六）单细胞的细胞类型的注释
系列文章目录 单细胞测序流程(一)简介与数据下载 单细胞测序流程(二)数据整理 单细胞测序流程(三)质控和数据过滤--Seurat包分析,小提琴图和基因离差散点图 单细胞测序流程(四)主成分分析--P ...
- 用GraphPad Prism直出单细胞亚群堆叠柱状图
单细胞亚群分析后,虽然可以用R语言(如ggplot2包)出单细胞亚群的堆叠柱状图,今天想介绍一种零代码的Prism出图方式. 1.准备好亚群数据. R中,通过 table(test.seu$cellt ...
- 单细胞测序流程（七）单细胞的细胞类型轨迹分析
系列文章目录 单细胞测序流程(一)简介与数据下载 单细胞测序流程(二)数据整理 单细胞测序流程(三)质控和数据过滤--Seurat包分析,小提琴图和基因离差散点图 单细胞测序流程(四)主成分分析--P ...
- 三种方式细胞评分对比
三种方式细胞评分对比 library(patchwork) library(ggplot2) library(ggalluvial) library(svglite) library(Seurat) ...
- spss分析qpcr数据_实时荧光定量PCR的三种数据分析方法比较.doc
窑128窑热带病与寄生虫学 圆园12 年第 10 卷第 3 期 允燥怎则灶葬造 燥枣 栽则燥责蚤糟葬造 阅蚤泽藻葬泽藻泽 葬灶凿 孕葬则葬泽蚤贼燥造燥早赠 圆园12援 V燥l 10. 晕燥 3 doi ...
- 大型单细胞数据分析解决方案
为什么要做大型单细胞数据分析 因为单细胞数据在呈指数增长,遇到大数据集只是早晚的问题.曾经我们困惑一个物种的基因组那么大,如果给很多物种都测基因组的话,拿什么来存储这些数据?随着单细胞技术的成熟,测序 ...
最新文章
- 【转】初等数论 ——原根、指标及其应用
- VMware Workstation 9下基于Ubuntu 12.10服务器版本的Hadoop集群的配置
- Android之提示Could not find com.android.support:appcompat-v7:25.3.1.
- 《Android游戏开发详解》一2.17 对象是独立的
- Debian for ARM install python 3.5.x
- 家用中央空调设计浅议
- iframe 透明参数
- java 两张图片合并_java实现把两张图片合并(Graphics2D)
- python中的成员运算符用于判断指定_Python中的成员运算符用于判断指定序列中是否包含某个值...
- winrar破解注册
- vue vue-quill-editor 富文本 改变图片大小
- S32K3 MCAL WDG看门狗配置详解 基于EBtresos
- linux u盘文件乱码,轻松解决Linux下U盘乱码的方法
- 12.20 沙牛家书 《不负牛市不负沙》
- 表情包小程序完整版源码 后台API+前端
- 解决小程序开发生成B类小程序码scene参数长度受限的问题
- iPhone5设置铃声方法教程
- kali httrack-网站克隆工具
- 国风频频出圈！品牌如何借势发力？小红书数据查询3招玩转国风
- JS中setter/getter理解