保姆式GEO数据挖掘演示--重现9分文章
不知不觉GEO教程系列走过了3讲，大家的反馈是部分代码太干，干货大家也不稀罕，唉！！！ 不过，付费了就是大爷，我们还是得伺候着，所以再次把所有代码标上注释！！！这次应该是都能看懂了。
写在前面
模拟1000行代码不如实操训练，重现文章中的数据才是学习GEO数据挖掘的最好途径，基于以上精神，我们就来重现一下高分文章的数据挖掘过程。
至于为什么选择这篇文章，是因为我还是个GEO数据挖掘的小白，这篇文章在生信技能树旗下公众号“生信菜鸟团”有发布教程（GEO数据挖掘-第二期-三阴性乳腺癌（TNBC））。但即使是如此，对于特级小白的我，还需要一步一步的去理解和构建对代码的感觉，又作为一个科研人员，只有代码训练，没有学习到科研思路也算是憾事，因此本教程结合文献导读和GEO数据挖掘保姆式教程（几句每一句代码的前世今生和重要节点都会解释），给我自己一个完整的从idea提出到数据获取一条龙洗髓Buff，以期冲击一下我的任督二脉。不喜欢看文献导读的可以直接跳转到 4.GEO数据挖掘过程。
原文链接
https://europepmc.org/articles/pmc4873402
1. 摘要
三阴性乳腺癌（Triple-negative breast cancer，TNBC）是乳腺癌中预后最差的一类，且目前并无靶向治疗药物。为了探索TNBC靶向基因，检测乳腺肿瘤中磷酸酶的表达并判断其对in vitro 及 in vivo 中TNBC生长的影响。使用102例乳腺癌的Affymetrix microarray analysis数据，与雌激素受体阳性肿瘤（estrogen receptor (ER)-positive tumors）对比，我们发现146个在TNBC中差异表达的磷酸酶。其中，与ER-positive肿瘤对比，TNBC中有19个上调磷酸酶（0.66-fold; FDR=0.05）。使用小干扰RNA技术（small-interfering RNA, siRNA）下调TNBC细胞系及ER-positive乳腺癌细胞系17个过表达的磷酸酶，发现6种磷酸酶的缺失抑制TNBC生长的幅度大于ER-positive乳腺癌细胞系。phosphatase PTP4A3（也被称为PRL-3）是一个调控细胞周期G1期向S期进展的蛋白酪氨酸激酶，但在这里，phosphatase PTP4A3还选择性调控TNBC细胞凋亡。in vivo实验结果提示PTP4A3抑制剂可减弱TNBC生长。in silico分析结果提示PTP4A3可单独作为一个预测TNBC不良预后的因子，并为TNBC患者针对磷酸激酶或者其各自的信号通路开发新的靶向治疗策略奠定基础。
2. 知识积累
TNBC：三阴性乳腺癌是指不表达雌激素受体（estrogen receptor，ER）、孕激素受体（progesterone receptor，PR）及人表皮生长因子受体2HER-2/neu（Receptor tyrosine-protein kinase erBb-2，HER2/neu）的乳腺癌。
第13届St.Gallen国际乳腺癌会议共识，乳腺癌分型
乳腺癌病理分型（资料参考：National Comprehensive Cancer Network，NCCN）：Ductal（导管型）, includes medullary and micropapillary sybtypes（包括髓样和微乳头糖浆型）；Lobular（小叶型）；Mixed（混合型）；Metaplastic；tubular（管状）；mucinous（粘液性）；
中国乳腺癌指南，乳腺癌不同分子分型推荐治疗如下：
可见在ER/PR positive乳腺癌患者中，至少可以使用内分泌治疗，HER2阳性患者还可以使用赫赛汀(注射用曲妥珠单抗)等抗HER2靶向治疗，而三阴性乳腺癌目前只有常规化疗（环磷酰胺、多柔比星、多西他赛等）可以选择，本文正是基于这个背景，开展三阴性乳腺癌的潜在靶向治疗基因筛选。
以上资料参考来源：http://oncol.dxy.cn/article/534043
PTP4A3：蛋白酪氨酸磷酸酶IVA 3（Protein tyrosine phosphatase type IVA 3，PTP4A3），编码该蛋白的基因属于戊基化蛋白酪氨酸磷酸酶(PTPs)家族，而PTPs在cellular process（细胞周期等）具有重要作用。
in silico analysis：silico是指“硅”，让我想起《三体》中的碳基生物（人类）和硅基生物（某些外星人），in silico其实也可以这么理解，即脱离我们常规的湿实验（in vitro，in vivo，ex vivo等），用计算机模拟或者芯片等等进行的虚拟实验。
Ki-67：Antigen KI-67 (Ki-67)，是一个在细胞生长过程中起重要作用的蛋白（已被视为增殖的biomarker之一），与核糖体RNA转录有关，Ki-67的缺失可以阻滞核糖体RNA的合成。同时，Ki-67在细胞周期的G1、S（显著高表达）、G2和有丝分裂期出现，在G0期中无。因此，目前Ki-67已经被临床病理科用作为乳腺癌分子标记物之一。
PCA：主成分分析（principal component analysis，PCA），使用正交变换将一组可能相关变量(每个变量都具有不同的数值)的观测值转换为一组线性不相关变量的值，称为主成分（理论是不是看不懂，不用懂，知道拿来干啥就行）。即考察多个变量之间的相关性（分组是不是有效），具体想看PCA可以戳链接：一文看懂主成分分析，或者去看YouTube上看看StatQuest系列视频节目的PCA解读~敲喜欢这位小哥，开头都会唱歌哦~
3. 前言
TNBC（无ER，RP及HER2表达）是预后最差且无靶向治疗药物的乳腺癌，表达ER，RP或者HER2的乳腺癌可选择内分泌治疗或抗HER2等靶向治疗。除了ER和HER2信号通路，还有其他信号通路可调控肿瘤生长。我们主要的目的就是找到调控TNBCs生长的关键信号通路和分子。在我们的前期研究中发现，与ER-positive乳腺癌患者对比，ER-negative中过表达的激酶集合，且集合内许多激酶正是ER-negative乳腺癌细胞系生长所必须的（原文这里提到的是 a set of kinases overexpressed in ER-negative compared to ER-positive breast cancers）。磷酸化对于信号通路的调控非常重要。
在本研究中，我们假定磷酸酶在TNBC中的表达与ER-positive肿瘤组是有差别的。具体研究步骤如下：
使用转录组学分析数据（transcriptional profiling data）检测磷酸酶表达水平并找到高表达和低表达组合（TNBC对比ER-positive乳腺癌组）；
发现对TNBC生长至关重要的磷酸酶，并挑选一部分进行验证试验： PTP4A3, PPAP2B, CDC25B, and TIMM50；
PTP4A3在in vitro，in vivo实验中均被证实与TNBC生长相关，PTP4A3抑制剂可增加乳腺癌细胞G1期比例，从而减少细胞生长及小鼠移植瘤中Ki-67的表达；
计算PTP4A3对乳腺癌患者survival的影响。
由此，我们分析出本文中主要用到的GEO数据挖掘结果就是差异分析后得到的上调、下调基因集，以及生存分析。
4. GEO数据挖掘过程
Title: Phosphatase PTP4A3 promotes triple-negative breast cancer growth and predicts poor patient survival
GSE number: GSE76275
sample: 265 samples which included 67 non-TNBC and 198 TNBC
代码来源：生信菜鸟团 GEO数据挖掘-第二期-三阴性乳腺癌（TNBC）
总的思路：1.提取表达矩阵，按照我们所需的分组排序得到矩阵2（记录所有信息，包括临床数据）；2.对矩阵1中探针ID进行基因名注释，并去除没有注释上的探针ID，得到矩阵3（只有探针ID和注释基因名）；3.为了将矩阵2和3能对应上，则将矩阵2中没有注释上的探针ID一并删除（将临床信息与矩阵3匹配上）；4.
新建R project，命名为GEOOO7625（随意起名，建议为GSE76275会更加方便管理）。
常规修改镜像（下载速度会变快哦~）
rm(list = ls()) ## 魔幻操作，一键清空~
options()$repos
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
step 1 下载并载入GEO数据
## 1.获取GEO数据，将(F)改为(T)后即可运行，注意查看下载文件的大小，检查数据，如下图1
if(T){
library(GEOquery)
gset <- getGEO('GSE76275', destdir=".",
AnnotGPL = T, ## 注释文件
getGPL = T) ## 平台文件
save(gset,file='GSE76275_eSet.Rdata') ## 保存到本地
}
load('GSE76275_eSet.Rdata') ## 载入数据
图1
当然，如果你还不知道什么是GEO数据库，更加不知道如何查找数据并下载GEO数据，那么你可以点击这个教程（解读GEO数据存放规律及下载，一文就够 ）
简单来说，GEO就是NCBI旗下一个收录很多很多gene expression的公共数据库，你可以在上面下载别人上传好的各种数据，通过下载公共数据，自行挖掘数据。
step 2 提取所需的表达矩阵
gset是一个超级大的列表，里面拥有表达矩阵、分组信息、临床资料等等内容，因此我们需要从gset中提取中基因测序的表达信息，并附上分组已经注释信息，用于后续统计。
######### 导入数据后查看数据情况，行、列名，数据维度、目标数据的大致情况，比如我们想看不同分组的各个基因的表达量，那么你就要查看分组在哪一行/列，怎么分组的，基因是什么情况，只要在掌握已有数据架构的情况下，才能随心所欲的调整表达矩阵（以下代码可以自行增减，达到查看数据情况的目的即可）
b = gset[[1]] ## 降级提取b
exprSet=exprs(b) ## 获取表达矩阵（如下图2），发现已是log处理后的数据。通常表达矩阵的原始数字从0但好几百万都不等，需要进行归一化处理。14488875 elements
pdata = pData(b) ## 使用函数?pData获取样本临床信息（如性别、年龄、肿瘤分期等等） 265 obs. of 69
colnames(pdata) ## 查看列名，即查看所包含的临床信息类型，找到triple negative status在第67列
length(colnames(pdata)) ## 查看pdata列数 69
pdata[,67] ## 查看第67列triple negative status的数据情况,发现按照"TN"，"not TN"排列
group_list = as.character(pdata[, 67]) ## 将67列改成字符
dim(exprSet) ## 查看矩阵的维度 54675行 265列
table(group_list) ## 查看67列信息
not TN TN
67 198
图2
注意，这是已经归一化后的数据。
############# 调出所需的表达矩阵(正相思路)；移除不需要的数据，只留下想要的(反向操作思路)。
## 以下为正相思路
{
not_TN_expr = exprSet[, grep("not TN", group_list)]# 提出"not TN"所有行数据 3663225 elements
TN_expr = exprSet[, !(colnames(exprSet) %in% colnames(not_TN_expr))] # 提出TN数据 10825650
exprSet=cbind(not_TN_expr, TN_expr) ## 14488875 elements
## not_TN_expr：矩阵，行为探针ID，列为non_TN样本
## TN_expr：矩阵，行为探针ID，筛选出列为TN的样本
## exprSet变成了not_TN....TN....排序的矩阵，那么后面的group_list也应该相应排序，才能对上信息
}
############### 样本分组
group_list = c(rep('not_TN',ncol(not_TN_expr)), ## 1:265
rep('TN', ncol(TN_expr)))
## rep('not_TN',ncol(not_TN_expr))：连续打和矩阵列数not_TN_expr一样数量的non_TN
## rep('TN', ncol(TN_expr))：连续打和矩阵列数TN_expr一样数量的TN
save( exprSet, group_list, file = 'exprSet_by_group.Rdata')
######## 保存按照group_list分组的数据，我们得到了一个ID按照not TN, TN排序的大矩阵 exprSet！
其实这个数据本身就只有并且已经按照"TN"，"not TN"排列，所以对这个数据来说，不调整都是OK的，但还放上来这个步骤，是为了提醒大家对数据要仔细甄别；另外如果数据是只有"TN"，"not-TN"排列，这时候要把not—TN改成not_TN，因为"-"会让R语言模糊，不知道是不是一个减号还是一个符号，所以提醒大家在起名字的时候，千万不要用拥有功能的符号当做名字去用了，R语言会懵的。
step 3 对平台文件进行注释（gene ID注释上基因名）
## step 3. 对平台文件进行注释处理
################ 筛选探针
GPL = b@featureData@data # 先打开b=gset[[1]]的b查看一下，用@提取注释文件 54675 obs. of 21
colnames(GPL) # 看一下所有的列名，发现Gene symbol位于第三列，如下图3
ids = GPL[,c(1,3)] ## 只要探针ID(第一列)和基因名(第三列)，54675 obs of 2 variables(注意数据量变化)
ids = ids[ids[,2] != '',] ## 去NA注释(去除没有注释上的探针ID) 45118 obs of 2 variables
图3
############### 一个探针对应多个基因，请将它们分开~
library("plyr")
a<-strsplit(as.character(ids[,2]), "///") # 以///为分隔符，分开gene symbol，54675 elements，如图4
tmp <- mapply(cbind, ids[,1],a) # 探针ID（id第一列）和基因名（a)合并在一起 54675 elements
df <- ldply (tmp, data.frame) # 由于tmp是list，所以使用ldply循环将list变为数据框 58465 obs
ID2gene = df[,2:3] ## 提取探针ID和基因匹配的矩阵（如下图5），58456 obs of 2 variables
colnames(ID2gene) = c("id", "gene") # 修改列名为id,gene
dim(ID2gene) ## 58456 obs of 2 variables
save(ID2gene, file = 'ID2gene.Rdata') # 保存数据
######此时我们得到了探针ID对应上基因名的小矩阵ID2gene，即已经删除未被注释的ID
图4 图5
step 4 去除最初的表达矩阵exprSet中没有基因注释的探针
这时候我们需要好好看一眼ID2gene数据，会发现以下情况，一个探针对应多个基因名（这时候随机选择一个），还需要筛选同一个基因的不同表达量（通常我们只留下最大值）
############# 去除未注释上基因的探针ID
exprSet = exprSet[rownames(exprSet) %in% ID2gene[,1],] ## 11956270 elements
ID2gene = ID2gene[match(rownames(exprSet), ID2gene[,1]),] ## 45118 obs. of 2 variables
dim(exprSet) ## 检查匹配结果，如下图6
dim(ID2gene) ## 检查匹配结果，如下图6
tail(sort(table(ID2gene[,2])), n = 12L) ##
图6
######### 相同基因的表达数据取最大值
{
MAX = by(exprSet, ID2gene[,2],
function(x) rownames(x)[ which.max(rowMeans(x))])
MAX = as.character(MAX) ## 22836 elements
exprSet = exprSet[rownames(exprSet) %in% MAX,] ## 筛出已取最大值的表达矩阵 5765075 elements
rownames( exprSet ) = ID2gene[ match( rownames( exprSet ), ID2gene[ , 1 ] ), 2 ]
}
dim(exprSet)
### [1] 21755 265 原来为 45118 265
exprSet[1:5,1:5] ## 检查一下 下图7
save(exprSet, group_list, file = 'finalSet.Rdata') ## 保存宝贝数据！至此已经完成了数据调整。
图7
step 5 画PCA图
## 上代码
install.packages('ggfortify') ## 安装所需包，在这里ggfortify只是一个字符，所以加""
library(ggfortify) ## 加载包，在这里ggfortify作为一个包，是有功能性的，可以不加""
data = as.data.frame(t(exprSet)) ## 265 obs. of 21755 variables
data$group = group_list ## 265 obs. of 21756 variables
png('pca_plot.png', res=100) ## res=100就是图片的分辨率
autoplot(prcomp(data[, 1:(ncol(data)-1)]), data = data, colour = 'group') + theme_bw() ## 加灰色背景
dev.off() ## 保存png图，查看图8，完美重现PCA图
出了图，除了傻乐之外，还需要想一想，什么样的PCA算是合格的（思考脸），没想清楚，找了大大大大师兄帮忙看，说是大体有分开趋势，并且PCA相对来说不太需要验证，emm…
step 6--基因表达差异分析DEG
接上部分，我们得到了根据not TN, TN排序好、已将探针ID注释上基因名的矩阵，并且已经删除了未注释上的探针ID，你问我为啥要删？好咩，这个探针ID还没有名字，就更加没有信号通路和KEGG之类的注释了，也就是说，它现在还是一个神秘的基因且我们无法解读其功能---放在下游分析里面就是没有用！
BiocInstaller::biocLite('limma') ## 使用bioconductor安装包
library(limma) ## 载入包
design <- model.matrix(~0 + factor(group_list)) ## 提取矩阵，如下图9
colnames(design) = levels(factor(group_list))
rownames(design) = colnames(exprSet) ## 替换行名为样本ID
design ## 检查，如图10，只是得到了一个矩阵，还需要把表达量匹配上
contrast.matrix <- makeContrasts("TN-not_TN", levels = design)
> contrast.matrix ## 往下看，not_T用-1来表示，TN用1来表示
Contrasts
Levels TN-not_TN
not_TN -1
TN 1
load("~/Rtrain/GEOOO/ID2gene.Rdata") # 载入数据
####### 做DEG差异分析，下面的代码我也看不懂了，就跑吧...
{
fit <- lmFit(exprSet, design) ## lmFit看名字就是线性拟合（不就是一种统计方式么）
fit2 <- contrasts.fit(fit, contrast.matrix) ## 计算标准误差
fit2 <- eBayes( fit2 ) ## Empirical Bayes Statistics for Differential Expression
nrDEG = topTable( fit2, coef = 1, n = Inf )
write.table( nrDEG, file = "nrDEG.out")
}
head(nrDEG)
图9
图10
step 7 差异分析--heatmap
用热图来做大数据可视化最好用啦~我们还有初阶小白学热图的教程，具体可见：一招学会热图
step 8 差异分析-火山图
后续代码阅读原文即可！
■ ■ ■
生信基础知识大全系列：生信基础知识100讲
史上最强的生信自学环境准备课来啦！！ 7次改版，11节课程，14K的讲稿，30个夜晚打磨，100页PPT的课程。
如果需要组装自己的服务器；代办生物信息学服务器
如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等)，点我？
如果需要线下辅导及培训，看招学徒
如果需要个人电脑：个人计算机推荐
如果需要置办生物信息学书籍，看：生信人必备书单
如果需要实习岗位：实习职位发布
如果需要售后：点我
如果需要入门资料大全：点我