TNBC数据分析-GSE76275-GPL570

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

下面是sophie的投稿

数据集介绍

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

  • 芯片平台：[GPL570],  Affymetrix GeneChip Human Genome U133 Plus 2.0 Array https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570

样品列表：

共265个，在此就不一一列出了，详见链接：

GSM1974566 S1-H10
GSM1974567 S1-H14
GSM1974568 S1-H19
GSM1974569 S1-H20B
GSM1974570 S1-H22
GSM1974571 S1-H27
GSM1974572 S1-H28
GSM1974573 S1-H29
GSM1974574 S1-H2B
GSM1974575 S1-H31

  • 文章链接是：Comprehensive genomic analysis identifies novel subtypes and targets of triple-negative breast cancer. Clin Cancer Res 2015 Apr 1;21(7):1688-98. PMID: 25208879 https://www.ncbi.nlm.nih.gov/pubmed/25208879 Phosphatase PTP4A3 Promotes Triple-Negative Breast Cancer Growth and Predicts Poor Patient Survival. Cancer Res 2016 Apr 1;76(7):1942-53. PMID: 26921331 https://www.ncbi.nlm.nih.gov/pubmed/26921331

核心步骤

R包加载

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

获取表达量矩阵和分组信息

  • 主要是获取分组信息和判断表达矩阵是否需要log 在读取pd进行样本分组时，发现利用pd任何一列都无法正确区分TNBC和non-TNBC得到文献中给出的分组样本数，但是GEO提供了两种样本的分开版本，所以分别处理 GSE76275的两个子数据集：'GSE76124'和'GSE76274'。下边先处理'GSE76124'，有198个TNBC

获取并且检查表达量矩阵

  • 先获取表达矩阵

gse_number <- 'GSE76275'
#gset <- geoChina(gse_number)
#如果报错：The size of the connection buffer (131072) was not large enough to fit a complete line，则进行如下设置
Sys.setenv("VROOM_CONNECTION_SIZE"=131072*10)
gset <- getGEO(gse_number,
destdir = '.', #指向下载的位置-工作目录
getGPL = F)
a=gset[[1]]
dat=exprs(a)
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_text(angle=90),
plot.title = element_text(hjust = 0.5,size =15))+
xlab('')+
ylab('')+
ggtitle(title)
p1

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

根据生物学背景、研究目的和子数据集进行人为分组

# 1. 获取GSE76124的表达量矩阵
gse_number <- 'GSE76124'
gset <- getGEO(gse_number, 
               destdir = '.', #指向下载的位置-工作目录
               getGPL = F)
m=gset[[1]] 
exp=exprs(m) 
dim(exp)
exp[1:4,1:4] 
paste(gse_number,"exp",sep = "_")
GSE76124_exp= exp

# 2. 获取GSE76274表达量矩阵
gse_number <- 'GSE76274'
gset <- getGEO(gse_number, 
               destdir = '.', #指向下载的位置-工作目录
               getGPL = F)
m=gset[[1]] 
exp=exprs(m) 
dim(exp)
paste(gse_number,"exp",sep = "_")
GSE76274_exp= exp
dim(GSE76274_exp)

#3. 按照子数据集的列名生成分组信息
dim(GSE76124_exp)
dim(GSE76274_exp)
colnames_exp = c(colnames(GSE76124_exp),colnames(GSE76274_exp));length(colnames_exp)
names(group_list)=colnames_dat
table(group_list)
identical(colnames(b),colnames_dat)
#所以表达矩阵的样本是按照先TNBC后non-TNBC排列的.
group_list = as.factor(c(rep("TNBC",198),rep("non-TNBC",67)))
group_list = factor(group_list,levels = c("non-TNBC","TNBC"))

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

probe_id 和symbol的转换至表达矩阵

获取芯片注释信息

library(stringr)
ids=idmap('GPL570')
#超级好用的函数，首选，如果不行再尝试其他

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

> head(ids)
        probe_id symbol
193731   1053_at   RFC2
193732    117_at  HSPA6
193733    121_at   PAX8
193734 1255_g_at GUCA1A
193735   1316_at   THRA
193736   1320_at PTPN21

探针基因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]
table(group_list)
save(dat,pd,ids,group_list,gse_number,gpl_number,file = "step1-output.Rdata")

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

> dat[1:4,1:4]  #保留每个基因ID第一次出现的信息
       GSM1974566 GSM1974567 GSM1974568 GSM1974569
ZZZ3     9.816587   9.610849   9.567354  10.329330
ZZEF1    8.537570   8.836201   8.886380   8.853307
ZYX      9.405403   8.879054   8.072259   8.111867
ZYG11B   9.728289   9.608074  10.726756  10.250213

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

>table(group_list)
group_list
non-TNBC     TNBC 
      67      198 

保存为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,这个矩阵是不含有分组信息的
dim(exp)
exp[,12488]
# 画图，主成分分析图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,
main = gse_number
)
p4

出图如下：

标准步骤之limma差异分析

library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,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
COL18A1-AS1  -1.0422   6.215 -40.81 2.471e-116 4.988e-112 252.9
FZD9          1.5495   7.108  27.14  1.612e-78  1.627e-74 167.8
C5orf49      -0.6688   4.795 -24.36  1.237e-69  8.327e-66 147.7
GDNF         -0.5441   7.131 -23.23  6.998e-66  3.532e-62 139.2
RIBC1        -0.9200   6.953 -23.05  2.701e-65  1.091e-61 137.8
VSIG10L2     -1.0818   6.296 -22.07  5.243e-62  1.764e-58 130.4
SLC25A21-AS1 -0.9174   5.232 -21.93  1.606e-61  4.633e-58 129.3
CCDC170      -1.0647   6.190 -21.34  1.700e-59  4.290e-56 124.6
GTSE1         1.0516   9.499  20.42  2.324e-56  5.213e-53 117.5
BCAS1        -1.8479   6.888 -20.37  3.473e-56  7.011e-53 117.1

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

代码如下：

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 >2,'up',
ifelse( df$logFC < -2,'down','stable') )
)
#统计上下调基因数量
table(df$g)
#给绘制火山图用的数据新增一列symbol
df$name=rownames(df)
head(df)
logFC_t = 2
#设置可循环使用的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
library(ggplot2)
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

出图如下：

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

转录组产品线

公共数据库产品线

赞 (0)

相关推荐

  • [R] 如何绘制各样本的pathway丰度热图？

    前言 一般而言,我们做完pathway富集分析,就做下气泡图或bar图来进行展示,但它们实际上只考虑了富集因子和Pvalue.如果我们不关注这两个因素,而是在乎样本本身的pathway丰度呢? 对于K ...

  • TNBC数据分析-GSE27447-GPL6244

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

  • 三阴性乳腺癌表达数据分析笔记之TNBC定义

    学徒和学员已经陆续出师,是时候把生信技能树的舞台交给后辈了! 下面是学徒写的<GEO数据挖掘课程>的配套笔记(第5篇) B站课程<三阴性乳腺癌表达矩阵探索>笔记之文献解读 三阴 ...

  • 人员数据分析的CRISP-DM模型

    如何证明人力资源实践的有效性是重要且有价值,传统上,研究人员通过使用调查,访谈或观察收集数据来产生此类证据.借助这些数据,他们获得了对劳动力的洞察力,并制定了切实可行的干预措施以改善结果. 技术进步导 ...

  • HR数据分析--员工绩效指标

    员工绩效指标是跟踪员工绩效的关键,正确地实施它们是棘手的.但是,如果做得正确,员工绩效指标将使组织和员工都受益.我们在下面列出了最重要的指标,并提供了每个指标的一些实际示例. 员工绩效指标多种多样.我 ...

  • 为HR数据分析建立业务假设?

    为业务人员分析制定业务问题和发展假设,以确保你在分析主题中增加业务价值,研究如何构建业务问题,业务问题是否与实际定义假设相关. 什么是假设?假设是:基于有限证据做出的假设或建议的解释作为进一步调查的起 ...

  • HR数据分析中常用的21个数据源

    我们通常听到的一个问题是"什么可以用于分析的数据源?" 在本文中,我们将列出HR和更广泛业务中的许多常见数据源,这些数据源将有助于您进行人员分析. HR数据源可以分为3类: 一.H ...

  • 人力资源数据分析

    最近几天,支付宝.抖音.酷狗.喜马拉雅等公司相继发布2019年个人使用报告,发现自己的所作所为都在上面展现的一览无余没有死角,一方面感到数据分析的可怕,另外一方面在想是否可以利用数据在促进工作的提升, ...

  • 人力资源数据分析10条黄金法则

    根据德勤关于全球人力资本趋势的报告,人力资源数据分析革命正在加速.完全有能力应用人力资源分析的组织从4%增至8%.感觉有些能力应用人力资源分析的组织从24%增加到32%.我很高兴看到我的国家(荷兰)在 ...

  • Python数据分析库有哪些？常见分类！

    众所周知,Python前景好.需求量大.薪资高.就业岗位多,除了基本的开发工作之外,还可以从事人工智能.数据分析.网络爬虫等岗位.那么说起数据分析,你知道Python常用数据分析库有哪些吗?我们一起来 ...