# 01-TCGA-GBM count表达数据 ---------------------------------------------------- # 01-1 下载原始count数据 -------------------------------------------------------------- # 从UCSC数据库下载TCGA-GBM队列 countExpr <- read.csv("TCGA-GBM.htseq_counts.csv") # 01-2 ID转换 ------------------------------------------------------------- library(stringr) library(org.Hs.eg.db) ### 提取出SYMBOL信息 g2s <- unique(toTable(org.Hs.egSYMBOL)) ### 提取出ENSEMBL信息 g2e <- unique(toTable(org.Hs.egENSEMBL)) ### 整合两个数据，得到ENSEMBL和SYMBOL的对应关系 s2e <- merge(g2e,g2s,by='gene_id') s2e <- s2e[,-1] ## 去除gene_id列 colnames(countExpr)[1] <- "ensembl_id" # 将ensembl_id列.和.后面的数字去除(截取1-15位) countExpr$ensembl_id <- substr(countExpr$ensembl_id,1,15) countExpr <- merge(s2e,countExpr,by = "ensembl_id") countExpr <- countExpr[,-1] ## 查看是否存在表达量全是0的行 View(countExpr[which(rowSums(countExpr[,2:174]) == 0),]) countExpr <- countExpr[which(rowSums(countExpr[,2:174]) > 0),] ## 查看symbol是否存在重复 symbol_list <- countExpr$symbol symbol_list_duplic <- symbol_list[duplicated(countExpr$symbol)] dupli_exp <- countExpr[which(countExpr$symbol %in% symbol_list_duplic),] ## 根据symbol列去重,取均值 countExpr <- aggregate(countExpr[,2:174], by = list(symbol = countExpr$symbol), mean) rownames(countExpr) <- countExpr[,1] countExpr <- countExpr[,-1] # 01-3 countExpr_tumor ---------------------------------------------------------- countExpr_tumor <- as.data.frame(t(countExpr)) # 添加分组(第14和15位的数字,01-09为肿瘤组织样本) countExpr_tumor$Group <- factor(ifelse(as.numeric(substr(rownames(countExpr_tumor),14,15)) < 10,'Tumor','Normal')) # 查看一下肿瘤和正常样本的数量 list(countExpr_tumor$Group) # 提取肿瘤组织样本 countExpr_tumor <- countExpr_tumor[which(countExpr_tumor$Group == "Tumor"),] # 删除Group列 countExpr_tumor <- countExpr_tumor[,-32966] # 将.替换成- rownames(countExpr_tumor) <- gsub("\\.","-",rownames(countExpr_tumor)) countExpr_tumor <- as.data.frame(t(countExpr_tumor)) # mRNA表达矩阵 library(AnnoProbe) genes <- rownames(countExpr_tumor) ID_type = "SYMBOL" symbolID <- annoGene(genes,ID_type) mrna <- symbolID$SYMBOL[which(symbolID$biotypes == "protein_coding")] # N = 18926 countExpr_tumor <- countExpr_tumor[mrna,] # 18926:168 # 02-TCGA-GBM fpkm表达数据 ---------------------------------------------------- # 02-1 下载原始fpkm数据 -------------------------------------------------------------- # 从UCSC数据库下载TCGA-GBM队列 fpkmExpr <- read.csv("TCGA-GBM.htseq_fpkm.csv") # 02-2 ID转换 ------------------------------------------------------------- library(stringr) library(org.Hs.eg.db) ### 提取出SYMBOL信息 g2s <- unique(toTable(org.Hs.egSYMBOL)) ### 提取出ENSEMBL信息 g2e <- unique(toTable(org.Hs.egENSEMBL)) ### 整合两个数据，得到ENSEMBL和SYMBOL的对应关系 s2e <- merge(g2e,g2s,by='gene_id') s2e <- s2e[,-1] ## 去除gene_id列 colnames(fpkmExpr)[1] <- "ensembl_id" # 将ensembl_id列.和.后面的数字去除(截取1-15位) fpkmExpr$ensembl_id <- substr(fpkmExpr$ensembl_id,1,15) fpkmExpr <- merge(s2e,fpkmExpr,by = "ensembl_id") fpkmExpr <- fpkmExpr[,-1] ## 查看是否存在表达量全是0的行 View(fpkmExpr[which(rowSums(fpkmExpr[,2:174]) == 0),]) fpkmExpr <- fpkmExpr[which(rowSums(fpkmExpr[,2:174]) > 0),] ## 查看symbol是否存在重复 symbol_list <- fpkmExpr$symbol symbol_list_duplic <- symbol_list[duplicated(fpkmExpr$symbol)] dupli_exp <- fpkmExpr[which(fpkmExpr$symbol %in% symbol_list_duplic),] ## 根据symbol列去重,取均值 fpkmExpr <- aggregate(fpkmExpr[,2:174], by = list(symbol = fpkmExpr$symbol), mean) rownames(fpkmExpr) <- fpkmExpr[,1] fpkmExpr <- fpkmExpr[,-1] # 02-3 fpkmExpr_tumor ---------------------------------------------------------- fpkmExpr_tumor <- as.data.frame(t(fpkmExpr)) # 添加分组(第14和15位的数字,01-09为肿瘤组织样本) fpkmExpr_tumor$Group <- factor(ifelse(as.numeric(substr(rownames(fpkmExpr_tumor),14,15)) < 10,'Tumor','Normal')) # 查看一下肿瘤和正常样本的数量 list(fpkmExpr_tumor$Group) # 提取肿瘤组织样本 fpkmExpr_tumor <- fpkmExpr_tumor[which(fpkmExpr_tumor$Group == "Tumor"),] # 删除Group列 fpkmExpr_tumor <- fpkmExpr_tumor[,-32966] # 将.替换成- rownames(fpkmExpr_tumor) <- gsub("\\.","-",rownames(fpkmExpr_tumor)) fpkmExpr_tumor <- as.data.frame(t(fpkmExpr_tumor)) # mRNA表达矩阵 library(AnnoProbe) genes <- rownames(fpkmExpr_tumor) ID_type = "SYMBOL" symbolID <- annoGene(genes,ID_type) mrna <- symbolID$SYMBOL[which(symbolID$biotypes == "protein_coding")] # N = 18926 fpkmExpr_tumor <- fpkmExpr_tumor[mrna,] # 18926:168 # 03-生存信息 ----------------------------------------------------------------- # 从UCSC Xena下载GBM患者的生存信息 survDat <- read.csv("TCGA-GBM.survival.csv") samples <- colnames(fpkmExpr_tumor) survDat <- survDat[which(survDat$sample %in% samples),] survDat <- survDat[,-3] # 共167个样本存在生存信息 samples <- survDat$sample fpkmExpr_tumor <- fpkmExpr_tumor[,samples] # 18926:167 countExpr_tumor <- countExpr_tumor[,samples] # 18926:167 # 04-获取氧化应激基因集 ------------------------------------------------------------ # osRGenes <- read.csv("osGenes.csv") osRGenes <- read.csv("GeneCards-SearchResults.csv") # 仅保留蛋白编码基因 osRGenes <- osRGenes[which(osRGenes$Category == "Protein Coding"),] # 仅保留Relevance.score > 15的基因 osRGenes <- osRGenes[which(osRGenes$Relevance.score > 15),] osRGenes <- osRGenes$Gene.Symbol # N = 168 write.table(osRGenes,file = "osRGenes.xls", quote = F, sep = "\t", row.names = F) # 05-一致性聚类 ---------------------------------------------------------------- library(ConsensusClusterPlus) library(Biobase) library(pheatmap) # 从fpkmExpr_tumor表达矩阵中提取ARGs表达量 fpkmExpr_tumor_OSRGs <- fpkmExpr_tumor[osRGenes,] fpkmExpr_tumor_OSRGs <- na.omit(fpkmExpr_tumor_OSRGs) # 标准化 d <- sweep(fpkmExpr_tumor_OSRGs,1,apply(fpkmExpr_tumor_OSRGs,1,median,na.rm = T)) # 注:列为样本,行为基因 cluster <- ConsensusClusterPlus(d = as.matrix(d), maxK = 5, pItem = 0.8, pFeature = 1, reps = 500, clusterAlg = "pam", distance = "pearson", seed = 1262118388.71279, innerLinkage = "complete", finalLinkage = "complete", corUse = "pairwise.complete.obs", plot = 'pdf', title = "cluster") icl <- calcICL(cluster,title = "cluster",plot = "pdf") # 选择最适k值 Kvec = 2:5 x1 = 0.1; x2 = 0.9 PAC = rep(NA,length(Kvec)) names(PAC) = paste("K=",Kvec,sep="") # from 2 to maxK for(i in Kvec){ M = cluster[[i]]$consensusMatrix Fn = ecdf(M[lower.tri(M)]) PAC[i-1] = Fn(x2) - Fn(x1) } optK = Kvec[which.min(PAC)] optK # k = 2 library(tidyverse) table(cluster[[2]]$consensusClass) Cluster <- cluster[[2]]$consensusClass identical(names(Cluster),colnames(fpkmExpr_tumor_OSRGs)) k2 <- cluster[[2]] k2 <- as.data.frame(k2[["consensusClass"]]) k2$sample <- rownames(k2) colnames(k2)[1] <- "group" k2$group <- gsub("1","Cluster1",k2$group) k2$group <- gsub("2","Cluster2",k2$group) # 06-生存分析(基于一致性聚类结果) ------------------------------------------------------- library("survival") library("survminer") dat1 <- merge(k2,survDat,by = "sample") dat1$OS.time <- as.numeric(dat1$OS.time) # 将生存时间由天变成年 dat1$OS.time <- dat1$OS.time / 365 dat1$OS <- as.numeric(dat1$OS) dat1$group <- factor(dat1$group) # 绘制KM曲线 fit <- survfit(Surv(OS.time,OS) ~ group,data = dat1) Kaplan_Meier <- ggsurvplot(fit, pval = TRUE, conf.int = T, legend.labs = c("cluster1","cluster2"), legend.title = "cluster", risk.table = TRUE, risk.table.col = "strata", linetype = "strata", ggtheme = theme_bw(), palette = c("#00AFBB","#FC4E07"))+ xlab("Time in years") Kaplan_Meier # Kaplan_Meier 7:6 # 07-特异性差异基因的鉴定(基于一致性聚类结果) ----------------------------------------------------------- library(DESeq2) # 去除表达矩阵中的小数,转换成整数 countExpr_tumor1 <- as.data.frame(apply(countExpr_tumor,2,function(x) as.integer(x))) rownames(countExpr_tumor1) <- rownames(countExpr_tumor) # 提取c1和c2簇的样本名 c1samples <- k2$sample[which(k2$group == "Cluster1")] # N = 93 c2samples <- k2$sample[which(k2$group == "Cluster2")] # N = 74 # 样本分组 condition <- factor(ifelse(colnames(countExpr_tumor1) %in% c1samples,'Cluster1','Cluster2')) colData <- data.frame(row.names = colnames(countExpr_tumor1),condition) # 构建dds矩阵 dds <- DESeqDataSetFromMatrix(as.matrix(countExpr_tumor1), colData = colData, design = ~ condition) # 对dds矩阵进行标准化 dds <- DESeq(dds) DEGs <- results(dds) DEGs <- na.omit(DEGs) DEGs <- as.data.frame(DEGs) # 提取差异分析结果 logFC_cutoff <- 1 DEGs$change = as.factor(ifelse(DEGs$padj < 0.05 & abs(DEGs$log2FoldChange) > logFC_cutoff, ifelse(DEGs$log2FoldChange > logFC_cutoff ,'Up','Down'),'no significance')) sig_DEGs <- subset(DEGs,DEGs$padj < 0.05 & abs(DEGs$log2FoldChange) > 1) DEGs_list <- rownames(sig_DEGs) # N = 54 write.table(DEGs,file = "DEGs.xls", quote = F, sep = "\t", row.names = T) # 火山图 library(ggplot2) library(ggthemes) library(RColorBrewer) library(scales) library(ggrepel) library(dplyr) data <- DEGs data$symbol <- rownames(data) ggplot(data,aes(x=log2FoldChange,y=-log10(padj)))+ geom_point(aes(color=change),size=1,alpha=1,na.rm = T)+ geom_hline(yintercept =-log10(0.05),color="#990000", linetype="dashed")+ xlim(-2,2)+ geom_vline(xintercept = -1,color="#990000",linetype="dashed")+ geom_vline(xintercept = 1,color="#990000",linetype="dashed")+ theme_bw(base_size = 14)+ scale_color_manual(values=c("orange","#00B2FF","red"))+ xlab(expression(log[2]("FoldChange")))+ ylab(expression(-log[10]("padj")))+ theme(legend.title = element_blank())+ ggtitle(label = "Volcano Plot", subtitle = "Colored by fold-change direction") # volcano 6:6 # 热图 DEGs_exp <- fpkmExpr_tumor[DEGs_list,] DEGs_exp <- log2(DEGs_exp+1) # 将c1样本移到前面,c2样本移到后面 DEGs_exp <- DEGs_exp[,c(which(colnames(DEGs_exp) %in% c1samples),which(colnames(DEGs_exp) %in% c2samples))] # 将down基因放到上面,up基因放在下面 downGenes <- rownames(sig_DEGs)[which(sig_DEGs$change == "Down")] upGenes <- rownames(sig_DEGs)[which(sig_DEGs$change == "Up")] DEGs_exp <- DEGs_exp[c(which(rownames(DEGs_exp) %in% downGenes),which(rownames(DEGs_exp) %in% upGenes)),] ## 注释信息:列(样本) df <- k2 df <- df[,-2,drop = FALSE] colnames(df)[1] <- 'Group' ## 注释信息:行(基因) df1 <- sig_DEGs[,7,drop = F] ## 画图 library(pheatmap) bk = unique(c(seq(-4,4,length = 100))) # 注释条颜色 ann_color <- list(Group = c(Cluster1 = "#00AFBB",Cluster2 = "#FC4E07")) pheatmap(DEGs_exp, cellwidth = 1, cellheight = 6.5, breaks = bk, annotation_col = df, annotation_row = df1, annotation_colors = ann_color, scale = 'row', fontsize_row = 7, cluster_row = F, cluster_col = F, show_rownames = T, show_colnames = F, border_color = NA, drop_levels = T, color = colorRampPalette(c("navy", "white", "firebrick3"))(100),) # heatmap 6:6 # 08-富集分析 ----------------------------------------------------------------- library('BiocManager') library('AnnotationHub') library('org.Hs.eg.db') library('clusterProfiler') library('GOplot') library('ggplot2') library('stringr') ## 获取DEGs的ENTREZID,即将SYMBOL转换成ENTREZID targetID <- bitr(DEGs_list,fromType = 'SYMBOL',toType = 'ENTREZID','org.Hs.eg.db') # 08-1 GO分析 --------------------------------------------------------------- ALL <- enrichGO(gene = targetID$ENTREZID, keyType = 'ENTREZID', ont = "ALL", OrgDb = org.Hs.eg.db, pAdjustMethod = 'BH', pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = T) head(ALL) write.table(ALL,file = "GO.xls",sep = "\t",quote = F,row.names = T) # 绘图 GO_res <- ALL@result # 1.bp数据整理 bp <- GO_res[which(GO_res$ONTOLOGY == "BP"),] # 按照校正后的p值从低到高排序 bp <- bp[order(bp$p.adjust,decreasing = F),] # 提取前10行 bp <- bp[1:10,] # count值升序排列 newbp <- arrange(bp,bp[,10]) # 2.cc数据整理 cc <- GO_res[which(GO_res$ONTOLOGY == "CC"),] # 按照校正后的p值从低到高排序 cc <- cc[order(cc$p.adjust,decreasing = F),] # 提取前10行 cc <- cc[1:10,] # count值升序排列 newcc <- arrange(cc,cc[,10]) # 3.mf数据整理 mf <- GO_res[which(GO_res$ONTOLOGY == "MF"),] # 按照校正后的p值从低到高排序 mf <- mf[order(mf$p.adjust,decreasing = F),] # 提取前10行 mf <- mf[1:10,] # count值升序排列 newmf <- arrange(mf,mf[,10]) # 整合到一起 go_enrich_df <- data.frame(ID = c(newbp$ID,newcc$ID,newmf$ID), Description = c(newbp$Description, newcc$Description, newmf$Description), GeneNumber = c(newbp$Count,newcc$Count,newmf$Count), type = factor(c(rep("Biological process",10), rep("Cellular component",10),rep("Molecular function",10)), levels = c("Molecular function","Cellular component","Biological process"))) go_enrich_df$Description <- factor(go_enrich_df$Description, levels = go_enrich_df$Description) pdf("DEGs_GO.pdf",height = 6,width = 7) ggplot(go_enrich_df,aes(x = Description,y = GeneNumber,fill = type))+ geom_bar(stat="identity",width = 0.8)+ theme_bw()+ ylab("Gene number")+ xlab(NULL)+ scale_x_discrete(limits = rev(levels(go_enrich_df$Description)))+ theme(axis.text=element_text(face = "plain", color="#000000"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ theme(legend.title=element_blank(),legend.position = "top")+ scale_fill_discrete(name="GO category") dev.off() # 08-2 KEGG分析 ------------------------------------------------------------- KEGG <- enrichKEGG( gene = targetID$ENTREZID, organism = "hsa", keyType = "kegg", pAdjustMethod = 'BH', pvalueCutoff = 0.05, qvalueCutoff = 0.2) KEGG_res <- KEGG@result write.table(KEGG_res,file = "KEGG.xls",sep = "\t",quote = F,row.names = T) # 作图 library(wesanderson) library(magrittr) library(tidyverse) library(clusterProfiler) kegg <- read.delim("KEGG.xls") %>% dplyr::select(2,9,6) %>% arrange(desc(p.adjust)) kegg <- kegg[order(kegg$p.adjust,decreasing = F),] kegg <- kegg[c(1:10),] kegg <- kegg[order(kegg$Count,decreasing = F),] kegg$rank <- c(1:10) kegg$Description <- factor(kegg$Description,levels=kegg$Description,ordered = T) ggplot(kegg,aes(Count,Description))+ geom_point(aes(size = Count,color=-1*log10(p.adjust)))+ scale_color_gradient(low = "green",high = "red")+ labs(color = expression(-log[10](p.adjust)), size="Gene number",x="Gene number",y="Pathway name", title="Pathway enrichment")+ theme(plot.title=element_text(hjust=0.5)) # KEGG 6:8 # 09-PPI网络筛选关键基因 ---------------------------------------------------------- # 输出差异基因列表至本地 write.table(DEGs_list, 'DEGs_list.xls', col.names = T, row.names = T, quote = F, sep="\t") # 选择Degree top10的基因作为关键基因 hubgenes <- c("CSF2","CSF3","LCN2","CCL7","CXCL6", "MMP8","CCR8","TNFSF11","ORM1","IL22RA2") # 绘制degree分布直方图 # 上传文件 degreeDat <- read.csv("degree.csv") # 这是按照degree值从大到小排序好的文件 order <- sort(degreeDat$degree,index.return = TRUE,decreasing = F) degreeDat$hubgene <- factor(degreeDat$hubgene,levels = degreeDat$hubgene[order$ix]) ggplot(data = degreeDat,aes(degree,hubgene))+ geom_bar(stat = "identity", width = 0.8,colour="black",size=0.25,fill="#FC4E07",alpha=0.7)+ geom_text(aes(label = degree), hjust = 1, nudge_x = -0.5)+ theme_minimal()+ ylab("")+ ggtitle("Top 10")+ theme(plot.title = element_text(hjust = 0.05)) # degree_top10 5:5 # 10-肿瘤微环境评分(ESTIMATE) ---------------------------------------------------------- library(utils) library(estimate) library(tidyr) library(ggplot2) library(ggsignif) library(ggpubr) library(RColorBrewer) # 10-1 数据处理 ----------------------------------------------------------------- estimateDat <- fpkmExpr_tumor # 基因在行,样本在列 write.table(estimateDat,'estimateDat.txt', col.names = T, row.names = T, quote = F, sep="\t") # 10-2 ESTIMATE计算 ------------------------------------------------------------- filterCommonGenes(input.f = 'estimateDat.txt', output.f = 'estimateDat.gct', id = 'GeneSymbol') estimateScore('estimateDat.gct', 'estimatePurity.gct') scores <- read.table('estimatePurity.gct', skip = 2, header = T) rownames(scores) <- scores[,1] scores <- as.data.frame(t(scores[,3:ncol(scores)])) ## 输出文件 write.table(scores, '肿瘤微环境评分.xls', col.names = T, row.names = T, quote = F, sep="\t") # 10-3 差异分析 ----------------------------------------------------------------- # 加入分组信息 scores$sample <- rownames(scores) rownames(scores) <- NULL ## 将'sample'列'.'替换成'-' scores$sample <- gsub('\\.', '-', scores$sample) ## 加入分组信息 scores <- merge(k2,scores,by = 'sample') colnames(scores)[2] <- 'Group' rownames(scores) <- scores[,1] scores <- scores[,-1] # 10-3-1 stormalScore差异分析 ------------------------------------------------- scores_stormal <- scores[,c(1,2)] # 绘图 ggplot(data = scores_stormal)+ geom_boxplot(mapping = aes(x = Group,y = StromalScore,colour = Group), alpha = 0.5,size=1.5,width = 0.6)+ geom_jitter(mapping = aes(x = Group,y = StromalScore,colour = Group), alpha = 0.3,size = 3)+ scale_color_manual(limits = c("Cluster1","Cluster2"), values = c("#00AFBB","#FC4E07"))+ geom_signif(mapping = aes(x = Group,y = StromalScore), comparisons = list(c("Cluster1","Cluster2")), map_signif_level = T, # T显示显著性，F显示p value size = 1, # 修改线的粗细 textsize = 4, # 修改显著性标记的大小 test = "t.test")+ # 检验的类型 theme_light()+ theme(legend.position ="none")+ labs(x = "",y="Stromal score")+ theme(plot.title = element_text(size = 15, colour = "black", hjust = 0.5), axis.title.y = element_text(size = 15, # family = "myFont", color = "black", face = "bold", vjust = 1.9, hjust = 0.5, angle = 90), legend.title = element_text(color = "black", # 修改图例的标题 size = 15, face = "bold"), legend.text = element_text(color="black", # 设置图例标签文字 size = 10, face = "bold"), axis.text.x = element_text(size = 13, color = "black", face = "bold", vjust = 0.5, hjust = 0.5, angle = 0), axis.text.y = element_text(size = 13, color = "black", face = "bold", vjust = 0.5, hjust = 0.5, angle = 0) ) # scores_stormal 5:5 # 10-3-2 immuneScore差异分析 ------------------------------------------------- scores_immune <- scores[,c(1,3)] # 绘图 ggplot(data = scores_immune)+ geom_boxplot(mapping = aes(x = Group,y = ImmuneScore,colour = Group), alpha = 0.5,size=1.5,width = 0.6)+ geom_jitter(mapping = aes(x = Group,y = ImmuneScore,colour = Group), alpha = 0.3,size = 3)+ scale_color_manual(limits = c("Cluster1","Cluster2"), values = c("#00AFBB","#FC4E07"))+ geom_signif(mapping = aes(x = Group,y = ImmuneScore), comparisons = list(c("Cluster1","Cluster2")), map_signif_level = T, # T显示显著性，F显示p value size = 1, # 修改线的粗细 textsize = 4, # 修改显著性标记的大小 test = "t.test")+ # 检验的类型 theme_light()+ theme(legend.position ="none")+ labs(x = "",y="Immune score")+ theme(plot.title = element_text(size = 15, colour = "black", hjust = 0.5), axis.title.y = element_text(size = 15, # family = "myFont", color = "black", face = "bold", vjust = 1.9, hjust = 0.5, angle = 90), legend.title = element_text(color = "black", # 修改图例的标题 size = 15, face = "bold"), legend.text = element_text(color="black", # 设置图例标签文字 size = 10, face = "bold"), axis.text.x = element_text(size = 13, color = "black", face = "bold", vjust = 0.5, hjust = 0.5, angle = 0), axis.text.y = element_text(size = 13, color = "black", face = "bold", vjust = 0.5, hjust = 0.5, angle = 0) ) # scores_immune # 10-3-3 ESTIMATEScore差异分析 ------------------------------------------------- scores_estimate <- scores[,c(1,4)] # 绘图 ggplot(data = scores_estimate)+ geom_boxplot(mapping = aes(x = Group,y = ESTIMATEScore,colour = Group), alpha = 0.5,size=1.5,width = 0.6)+ geom_jitter(mapping = aes(x = Group,y = ESTIMATEScore,colour = Group), alpha = 0.3,size = 3)+ scale_color_manual(limits = c("Cluster1","Cluster2"), values = c("#00AFBB","#FC4E07"))+ geom_signif(mapping = aes(x = Group,y = ESTIMATEScore), comparisons = list(c("Cluster1","Cluster2")), map_signif_level = T, # T显示显著性，F显示p value size = 1, # 修改线的粗细 textsize = 4, # 修改显著性标记的大小 test = "t.test")+ # 检验的类型 theme_light()+ theme(legend.position ="none")+ labs(x = "",y="ESTIMATE score")+ theme(plot.title = element_text(size = 15, colour = "black", hjust = 0.5), axis.title.y = element_text(size = 15, # family = "myFont", color = "black", face = "bold", vjust = 1.9, hjust = 0.5, angle = 90), legend.title = element_text(color = "black", # 修改图例的标题 size = 15, face = "bold"), legend.text = element_text(color="black", # 设置图例标签文字 size = 10, face = "bold"), axis.text.x = element_text(size = 13, color = "black", face = "bold", vjust = 0.5, hjust = 0.5, angle = 0), axis.text.y = element_text(size = 13, color = "black", face = "bold", vjust = 0.5, hjust = 0.5, angle = 0) ) # scores_estimate 5:5 # 10-3-4 tumorpurity差异分析 ------------------------------------------------------------- scores_tumorpurity <- scores[,c(1,5)] # 绘图 ggplot(data = scores_tumorpurity)+ geom_boxplot(mapping = aes(x = Group,y = TumorPurity,colour = Group), alpha = 0.5,size=1.5,width = 0.6)+ geom_jitter(mapping = aes(x = Group,y = TumorPurity,colour = Group), alpha = 0.3,size = 3)+ scale_color_manual(limits = c("Cluster1","Cluster2"), values = c("#00AFBB","#FC4E07"))+ geom_signif(mapping = aes(x = Group,y = TumorPurity), comparisons = list(c("Cluster1","Cluster2")), map_signif_level = T, # T显示显著性，F显示p value size = 1, # 修改线的粗细 textsize = 4, # 修改显著性标记的大小 test = "t.test")+ # 检验的类型 theme_light()+ theme(legend.position ="none")+ labs(x = "",y="Tumor Purity")+ theme(plot.title = element_text(size = 15, colour = "black", hjust = 0.5), axis.title.y = element_text(size = 15, # family = "myFont", color = "black", face = "bold", vjust = 1.9, hjust = 0.5, angle = 90), legend.title = element_text(color = "black", # 修改图例的标题 size = 15, face = "bold"), legend.text = element_text(color="black", # 设置图例标签文字 size = 10, face = "bold"), axis.text.x = element_text(size = 13, color = "black", face = "bold", vjust = 0.5, hjust = 0.5, angle = 0), axis.text.y = element_text(size = 13, color = "black", face = "bold", vjust = 0.5, hjust = 0.5, angle = 0) ) # tumorpurity 5:5 # 11-关键基因表达量与临床特征的相关性分析 --------------------------------------------------- # 11-1 TCGA-GBM临床信息 ----------------------------------------------------- # 从UCSC Xena数据库下载 clinical <- read.csv("TCGA-GBM.GDC_phenotype.csv") # 选择3个临床指标 clinical <- clinical[,c("submitter_id.samples","gender.demographic","age_at_initial_pathologic_diagnosis","karnofsky_performance_score")] colnames(clinical)[1] <- "sample" # 从fpkm表达矩阵中提取具有临床指标的样本 dat4 <- as.data.frame(t(fpkmExpr_tumor)) dat4$sample <- rownames(dat4) # 将sample列移到前面 dat4 <- dat4[,c(18927,1:18926)] # 将表达矩阵与临床指标进行合并 dat4 <- merge(clinical,dat4,by = "sample") clinical <- dat4[,c(1:4)] dat4 <- dat4[,-c(2:4)] rownames(dat4) <- dat4$sample dat4 <- dat4[,-1] # 删除缺少临床信息的样本 clinical <- na.omit(clinical) samples <- clinical$sample # 共124个样本具完整的临床信息 dat4 <- dat4[samples,] colnames(clinical) <- c("sample","Gender","Age","KPS") clinical$Gender <- gsub("male","Male",clinical$Gender) clinical$Gender <- gsub("feMale","Female",clinical$Gender) clinical$Age <- ifelse(clinical$Age >= 50,"≥50","<50") clinical$KPS <- ifelse(clinical$KPS >= 70,"≥70","<70") # 从dat4中提取10个关键基因的表达量 dat4 <- dat4[,hubgenes] # dat5 <- log2(dat4+1) dat5 <- as.data.frame(t(dat4)) dat5$Genes <- rownames(dat5) # 将宽数据转换成长数据 library(tidyr) dat5 <- gather(data = dat5,key = sample,value = expr,-c("Genes")) # 11-2 Gender ---------------------------------------------------------------- clinical_Gender <- clinical[,c(1,2)] malesamples <- clinical_Gender$sample[which(clinical_Gender$Gender == "Male")] clinical_Gender <- dat5 clinical_Gender$Gender <- ifelse(clinical_Gender$sample %in% malesamples,"Male","Female") # 绘图 library(ggplot2) library(tidyr) library(ggpubr) library(RColorBrewer) ggplot(clinical_Gender,aes(Genes,expr,fill = Gender)) + geom_boxplot(outlier.shape = 21,color = "black") + theme_bw()+ ylim(0,6)+ guides(fill = guide_legend())+ labs(x = NULL, y = "Expression")+ theme(legend.position = "top")+ theme(axis.text.x = element_text(angle = 45,hjust = 1))+ scale_fill_manual(values = c("red","skyblue2"))+ stat_compare_means(data = clinical_Gender, mapping = aes(group = Gender), label ="p.format", method = "t.test", hide.ns = F) # clinical_Gender 6:9 # 11-3 Age ---------------------------------------------------------------- clinical_Age <- clinical[,c(1,3)] samples <- clinical_Age$sample[which(clinical_Age$Age == "≥50")] clinical_Age <- dat5 clinical_Age$Age <- ifelse(clinical_Age$sample %in% samples,"≥50","<50") # 绘图 library(ggplot2) library(tidyr) library(ggpubr) library(RColorBrewer) ggplot(clinical_Age,aes(Genes,expr,fill = Age)) + geom_boxplot(outlier.shape = 21,color = "black") + theme_bw()+ ylim(0,6)+ guides(fill = guide_legend())+ labs(x = NULL, y = "Expression")+ theme(legend.position = "top")+ theme(axis.text.x = element_text(angle = 45,hjust = 1))+ scale_fill_manual(values = c("red","skyblue2"))+ stat_compare_means(data = clinical_Age, mapping = aes(group = Age), label ="p.format", method = "t.test", hide.ns = F) # clinical_Age 6:9 # 11-4 KPS ---------------------------------------------------------------- clinical_KPS <- clinical[,c(1,4)] samples <- clinical_KPS$sample[which(clinical_KPS$KPS == "≥70")] clinical_KPS <- dat5 clinical_KPS$KPS <- ifelse(clinical_KPS$sample %in% samples,"≥70","<70") # 绘图 library(ggplot2) library(tidyr) library(ggpubr) library(RColorBrewer) ggplot(clinical_KPS,aes(Genes,expr,fill = KPS)) + geom_boxplot(outlier.shape = 21,color = "black") + theme_bw()+ ylim(0,6)+ guides(fill = guide_legend())+ labs(x = NULL, y = "Expression")+ theme(legend.position = "top")+ theme(axis.text.x = element_text(angle = 45,hjust = 1))+ scale_fill_manual(values = c("red","skyblue2"))+ stat_compare_means(data = clinical_KPS, mapping = aes(group = KPS), label ="p.format", method = "t.test", hide.ns = F) # clinical_KPS 6:9 # 12-免疫检查点分子与关键基因表达相关性分析 --------------------------------------------------------- # 1.准备免疫检查点分子表达文件 ICMs <- read.table("/data/nas1/zhiyq/pipeline/CheckPoint/checkpoint.txt") ICMs <- ICMs$V1 icmsDat <- fpkmExpr_tumor[ICMs,] # 缺少2个基因,共45个免疫检查点记忆 icmsDat <- na.omit(icmsDat) icmsDat <- log2(icmsDat+1) icmsDat <- as.data.frame(t(icmsDat)) # 2.准备关键基因表达文件 hubgenesDat <- fpkmExpr_tumor[hubgenes,] hubgenesDat <- log2(hubgenesDat+1) hubgenesDat <- as.data.frame(t(hubgenesDat)) # 计算相关性和显著性p值 library(psych) library(ggplot2) library(aplot) library(ggpubr) library(ggrepel) x <- icmsDat y <- hubgenesDat d <- corr.test(x,y,use = "complete",method = 'spearman') r <- d$r p <- d$p cordat <- as.data.frame(r) pdat <- as.data.frame(p) write.table(cordat,file = "cor_hubgene2icms.xls", quote = F, sep = "\t", row.names = T) write.table(pdat,file = "pvalue_hubgene2icms.xls", quote = F, sep = "\t", row.names = T) # 绘制相关性热图 library(ggplot2) library(RColorBrewer) unloadNamespace("WGCNA") library(Hmisc) library(corrplot) library(ggcorrplot) library(ggpmisc) library(pheatmap) pheatmap(cordat, cellwidth = 27, cellheight = 9.5, scale = 'none', fontsize_row = 8, cluster_row = F, cluster_col = F, show_rownames = T, show_colnames = T, fontsize_number = 5, border_color = "white", drop_levels = T, fondsize_row = 6, fondsize_col = 6, display_numbers = TRUE, number_color = "black", color = colorRampPalette(c("blue","cyan","white","yellow","red"))(12),) # corrplot 7:7 # 13-关键基因与免疫因子相关性分析 --------------------------------------------------------------- # 从TISIDB数据库下载免疫因子 immune_factor <- read.csv("immune_factor.csv") # 13-1 chemokine与关键基因相关性 -------------------------------------------------- chemokines <- immune_factor$immune_factor[which(immune_factor$type == "chemokine")] # N = 41 chemokines_expr <- fpkmExpr_tumor[chemokines,] chemokines_expr <- na.omit(chemokines_expr) chemokines_expr <- as.data.frame(t(chemokines_expr)) # 计算相关性和显著性p值 library(psych) library(ggplot2) library(aplot) library(ggpubr) library(ggrepel) x <- chemokines_expr y <- hubgenesDat d <- corr.test(x,y,method = 'spearman') r <- d$r # 相关性 p <- d$p # 显著性p-value write.table(r,file = "cor_chemokines.xls", quote = F, sep = "\t", row.names = T) write.table(p,file = "pvalue_chemokines.xls", quote = F, sep = "\t", row.names = T) # 写一个函数 myfun <- function(pval) { stars = "" if(pval <= 0.001) stars = "***" if(pval > 0.001 & pval <= 0.01) stars = "**" if(pval > 0.01 & pval <= 0.05) stars = "*" if(pval > 0.05 & pval <= 0.1) stars = "" stars } # 准备输入文件 heatmap <- melt(r) %>% rename(replace=c("X1"="chemokines","X2"="hubgenes","value"="cor")) %>% mutate(pvalue=melt(p)[,3]) %>% mutate(signif = sapply(pvalue, function(x) myfun(x))) # 绘图 ggplot(heatmap,aes(chemokines,hubgenes))+ geom_tile(aes(fill = cor))+ geom_text(aes(label = signif),color = "black",size = 2)+ scale_fill_gradient2(low = '#3C8DAD',high = '#FF6767',mid = 'white', limit = c(-1,1), name = paste0("* p < 0.05","\n\n","** p < 0.01","\n\n","*** p < 0.001","\n\n","Correlation"))+ labs(x=NULL,y=NULL)+ theme(axis.text.x = element_text(size = 7,angle = 90,hjust = 1,color = "black"), axis.text.y = element_text(size = 7,color = "black"), axis.ticks.y = element_blank(), panel.background = element_blank()) # cor_chemokines 5:8 # 13-2 immunosuppressive factors与关键基因相关性 -------------------------------------------------- immunosuppressive <- immune_factor$immune_factor[which(immune_factor$type == "Immunoinhibitor")] # N = 24 immunosuppressive_expr <- fpkmExpr_tumor[immunosuppressive,] immunosuppressive_expr <- na.omit(immunosuppressive_expr) immunosuppressive_expr <- as.data.frame(t(immunosuppressive_expr)) # 计算相关性和显著性p值 library(psych) library(ggplot2) library(aplot) library(ggpubr) library(ggrepel) x <- immunosuppressive_expr y <- hubgenesDat d <- corr.test(x,y,method = 'spearman') r <- d$r # 相关性 p <- d$p # 显著性p-value write.table(r,file = "cor_immunosuppressive.xls", quote = F, sep = "\t", row.names = T) write.table(p,file = "pvalue_immunosuppressive.xls", quote = F, sep = "\t", row.names = T) # 写一个函数 myfun <- function(pval) { stars = "" if(pval <= 0.001) stars = "***" if(pval > 0.001 & pval <= 0.01) stars = "**" if(pval > 0.01 & pval <= 0.05) stars = "*" if(pval > 0.05 & pval <= 0.1) stars = "" stars } # 准备输入文件 heatmap <- melt(r) %>% rename(replace=c("X1"="immunosuppressive","X2"="hubgenes","value"="cor")) %>% mutate(pvalue=melt(p)[,3]) %>% mutate(signif = sapply(pvalue, function(x) myfun(x))) # 绘图 ggplot(heatmap,aes(immunosuppressive,hubgenes))+ geom_tile(aes(fill = cor))+ geom_text(aes(label = signif),color = "black",size = 2)+ scale_fill_gradient2(low = '#3C8DAD',high = '#FF6767',mid = 'white', limit = c(-1,1), name = paste0("* p < 0.05","\n\n","** p < 0.01","\n\n","*** p < 0.001","\n\n","Correlation"))+ labs(x=NULL,y=NULL)+ theme(axis.text.x = element_text(size = 7,angle = 90,hjust = 1,color = "black"), axis.text.y = element_text(size = 7,color = "black"), axis.ticks.y = element_blank(), panel.background = element_blank()) # cor_immunosuppressive 5:8 # 13-3 immunostimulatory factors与关键基因相关性 -------------------------------------------------- immunostimulatory <- immune_factor$immune_factor[which(immune_factor$type == "Immunostimulator")] # N = 46 immunostimulatory_expr <- fpkmExpr_tumor[immunostimulatory,] immunostimulatory_expr <- na.omit(immunostimulatory_expr) immunostimulatory_expr <- as.data.frame(t(immunostimulatory_expr)) # 计算相关性和显著性p值 library(psych) library(ggplot2) library(aplot) library(ggpubr) library(ggrepel) x <- immunostimulatory_expr y <- hubgenesDat d <- corr.test(x,y,method = 'spearman') r <- d$r # 相关性 p <- d$p # 显著性p-value write.table(r,file = "cor_immunostimulatory.xls", quote = F, sep = "\t", row.names = T) write.table(p,file = "pvalue_immunostimulatory.xls", quote = F, sep = "\t", row.names = T) # 写一个函数 myfun <- function(pval) { stars = "" if(pval <= 0.001) stars = "***" if(pval > 0.001 & pval <= 0.01) stars = "**" if(pval > 0.01 & pval <= 0.05) stars = "*" if(pval > 0.05 & pval <= 0.1) stars = "" stars } # 准备输入文件 heatmap <- melt(r) %>% rename(replace=c("X1"="immunostimulatory","X2"="hubgenes","value"="cor")) %>% mutate(pvalue=melt(p)[,3]) %>% mutate(signif = sapply(pvalue, function(x) myfun(x))) # 绘图 ggplot(heatmap,aes(immunostimulatory,hubgenes))+ geom_tile(aes(fill = cor))+ geom_text(aes(label = signif),color = "black",size = 2)+ scale_fill_gradient2(low = '#3C8DAD',high = '#FF6767',mid = 'white', limit = c(-1,1), name = paste0("* p < 0.05","\n\n","** p < 0.01","\n\n","*** p < 0.001","\n\n","Correlation"))+ labs(x=NULL,y=NULL)+ theme(axis.text.x = element_text(size = 7,angle = 90,hjust = 1,color = "black"), axis.text.y = element_text(size = 7,color = "black"), axis.ticks.y = element_blank(), panel.background = element_blank()) # cor_immunostimulatory 5:8 # 13-4 immune receptors与关键基因相关性 -------------------------------------------------- receptors <- immune_factor$immune_factor[which(immune_factor$type == "receptor")] # N = 18 receptors_expr <- fpkmExpr_tumor[receptors,] receptors_expr <- na.omit(receptors_expr) receptors_expr <- as.data.frame(t(receptors_expr)) # 计算相关性和显著性p值 library(psych) library(ggplot2) library(aplot) library(ggpubr) library(ggrepel) x <- receptors_expr y <- hubgenesDat d <- corr.test(x,y,method = 'spearman') r <- d$r # 相关性 p <- d$p # 显著性p-value write.table(r,file = "cor_receptors.xls", quote = F, sep = "\t", row.names = T) write.table(p,file = "pvalue_receptors.xls", quote = F, sep = "\t", row.names = T) # 写一个函数 myfun <- function(pval) { stars = "" if(pval <= 0.001) stars = "***" if(pval > 0.001 & pval <= 0.01) stars = "**" if(pval > 0.01 & pval <= 0.05) stars = "*" if(pval > 0.05 & pval <= 0.1) stars = "" stars } # 准备输入文件 heatmap <- melt(r) %>% rename(replace=c("X1"="receptors","X2"="hubgenes","value"="cor")) %>% mutate(pvalue=melt(p)[,3]) %>% mutate(signif = sapply(pvalue, function(x) myfun(x))) # 绘图 ggplot(heatmap,aes(receptors,hubgenes))+ geom_tile(aes(fill = cor))+ geom_text(aes(label = signif),color = "black",size = 2)+ scale_fill_gradient2(low = '#3C8DAD',high = '#FF6767',mid = 'white', limit = c(-1,1), name = paste0("* p < 0.05","\n\n","** p < 0.01","\n\n","*** p < 0.001","\n\n","Correlation"))+ labs(x=NULL,y=NULL)+ theme(axis.text.x = element_text(size = 7,angle = 90,hjust = 1,color = "black"), axis.text.y = element_text(size = 7,color = "black"), axis.ticks.y = element_blank(), panel.background = element_blank()) # cor_receptors 5:8 write.table(countExpr_tumor,file = "countExpr.xls", quote = F, sep = "\t", row.names = T) write.table(fpkmExpr_tumor,file = "fpkmExpr.xls", quote = F, sep = "\t", row.names = T)