HNSCC数据分析-GSE6631-HG_U95Av2-GPL8300

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

下面是sophie的投稿

数据集介绍

  • GEO链接:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6631

  • 芯片平台:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL8300 [HG_U95Av2] Affymetrix Human Genome U95 Version 2 Array

样品列表:

GSM153813 Normal mucosa 1
GSM153814 Cancer 1
GSM153815 Normal mucosa 2
GSM153816 Cancer 2
GSM153817 Normal mucosa 3
GSM153818 Cancer 3
GSM153819 Normal mucosa 4
...

  • 文章链接是:Selection and validation of differentially expressed genes in head and neck cancer. Cell Mol Life Sci 2004 Jun;61(11):1372-83. PMID: 15170515 https://www.ncbi.nlm.nih.gov/pubmed/15170515

核心步骤

R包加载

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

获取并且检查表达量矩阵

  • 主要是得是否需要log

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

# 检查,判断需不需要取log
dat[1:4,1:4]
dat = log2(dat)
boxplot(dat[,1:4],las=2)
library(limma)
dat=normalizeBetweenArrays(dat)
dat[1:4,1:4]
# 画图,使用ggplot需宽数据变长数据
class(dat)
data <- as.data.frame(dat)
data <- melt(data)
head(data)
title <- paste (gse_number, "/", gpl_number, 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_text(angle=90),
plot.title = element_text(hjust = 0.5,size =15))+
xlab('')+
ylab('')+
ggtitle(title)
p1
> dat[1:4,1:4] 
           GSM153813 GSM153814 GSM153815 GSM153816
1000_at    8.882447  9.080824  9.004535  8.681835
1001_at    6.815514  6.420102  6.472333  6.402285
1002_f_at  5.381824  5.200093  5.700040  5.723010
1003_s_at  7.219051  7.129450  7.361891  7.504878

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

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

library(stringr)
pd=pData(a)
group_list=ifelse(grepl('Normal',pd$title),'control','HNSCC')
group_list = factor(group_list,levels = c('control','HNSCC'));table(group_list)

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

probe_id 和symbol的转换至表达矩阵

获取芯片注释信息

代码如下:

ids = idmap("GPL8300")
dim(ids)

可以看到此芯片的探针与基因ID或者symbol的对应关系如下所示:

> head(ids)
        probe_id  symbol
553408   1000_at   MAPK3
553409   1001_at    TIE1
553410 1002_f_at CYP2C19
553411 1003_s_at   CXCR5
553412   1004_at   CXCR5
553413   1005_at   DUSP1

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

代码如下:

library(tidyr)
library(dplyr)
library(stringr)
#接下来,使探针与基因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
table(group_list)

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

> dat[1:4,1:4]  #保留每个基因ID第一次出现的信息
      GSM153813 GSM153814 GSM153815 GSM153816
ZZZ3   6.194570  6.398336  5.749128  6.469465
ZZEF1  5.952447  5.816752  5.864333  5.925680
ZYX   10.959871 10.586963 10.338066 10.142952
ZWINT  5.722146  6.759732  6.574198  7.398797

以及最简单的2分组,如下所示:

>table(group_list)
group_list
control   HNSCC 
     22      22

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

save(dat,ids,group_list,gse_number,title,file = "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
)
p4

出图如下:

标准步骤之limma差异分析

代码如下:

library(limma)
#因为是配对样本,所以要生成一个有效的pairinfo和group,来进行配对分析,注意关注两个参数的设置逻辑~
m=vector()
for (i in 1:22) {
m=c(m,rep(i,times=2))
}
pairinfo = factor(m)
design <- model.matrix(~group_list+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
MAL    -3.5154  10.802 -13.802 8.572e-13 7.359e-09 18.88
KRT4   -3.8327  11.516 -12.814 4.060e-12 1.743e-08 17.48
SASH1  -1.1653   9.229 -11.518 3.599e-11 1.030e-07 15.48
TGM3   -3.5633   9.432 -11.011 8.873e-11 1.904e-07 14.64
PLOD3   0.5634   8.631   9.822 8.268e-10 1.420e-06 12.54
ITPR2  -0.6529   7.214  -9.703 1.043e-09 1.493e-06 12.32
KLK13  -2.0788   8.963  -9.510 1.527e-09 1.873e-06 11.96
GPD1L  -1.1166   6.864  -9.233 2.660e-09 2.686e-06 11.43
SPINK5 -2.9586   9.547  -9.205 2.816e-09 2.686e-06 11.38
SPP1    2.9068   8.126   9.146 3.174e-09 2.725e-06 11.26

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

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 >0.5,'up',
ifelse( df$logFC < -0.5,'down','stable') )
)
#统计上下调基因数量
table(df$g)
#给绘制火山图用的数据新增一列symbol
df$name=rownames(df)
head(df)
logFC_t = 0.5
#设置可循环使用的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这个生物学功能数据库的注释情况!

  • 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

出图如下:

image-20210814110731350
 

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

转录组产品线

公共数据库产品线

(0)

相关推荐