HNSCC数据分析-GSE13399-GPL7540 / 开普饭

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

下面是sophie的投稿

数据集介绍

GEO链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE13399
芯片平台：GPL7540, Agilent-scanner-UNC-custom-4X44K-without-Virus [Probe Name version]

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

GSM338295 HNO4-0331B (self-hyb) GSM338296 TONSIL3 (self-hyb) GSM338297 HNO2-0254 (self-hyb) GSM338298 TONSIL4 (self-hyb) GSM338299 HNO4-357B (self-hyb) GSM338300 TONSIL6 (self-hyb) GSM338301 HNO4-604B (self-hyb) GSM338302 TONSIL8 (self-hyb) ...

文章链接是：Standardized WithIn Class Sum of Squares to evaluate methodologies and dataset elements. PLoS One 2010 Mar 26;5(3):e9905. PMID: 20360852(https://www.ncbi.nlm.nih.gov/pubmed/20360852)

核心步骤

R包加载

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

获取并且检查表达量矩阵

主要需要判断是否需要log

# 获取表达量矩阵 gse_number <- 'GSE13399' gset <- geoChina(gse_number) a=gset[[1]] dat=exprs(a) dim(dat)


# 检查，判断需不需要取log

dat[1:4,1:4]

#dat=log2(dat)

boxplot(dat[,1:4],las=2)

library(limma)

dat=normalizeBetweenArrays(dat)

# 画图，使用ggplot需宽数据变长数据 class(dat) data <- as.data.frame(dat) data <- melt(data) head(data) gse_number = "GSE13399" 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_list1=ifelse(grepl('normal',pd$source_name_ch1[1:16]),'control','HNSCC') group_list2=ifelse(grepl('normal',pd$source_name_ch2[17:32]),'control','HNSCC') group_list = factor(c(group_list1,group_list2),levels = c('control','HNSCC'));table(group_list)

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

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

获取芯片注释信息

#代码如下，没有对应R包，需要利用soft文件注释 gpl_number='GPL7540' a = getGEO(gpl_number,destdir = ".") b = a@dataTable@table colnames(b) ids = b[,c("ID","Gene Symbol")]#提取需要的列 colnames(ids) = c("probe_id","symbol")#改成与前边一致的列名 #将ids与前边的对应起来。 ids = ids[ids$symbol!="" & !str_detect(ids$symbol,"///"),] head(ids)> head(ids) 31 AB001451_1_2049 SHC2 32 AB001451_1_2297 SHC2 33 AB002301_1_6436 MAST4 34 AB002301_1_6479 MAST4 35 AB002308_3_6109 SEC16A 36 AB002308_3_6172 SEC16A

探针基因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第一次出现的信息 GSM338295 GSM338296 GSM338297 GSM338298 ZZZ3 4.614219 5.014920 4.803823 3.792137 ZZEF1 6.052270 6.514351 7.237988 6.830712 ZYX 5.146503 6.080845 5.473986 5.356575 ZYG11B 2.817818 4.144349 3.795590 2.743272

以及最简单的2分组，如下所示：

> table(group_list) group_list control HNSCC 16 16

保存为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

ggsave(plot = p2,filename = "p2.png",width = 10,height = 10,units = "cm")

# 下面是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

ggsave(plot = p3,filename = "p3.png",width = 10,height = 10,units = "cm")

# 画图，样品相关性热图p4(不能有含NA的行) colD=data.frame(Group=group_list) exprSet=t(exp) exprSet=na.omit(exprSet) rownames(colD)=colnames(exprSet)#问题-exprSet设置成转置后的exp p4 <- pheatmap::pheatmap(cor(exprSet),#热图对样本-列操作 annotation_col = colD, show_rownames = F, show_colnames = F, main = gse_number ) p4

出图如下:

标准步骤之limma差异分析

配对样本进行配对分析,代码如下：

#差异分析 library(limma) #因为是配对样本，所以要生成一个有效的pairinfo和group，来进行配对分析，注意关注两个参数的设置逻辑~ m=vector() for (i in 1:16) { m=c(m,rep(i,times=2)) } pairinfo = factor(m) group = factor(group_list,levels = c("control","HNSCC")) design <- model.matrix(~group+pairinfo) exprSet=as.data.frame(t(exp)) dim(exprSet) library(tidyr) fit <- lmFit(exprSet,design) fit <- eBayes(fit) options(digits = 4) #设置全局的数字有效位数为4 deg = topTable(fit,coef=2,number=Inf,adjust='fdr')

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

> deg[1:10,] logFC AveExpr t P.Value adj.P.Val B ZBED2 -1.988 8.673 -10.868 1.668e-09 2.905e-05 11.783 IL23A -2.116 7.911 -9.828 8.353e-09 5.019e-05 10.322 PEG3 -2.201 2.519 -9.667 1.083e-08 5.019e-05 10.084 LOC645301|SOX4 2.173 6.523 9.567 1.273e-08 5.019e-05 9.935 RNASE7 -2.069 6.054 -9.357 1.799e-08 5.019e-05 9.616 EGFL6 2.017 7.197 9.331 1.878e-08 5.019e-05 9.576 HOXD11 3.265 3.831 9.246 2.164e-08 5.019e-05 9.445 ZXDA|ZXDB -1.623 5.495 -9.144 2.569e-08 5.019e-05 9.287 SLC22A9 -1.366 4.192 -9.112 2.713e-08 5.019e-05 9.236 SOX4 2.368 6.446 9.027 3.132e-08 5.019e-05 9.102

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

代码如下：

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(dat[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数据分析-GSE13399-GPL7540

数据集介绍

核心步骤

R包加载

获取并且检查表达量矩阵

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

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

获取芯片注释信息

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

标准步骤之质控

标准步骤之limma差异分析

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

symbol 和 ENTREZID转换

进行gsea富集

富集可视化

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

转录组产品线

公共数据库产品线

HNSCC数据分析-GSE13399-GPL7540

数据集介绍

核心步骤

R包加载

获取并且检查表达量矩阵

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

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

获取芯片注释信息

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

标准步骤之质控

标准步骤之limma差异分析

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

symbol 和 ENTREZID转换

进行gsea富集

富集可视化

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

转录组产品线

公共数据库产品线

相关推荐