GEO挖掘实战一、初步探索数据
2020-09-09 本文已影响0人
小贝学生信
「生信技能树」三阴性乳腺癌表达矩阵探索 系列笔记
GEO挖掘实战一、初步探索数据 - 简书
GEO挖掘实战二、差异分析及富集分析 - 简书
GEO挖掘实战三、GSVA - 简书
GEO挖掘实战四、TNBC相关探索 - 简书
- 示例代码
链接：https://pan.baidu.com/s/1r3QlBXrtGON4wMLbOn27KQ
提取码：ajxa
0、阅读文献
主要复现文献：https://pubmed.ncbi.nlm.nih.gov/30175120/
- 主要研究nonTNBC与TNBC的关系
TNBC，即三阴性乳腺癌是指雌激素受体（ER）、孕激素受体（PR）和人表皮生长因子受体（HER2）均阴性的一种特殊类型乳腺癌。TNBC约占所有乳腺癌的15%，其许多生物学特性和基底细胞样型乳腺癌（Basal-like breast cancer）相似，但两者之间存在某些基因表达谱和免疫表型上的差异，因此亦不能完全等同。
- 获取GSE ID：76275
GEO链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76275
1、下载GEO数据
1.1、设置下载镜像源
#设置一般R包下载镜像源
options()$repos
options("repos"="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
#设置bioconductor包下载镜像源
options()$BioC_mirror
options("repos"="https://mirrors.ustc.edu.cn/bioc/")
1.2、R包与数据下载
if (!require("BiocManager"))
install.packages("BiocManager")
if (!require("GEOquery"))
BiocManager::install("GEOquery")
rawdata <- 'GSE76275_gset.Rdata'
if(!file.exists(rawdata)) {
#判断当前工作目录是否存在rawdata。不存在则下载，并保存
gest <- getGEO("GSE76275", destdir = ".",
AnnotGPL = T, #注释文件，可选
getGPL = T) #平台文件，可选
save(gest, file = rawdata)
}
rm(list = ls())
此外还有其它下载方法，如进入https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76275界面时
1
2、提取表达矩阵与分组信息
load('GSE76275_gset.Rdata')
class(gest) #为一个list
length(gest) #查看list只有一个内容
class(gest[[1]]) #注意list[1]还是一个list
?ExpressionSet
gset <- gest[[1]]
- ExpressionSet contains high-throughput assays and experimental metadata.作为高级对象，内容很多。这里我们关键掌握两个相关函数即可
exprspData，具体如下。
#取表达矩阵 exprs
exprs(gset)[1:4,1:4] #行为探针，列为样本
exp <-exprs(gset)
dim(exp) #265个样本，54675个探针
#取样本信息/实验设计 pData
meta <- pData(gset)
colnames(meta)
head(meta,3)
#通过观察，找到分组信息列
table(meta$characteristics_ch1.1)
meta$characteristics_ch1.1=='triple-negative status: not TN'
#利用ifelse语句生成二分类分组信息(是否为三阴性乳腺癌)
group_list <- ifelse(meta$characteristics_ch1.1=='triple-negative status: not TN',
'noTNBC', 'TNBC')
save(exp, group_list, file = "exp_group.Rdata")
rm(list = ls())
2
如上图，芯片探针表达量值均在1-20左右，可以说明表达量已经经过log处理。在后续差异分析时可采用limma/edgeR，而不能使用Deseq2，因为后者的输入表达矩阵需要为raw counts。
3、数据预分析（PCA、heatmap）
load('exp_group.Rdata')
3.1、PCA （Principal Components Analysis，主成分分析）
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
exp[1:4,1:4]
dat = log(exp+1) #防止为0
t(dat)[1:4,1:4] #转置矩阵：行名为样本，列名为基因
dat.pca <- PCA(t(dat), graph = FALSE)
pca.plot=fviz_pca_ind(dat.pca, geom.ind = "point",
col.ind = group_list,
addEllipses = TRUE,
legend.title = "Groups")
pca.plot #从结果来看，两组区分度挺理想的
3.1
3.2、heatmap
这里我们选取变化（离散程度）最大的前1000个探针做热图
cg <- names(head(sort(apply(exp,1,sd), decreasing = T),1000)) #choose_gene
cm <- exp[cg,] #choose_matrics
cm[1:4,1:4]
cm <- t(scale(t(cm))) #标准化处理，更好观察探针在不同组样本间差异
cm[cm>2]=2; cm[cm< -2]= -2 #避免极端值的影响
cm[1:4,1:4]
group_dat <- data.frame(group=group_list,row.names = colnames(exp))
#热图的特定分组格式
head(group_dat)
3.2-1
library(pheatmap)
pheatmap(cm, show_rownames = F,
show_colnames = F,
annotation_col = group_dat)
3.2-2
如上图两组基本也可以通过分开。但是在左侧里部分noTNBC与TNBC还是混在一起。这是主要因为TNBC主要还是基于生/病理指标而不是基因水平定义的，存在一定差异可以理解。
