HNSCC数据分析-GSE107591-GPL6244 / 开普饭

五月份的学徒专注于GEO数据库里面的表达量芯片数据处理，主要的难点是表达量矩阵获取和探针的基因名字转换，合理的分组后就是标准的差异分析，富集分析。主要是参考我八年前的笔记：

下面是sophie的投稿

数据集介绍

GEO链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE107591
芯片平台：[GPL6244] (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6244), GPL6244[HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version]

样品列表如下(样本太多并未一一列出，详见链接)

GSM2871752 68N GSM2871753 73N GSM2871754 73T GSM2871755 74N GSM2871756 74T GSM2871757 55T ...

文献链接是：

The oncogenic role of circPVT1 in head and neck squamous cell carcinoma is mediated through the mutant p53/YAP/TEAD transcription-competent complex. Genome Biol 2017 Dec 20;18(1):237. PMID: [29262850] (https://www.ncbi.nlm.nih.gov/pubmed/29262850)

TMPRSS2, a SARS-CoV-2 internalization protease is downregulated in head and neck cancer patients. J Exp Clin Cancer Res 2020 Sep 23;39(1):200. PMID: [32967703] (https://www.ncbi.nlm.nih.gov/pubmed/32967703)

核心步骤

R包加载

rm(list = ls()) library(AnnoProbe) library(GEOquery) library(ggplot2) library(ggstatsplot) library(reshape2) library(patchwork) library(stringr)

获取并且检查表达量矩阵

主要需要判断是否需要log

# 获取表达量矩阵 gse_number <- 'GSE107591' gset <- geoChina(gse_number) a=gset[[1]] dat=exprs(a) dim(dat) dat = dat[!duplicated(rownames(dat)),] # 检查，判断需不需要取log dat[1:4,1:4] boxplot(dat[,1:4],las=2) library(limma) dat=normalizeBetweenArrays(dat)

# 画图，使用ggplot需宽数据变长数据 class(dat) data <- as.data.frame(dat) data <- melt(data) head(data) title <- paste (gse_number, "/", a@annotation, sep ="") p1 <- ggplot(data,aes(x=variable,y=value))+ geom_boxplot()+ theme_ggstatsplot()+ theme(panel.grid = element_blank(), axis.text=element_text(size=10,face = 'bold'), axis.text.x=element_blank(), plot.title = element_text(hjust = 0.5,size =15))+ xlab('')+ ylab('')+ ggtitle(title) p1

可以看到，处理前后我们的表达量矩阵的表达量范围箱线图如下所示：

根据生物学背景及研究目的人为分组

pd=pData(a) #通过查看说明书知道取对象a里的临床信息用pData #查看pd发现是配对样本，但分组信息存在于两列中，因此，可以手动生成分组'control','HNSCC' library(stringr) group_list=ifelse(grepl('normal',pd$source_name_ch1),'control','HNSCC') group_list = factor(group_list,levels = c('control','HNSCC'));table(group_list)

为了演示方便，我们这里仅仅是区分"HNSCC"和“control”。

probe_id 和symbol的转换并对应至表达矩阵

获取芯片注释信息

#代码如下 gpl_number='GPL6244' ids = idmap(gpl_number)

可以看到，探针ID和symbol的对应关系如下：

> head(ids) probe_id symbol 448262 7896759 LINC01128 448263 7896761 SAMD11 448264 7896779 KLHL17 448265 7896798 PLEKHN1 448266 7896817 ISG15 448267 7896822 AGRN

探针基因ID对应以及去冗余

代码如下：

library(tidyr) library(dplyr) #接下来，使探针与基因symbol一一对应 ids=as.data.frame(ids) table(rownames(dat) %in% ids$probe_id) dat=dat[rownames(dat) %in% ids$probe_id,] ids=ids[match(rownames(dat),ids$probe_id),] ids$probe_id=as.character(ids$probe_id) rownames(dat)=ids$probe_id ids=ids[ids$probe_id %in% rownames(dat),] dat=dat[ids$probe_id,] #下一步是：基因symbol去冗余-按照表达量最大值筛选 ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行 ids=ids[order(ids$symbol,ids$median,decreasing = T),] #对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s #获得去冗余之后的dat/exp dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat #把ids的symbol这一列中的每一行给dat作为dat的行名 rownames(dat)=ids$symbol dat[1:4,1:4]

最后得到了表达量矩阵如下所示：

> dat[1:4,1:4] #保留每个基因ID第一次出现的信息 GSM2871752 GSM2871753 GSM2871754 GSM2871755 ZZZ3 7.934926 7.691796 7.776917 7.767660 ZZEF1 7.911873 7.824326 7.569077 7.933359 ZYX 8.215491 8.643215 9.206460 8.919280 ZYG11B 8.624613 7.726520 8.073491 7.742189```以及最简单的2分组，如下所示： > table(group_list) group_list control HNSCC 23 24

保存为R数据文件:step1-output.Rdata

标准步骤之质控

得到标准的3张图，包括主成分分析，高变基因的表达量热图，样品相关性热图

代码如下：

## 下面是画PCA的必须操作，需要看说明书。 exp = dat exp=t(exp)#画PCA图时要求是行名是样本名，列名时探针名，因此此时需要转换 exp=as.data.frame(exp)#将matrix转换为data.frame library("FactoMineR")#画主成分分析图需要加载这两个包 library("factoextra") dat.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的


# 画图，主成分分析图p2

this_title <- paste0(gse_number,'_PCA')

p2 <- fviz_pca_ind(dat.pca,

geom.ind = "point", # show points only (nbut not "text")

col.ind = group_list, # color by groups

palette = "Dark2",

addEllipses = TRUE, # Concentration ellipses

legend.title = "Groups")+

ggtitle(this_title)+

theme_ggstatsplot()+

theme(plot.title = element_text(size=15,hjust = 0.5))

p2 # 下面是1000_sd热图 library(pheatmap) cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个 n=t(scale(t(dat[cg,]))) n[n>2]=2 n[n< -2]= -2 n[1:4,1:4] ac=data.frame(Group=group_list) rownames(ac)=colnames(n) # 画图，高变基因的表达量热图p3 p3 <- pheatmap::pheatmap(n, show_colnames =F, show_rownames = F, main = gse_number, annotation_col=ac, breaks = seq(-3,3,length.out = 100))#因为已经手动设置了表达量最大值，所以，可以不用设置break p3 # 画图，样品相关性热图p4 colD=data.frame(Group=group_list) exprSet=t(exp) rownames(colD)=colnames(exprSet)#问题-exprSet设置成转置后的exp p4 <- pheatmap::pheatmap(cor(exprSet),#热图对样本-列操作 annotation_col = colD, show_rownames = F, show_colnames = F, main = gse_number, border_color = NA) p4

出图如下:

标准步骤之limma差异分析

#差异分析 library(limma) design=model.matrix(~factor( group_list )) exprSet=as.data.frame(t(exp)) dim(exprSet) fit=lmFit(exprSet,design) fit=eBayes(fit) options(digits = 4) #设置全局的数字有效位数为4 deg = topTable(fit,coef=2,adjust='BH', n=Inf) #结果检验 boxplot(dat[rownames(deg)[1],]~group_list) deg[1,] #核对无误

差异分析结果前10行如下所示:

> deg[1:10,] logFC AveExpr t P.Value adj.P.Val B SCIN -2.126 6.566 -9.652 7.727e-13 1.456e-08 18.78 ABCA8 -2.589 6.248 -8.253 8.962e-11 4.717e-07 14.33 PPARG -1.096 5.645 -8.189 1.117e-10 4.717e-07 14.12 PTN -1.677 6.727 -8.147 1.294e-10 4.717e-07 13.98 ATP6V0A4 -1.820 6.846 -8.128 1.385e-10 4.717e-07 13.92 SH3BGRL2 -2.549 8.491 -8.104 1.503e-10 4.717e-07 13.84 ANKRD20A11P -2.291 6.746 -7.942 2.643e-10 7.113e-07 13.31 CYP4F35P -2.794 5.491 -7.851 3.635e-10 8.241e-07 13.01 CTTNBP2 -1.808 6.107 -7.821 4.039e-10 8.241e-07 12.91 MMP12 3.475 6.073 7.798 4.375e-10 8.241e-07 12.83

有了差异分析就可以进行标准的可视化，包括火山图和上下调的差异基因热图

代码如下：

nrDEG=deg head(nrDEG) attach(nrDEG) plot(logFC,-log10(P.Value))#简单画图看一下 df=nrDEG df$v= -log10(P.Value) #df新增加一列'v',作为新的绘图参数，值为-log10(P.Value) #设定上下调基因 df$g=ifelse(df$P.Value>0.05,'stable', ifelse( df$logFC >1,'up', ifelse( df$logFC < -1,'down','stable') ) ) #统计上下调基因数量 table(df$g) #给绘制火山图用的数据新增一列symbol df$name=rownames(df) head(df) logFC_t = 1 #设置可循环使用的plot标题 this_tile <- paste0('Cutoff for logFC is ',round(logFC_t,3), '\nThe number of up gene is ',nrow(df[df$g == 'up',]) , '\nThe number of down gene is ',nrow(df[df$g == 'down',]) ) #画图，火山图p5 p5 <- ggplot(data = df, aes(x = logFC, y = -log10(P.Value))) + geom_point(alpha=0.6, size=1.5, aes(color=g)) + ylab("-log10(Pvalue)")+ scale_color_manual(values=c("#34bfb5", "#828586","#ff6633"))+ geom_vline(xintercept= 0,lty=4,col="grey",lwd=0.8) + xlim(-3, 3)+ theme_classic()+ ggtitle(this_tile )+ theme(plot.title = element_text(size=12,hjust = 0.5), legend.title = element_blank(), ) p5

#热图 library(pheatmap) x=deg$logFC names(x)=rownames(deg) cg=c(names(head(sort(x),100)), names(tail(sort(x),100)))#对x进行从小到大排列，取前100及后100，并取其对应的探针名，作为向量赋值给cg n=t(scale(t(exprSet[cg,]))) n[n>2]=2 n[n< -2]= -2 n[1:4,1:4] pheatmap(n,show_colnames =F,show_rownames = F) ac=data.frame(group=group_list) rownames(ac)=colnames(n) # 画图，上下调的差异基因热图p6 p6 <- pheatmap(n,show_colnames =F, show_rownames = F, cluster_cols = T, main = gse_number, annotation_col=ac) p6

出图如下:

标准步骤之生物学功能数据库注释

我们这里不根据任何武断的阈值来区分统计学显著的上下调基因，而是直接根据基因的变化情况排序进行gsea分析，而且仅仅是展示kegg这个生物学功能数据库的注释情况！

symbol 和 ENTREZID转换

gsea分析需要基因的ENTREZID，根据物种进行转换

# 加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join） deg$symbol=rownames(deg) library(org.Hs.eg.db) library(clusterProfiler) s2e <- bitr(deg$symbol, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)#人 #bitr()用于SYMBOL转ENTREZID #其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb dim(deg) dim(s2e) setdiff(deg$symbol,s2e$SYMBOL) DEG <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))

进行gsea富集

library(dplyr) library(ggplot2) geneList=DEG$logFC names(geneList)=DEG$ENTREZID geneList=sort(geneList,decreasing = T) head(geneList) library(clusterProfiler) kk_gse <- gseKEGG(geneList = geneList, organism = 'hsa',#按需替换 #nPerm = 1000, minGSSize = 10, pvalueCutoff = 0.9, verbose = FALSE) tmp=kk_gse@result dim(tmp) kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')#按需替换 #DOSE::setReadable():mapping geneID to gene Symbol tmp=kk@result dim(tmp) pro='comp1' write.csv(kk@result,paste0(pro,'_kegg.gsea.csv')) save(kk,file = 'gsea_kk.Rdata')

富集可视化

上面的kk这个变量就存储了kegg这个生物学功能数据库的gsea分析结果，我们进行简单可视化，代码如下：

# 展现前6个上调通路和6个下调通路 down_k <- kk_gse[tail(order(kk_gse$enrichmentScore,decreasing = F)),];down_k$group=-1 up_k <- kk_gse[head(order(kk_gse$enrichmentScore,decreasing = F)),];up_k$group=1

dat=rbind(up_k,down_k) colnames(dat) dat$pvalue = -log10(dat$pvalue) dat$pvalue=dat$pvalue*dat$group dat=dat[order(dat$pvalue,decreasing = F),] # gsea分析结果p7 p7<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + geom_bar(stat="identity") + scale_fill_gradient(low="#34bfb5",high="#ff6633",guide = FALSE) + scale_x_discrete(name ="Pathway names") + scale_y_continuous(name ="log10P-value") + coord_flip() + theme_ggstatsplot()+ theme(plot.title = element_text(size = 15,hjust = 0.5), axis.text = element_text(size = 12,face = 'bold'), panel.grid = element_blank())+ ggtitle("Pathway Enrichment") p7 #具体看上面条形图里面的每个通路的gsea分布情况p8 library(enrichplot) p8 <- gseaplot2(kk, geneSetID = rownames(down_k))+ gseaplot2(kk, geneSetID = rownames(up_k)) p8

HNSCC数据分析-GSE107591-GPL6244

数据集介绍

核心步骤

R包加载

获取并且检查表达量矩阵

根据生物学背景及研究目的人为分组

probe_id 和symbol的转换并对应至表达矩阵

获取芯片注释信息

探针基因ID对应以及去冗余

标准步骤之质控

标准步骤之limma差异分析

标准步骤之生物学功能数据库注释

symbol 和 ENTREZID转换

进行gsea富集

富集可视化

如果你也有类似的数据分析需求，却苦于不会写代码，可以考虑找我们的工程师帮忙哦！

转录组产品线

公共数据库产品线

HNSCC数据分析-GSE107591-GPL6244

数据集介绍

核心步骤

R包加载

获取并且检查表达量矩阵

根据生物学背景及研究目的人为分组

probe_id 和symbol的转换并对应至表达矩阵

获取芯片注释信息

探针基因ID对应以及去冗余

标准步骤之质控

标准步骤之limma差异分析

标准步骤之生物学功能数据库注释

symbol 和 ENTREZID转换

进行gsea富集

富集可视化

如果你也有类似的数据分析需求，却苦于不会写代码，可以考虑找我们的工程师帮忙哦！

转录组产品线

公共数据库产品线

相关推荐