胃癌单细胞数据集GSE163558复现(七)：评估恶性上皮细胞增殖及迁徙才智

图片

Hello小伙伴们全国好，我是生信本领树的小学徒”我才不吃蛋黄“。今天是胃癌单细胞数据集GSE163558复现系列第七期。第六期咱们笔据TCGA数据库中胃癌和往日胃组织之间的各异抒发基因，界说了每个上皮细胞的恶性和非恶性评分。本期，咱们将分析恶性上皮细胞G0-G4的Marker基因并绘画热图和小提琴图，此外，咱们还将使用AddModuleScore_UCell函数揣测细胞的增殖和迁徙评分。

细胞增殖是生物体的紧迫人命特征，细胞以分裂的形貌进行增殖，是生物体滋长、发育、繁衍以及遗传的基础。在肿瘤多阶段演进的早期阶段，肿瘤细胞在原发灶无穷增殖。增殖才智强是肿瘤细胞恶性进程高的紧迫标志之一。除了在原发灶增殖，肿瘤还不错发生搬动，即肿瘤细胞在隔离其发源部位的器官中滋长，搬动是无数肿瘤的最终且最致命的弘扬。本数据集有多个胃癌搬动样本：6名患者的10个簇新东谈主体组织样本，包括3个原发性肿瘤样本(PT)、1个周边非肿瘤样本(NT)和6个搬动样本(M)。搬动样本包括 2个肝脏搬动样本(Li)、2个淋奉迎搬动样本(LN)、1个腹膜搬动样本(P)和1个卵巢搬动样本(O)。肿瘤细胞远方搬动的经由如下：肿瘤细胞侵袭力增强；从原发部位打破血管/淋巴管，投入轮回；从血管再次打破投入组织定植于远方器官，最终在远方器官中增殖。在上述经由中，肿瘤细胞具有不同表型，并在肿瘤微环境中，与其周围免疫细胞和基质细胞互相作用，以支持瘤细胞滋长，并匡助瘤细胞走幸免疫系统的监视。为了更好的打破组织和血管内皮，肿瘤细胞会发生上皮-间充质搬动（Epithelial-Mesenchymal Transition，EMT），即上皮细胞向间充质细胞表型搬动的经由，在EMT发生经由中，肿瘤上皮细胞向间充质细胞搬动，弘扬为细胞时势和功能的双重篡改，时势上由多边形或鹅卵石状转形成细长的纺锤状或梭形，功能上细胞极性消散、细胞骨架篡改、细胞间去粘连化以及赢得侵袭领略才智等；这些表型篡改会使得细胞间黏附度缩短，迁徙领略特质增强。使用AddModuleScore_UCell函数揣测细胞的增殖和迁徙评分,不错协助咱们评估不同恶性上皮亚群的增殖、搬动后劲和恶性进程。

最初断根系统环境变量，修复责任目次，加载R包，读取恶性上皮Seurat数据:

rm(list=ls())getwd()setwd('6-TCGA_STAD/')library(tidyverse)library(tinyarray)library(data.table) library(Seurat)scRNA = readRDS('malignant.rds')

使用Seurat内置函数FindMarkers，以G1为对照，分析恶性上皮细胞各亚群（G0-4）Marker基因：

head(scRNA@meta.data)Idents(scRNA) = scRNA$celltypect = levels(scRNA@active.ident)ct1 = c("G3", "G2", "G4", "G0")all_markers = lapply(ct1, function(x){  # x = ct[1]  print(x)  markers <- FindMarkers(scRNA, group.by = "celltype",                         logfc.threshold = 0.1,                         ident.1 = x,                         ident.2 = ct[5])  #markers_sig <- subset(markers, p_val_adj < 0.1)  return(markers)})

对all_markers的list重定名，然后以“p_val_adj < 0.01”表率筛选各异抒发基因：

length(all_markers)names(all_markers) = ct1lapply(all_markers,nrow)all_markers_sig = lapply(all_markers, function(x){  markers_sig <- subset(x, p_val_adj < 0.01)})

在各异抒发基因的基础上进行富集分析，轮回绘画KEGG和GO迤逦调基因富集条形图：

plot = list()for (i in 1:length(all_markers_sig)){  deg = all_markers_sig[[i]]  deg$change = 'unknown'  deg[deg$avg_log2FC >2,]$change = 'up'  deg[deg$avg_log2FC < -2,]$change = 'down'  table(deg$change)  entrezIDs = bitr(rownames(deg), fromType = "SYMBOL", toType = "ENTREZID", OrgDb= "org.Hs.eg.db", drop = TRUE)  gene<- entrezIDs$ENTREZID  marker_new = deg[rownames(deg) %in% entrezIDs$SYMBOL,]  identical(rownames(marker_new) , entrezIDs$SYMBOL)  p = identical(rownames(marker_new) , entrezIDs$SYMBOL);p  if(!p) entrezIDs = entrezIDs[match(rownames(marker_new) ,entrezIDs$SYMBOL),]  marker_new$ENTREZID = entrezIDs$ENTREZID  a = double_enrich(marker_new,n = 5)  a  plot[[i]] <- a}

plot是包含多个富集分析条形图的list，咱们不错差异索取检察，比如检察G3联系于G1(plot[[1]])各异抒发基因的GO富集分析的条形图（G3_G1$gp）：

G3_G1 = plot[[1]]G3_G1$gp

图片

图片

图片

在这里，咱们取了300个细胞，绘画了'CD63','CLDN4','EGR1'等基因的热图：

Idents(scRNA)scRNA1 = scRNAscRNA1 <-  ScaleData(scRNA1,features = rownames(scRNA1))cells = subset(scRNA1,downsample=300)##取其中的300个细胞，为了图顺眼rownames(cells)gene_order <- c('CD63','CLDN4','EGR1','SRGN','VIM','LAPTM5','AGR2','MT1C','S100A6','PLCG2','SAT1','TSPYL2')gene <- factor(gene_order, levels = unique(gene_order))cells2 <- cells[rownames(cells) %in% gene,]df <- as.data.frame(AverageExpression(object = cells2)$RNA)df <- df[gene_order, ]df = na.omit(df)pheatmap(df, cluster_rows = FALSE, cluster_cols = FALSE,         show_colnames = TRUE, scale = "row",gaps_row = c(seq(3, 11,3)), gaps_col = c(1:4))ggsave('gene_heatmap.pdf',width = 12,height = 8)

图片

数据准备：

最初索取scRNA中基因抒发矩阵raw.data：

raw.data = as.matrix(scRNA@assays$RNA$counts)raw.data[1:6,1:6]length(colnames(raw.data))

索取CD44抒发数据框data.frame

a = scRNA@assays$RNA$datarownames(a)b = as.data.frame(a['CD44',]) colnames(b) = 'CD44'identical(rownames(b),colnames(a))

修复raw.data列名为celltype名，将raw.data赋值为data:

colnames(raw.data) = scRNA$celltypelibrary(ggcorrplot)library(ggthemes)data = raw.datacolnames(data)table(colnames(data))

图片

创建数据框dat：

rownames(data)dat = data.frame(expression = b,group = colnames(data))dat = as.data.frame(dat)

将dat$group修复为因子型变量，从头修复levels规定G0-G4，str函数检察数据框dat变量，na.omit函数删除NA值，数据准备结束后，使用ggplot绘画小提琴图：

dat$group=factor(dat$group, levels = c("G0","G1","G2","G3","G4"))str(dat$group)str(dat$CD44)dat = na.omit(dat)p = ggplot(dat = dat,mapping = aes(x = group,y = CD44)) +geom_violin(scale = "width",adjust =1,trim = TRUE,mapping = aes(fill = group)) +  theme_few() +scale_fill_manual(values = mycolors)+  geom_jitter(width = 0.35, size = 1.1, color = "black") + # 添加点，不错诊治width和size参数  theme(axis.text.x =element_text(size=20), axis.text.y=element_text(size=20))+  labs(x="",y="Expression Level",title = "CD44")+  theme(plot.title = element_text(hjust = 0.5)) +  theme(plot.title = element_text(size=25))+  NoLegend()p

图片

设定参考组，添加显赫性标记（星号）：

library(ggpubr)p+stat_compare_means(method = "anova", label.y = 3.5)+ # Add global p-value  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "G0") # Pairwise comparison against reference

图片

我方设定对比，添加显赫性标记（p值）：

compaired <- list( c("G0", "G1"), c("G1", "G2"), c("G2", "G3") ,c("G3", "G4"))p + stat_compare_means(comparisons=compaired,method = "t.test")

图片

我方设定对比，添加显赫性标记（星号）：

p+geom_signif(comparisons = compaired,step_increase = 0.1,map_signif_level = T,test = t.test)

图片

CD44干系基因抒发烧图(同上):

gene_order <- c('CD44','PROM1' ,'LGR5','SOX2','TFRC','CXCR4' ,'JAG1' )gene <- factor(gene_order, levels = unique(gene_order))cells2 <- cells[rownames(cells) %in% gene,]df <- as.data.frame(AverageExpression(object = cells2)$RNA)df <- df[gene_order, ]df = na.omit(df)pheatmap(df, cluster_rows = FALSE, cluster_cols = FALSE,         show_colnames = TRUE, scale = "row")

图片

最初，笔据原文，界说增殖和迁徙评分list：proliferation = c('MKI67','IGF1','ITGB2','PDGFC','JAG1','PHGDH')；migration = c('VIM','SNAI1','MMP9','AREG','ARID5B' ,'FAT1')，然后使用AddModuleScore_UCell函数揣测评分：

library(UCell)proliferation = c('MKI67','IGF1','ITGB2','PDGFC','JAG1','PHGDH')migration = c('VIM','SNAI1','MMP9','AREG','ARID5B' ,'FAT1')marker <- list(proliferation,migration)#将基因整成listnames(marker)[1] <- 'proliferation'names(marker)[2] <- 'migration'score <- AddModuleScore_UCell(scRNA,                              features=marker)

准备包含增殖/迁徙评分的数据框data（同上），然后绘画小提琴图：

raw.data = as.matrix(score@assays$RNA$counts)raw.data[1:6,1:6]length(colnames(raw.data))rownames(score)a = score$proliferation_UCellcolnames(raw.data) = score$celltypelibrary(ggcorrplot)library(ggthemes)data = raw.datacolnames(data)table(colnames(data))b = data.frame(expression = a,group = colnames(data))data = as.data.frame(b)dat$group=factor(dat$group, levels = c("G0","G1","G2","G3","G4"))str(data$group)str(data$expression)data = na.omit(data)p = ggplot(dat = data,mapping = aes(x = group,y = expression)) +geom_violin(scale = "width",adjust =1,trim = TRUE,mapping = aes(fill = group)) +  theme_few() +scale_fill_manual(values = mycolors)+  geom_jitter(width = 0.35, size = 1.1, color = "black") + # 添加点，不错诊治width和size参数  theme(axis.text.x =element_text(size=20), axis.text.y=element_text(size=20))+  labs(x="",y="Proliferation score")+  theme(plot.title = element_text(hjust = 0.5)) +  theme(plot.title = element_text(size=25))+  NoLegend()p

图片

我方设定对比，然后添加显赫性标记p值：

compaired <- list( c("G0", "G1"), c("G1", "G2"), c("G2", "G3") ,c("G3", "G4"))p + stat_compare_means(comparisons=compaired,method = "t.test")

图片

本期，咱们分析了恶性上皮细胞G0-G4的Marker基因并绘画热图和小提琴图，并使用AddModuleScore_UCell函数揣测细胞的增殖和迁徙评分。下一期，咱们将认真投入单细胞测序高档分析，使用monocle2进行拟时序分析（Pseudo-time analysis）。干货满满，宽容全国抓续追更，谢谢！

图片