GEO数据挖掘，表达芯片分析
举例：王同学近期拟通过生物信息学相关软件与数据库来探讨女性非抽烟者的非小细胞肺癌预后相关的显著性基因及潜在的治疗靶点，他在NCBI上查询到了1套芯片数据GSE19804。请帮助他完成该项目的设计与分析。
一、一般流程
1、找数据，找到GSE编号
2、下载数据：包括表达矩阵、临床信息、分组信息
3、数据探索：分组之间是否有差异，PCA，热图
4、limma差异分析及可视化：P值、logFC、火山图、热图
5、富集分析KEGG、GO
二、安装R包
这是表达芯片分析需要用到的R包，其中有些包需要根据分析对象的不同而替换
options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
cran_packages <- c('tidyr',
'tibble',
'dplyr',
'stringr',
'ggplot2',
'ggpubr',
'factoextra',
'FactoMineR',
'devtools',
'patchwork')
Biocductor_packages <- c('GEOquery',
'hgu133plus2.db',#根据需要替换
'ggnewscale',
"KEGG.db",
"limma",
"impute",
"GSEABase",
"GSVA",
"clusterProfiler",
"org.Hs.eg.db",#根据需要替换
"preprocessCore",
"enrichplot",
"ggplotify")
for (pkg in cran_packages){
if (! require(pkg,character.only=T) ) {
install.packages(pkg,ask = F,update = F)
require(pkg,character.only=T)
}
}
for (pkg in Biocductor_packages){
if (! require(pkg,character.only=T) ) {
BiocManager::install(pkg,ask = F,update = F)
require(pkg,character.only=T)
}
}
#前面的所有提示和报错都先不要管。主要看这里
for (pkg in c(Biocductor_packages,cran_packages)){
require(pkg,character.only=T)
}
#没有error就是成功！
install.packages('tweenr')
library(tweenr)
#哪个报错，就回去安装哪个。如果你没有安装xx包，却提示你xx包不存在，这也正常，是因为复杂的依赖关系，缺啥补啥。
if(!require(AnnoProbe))devtools::install_local("./AnnoProbe-master.zip",upgrade = F)
library(AnnoProbe)
三、下载数据
rm(list=ls())#清空页面
options(stringsAsFactors=F)
f='Gse19804_eSet.Rdata'
if (!requireNamespace("BiocManager",quietly = TRUE))
install.packages("BiocManager")#下载BiocManager包
if (!requireNamespace("GEOquery",quietly=TRUE))
BiocManager::install("GEOquery")#下载GEOquery包
library(GEOquery)
if(!file.exists(f)){
gset<-getGEO('GSE19804',destdir=".",
AnnotGPL = F,#注释文件
getGPL = F) #平台文件
save(gset,file=f)} #保存到本地
load("GSE19804_eSet.Rdata")
class(gset)
length(gset)
gset=gset[[1]]
exp<-exprs(gset)
View(gset)
exp[1:4,1:4]#查看数据，判断是否需要取log2
exp=log2(exp+1)
四、提取信息
提取临床信息：
pd<-pData(gset)
调整pd的行名顺序与exp列名完全一致：
p=identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
提取芯片平台编号：
gpl_number<-gset@annotation
save(f,pd,gpl_number,file="steploutput.Rdata")
分组，这里按照对照和患者分组，对照在前，患者在后：
library(stringr)
Group<-ifelse(pd$characteristics_ch1 == "tissue: paired normal adjacent","non_cancer", "cancer")
save(exp,Group,file="Expreset1.Rdata")
五、数据检验
1、箱线图
rm(list=ls())
options(stringsAsFactors = F)
load(file = "Expreset1.Rdata")
View(exp)
boxplot(exp[,1:100],cex=0.2,col="red")
箱线图查看数据是否有异常值，可以看出成分集中，位于4-8之间，总位数在6左右，此时不需要进行标准化与归一化
之前做的没有保存，这张是GEO2R上生成得，大体类似
2、PCA主成分分析
首先，需要对数据进行降维，54,675进行降维分类
rm(list=ls())
options(stringsAsFactors = F)
load(file = "Expreset1.Rdata")
exp[1:4,1:4]
exp<-t(exp)#矩阵转置，画PCA图要求行名是样本名，列名是探针名
exp[1:4,1:4]
exp<-as.data.frame(exp)#转换格式为data.frame
exp<-cbind(exp,Group)#cbind横向追加，即将分组信息追加到最后一列
需要安装两个包：
install.packages("FactoMineR")
library("FactoMineR")
install.packages("factoextra")
library("factoextra")
进行降维：
exp.pca<-PCA(exp[,-ncol(exp)],graph = FALSE)#进行降维
head(exp.pca$ind$coord)#查看降维结果
画图：
fviz_pca_ind(exp.pca,
geom.ind = "point",#只显示点
col.ind = exp$Group,#按分组划定不同颜色
palette = c("#00AFBB","#E7B800"),
addEllipses = TRUE,
legend.title = "Groups"
)
输出的结果显示，集中程度一般，处于平面维度差异不大，分布不是特别开，两者有重叠的部分。
备注：这里分组不是很明显，我后面有再试过其他的方法，发现还是这样，如果有人可以指导一下非常感激~
3、层次聚类分析
rm(list=ls())
options(stringsAsFactors = F)
load(file = "Expreset1.Rdata")
exp<-t(exp)
d<-dist(exp)#计算距离
fit.complete<-hclust(d,method='complete')
plot(fit.complete,hang= -1,cex=0.8)
可以看出，基本可以聚成两类，也就是肿瘤组和非肿瘤组。
六、探针ID转换为基因ID
1、探针转换为基因ID并合并
rm(list=ls())
options(stringsAsFactors = F)
load(file = "Expreset1.Rdata")
if (!requireNamespace("BiocManager",quietly = TRUE))
install.packages("BiocManager")
if (!requireNamespace("hgu133plus2.db",quietly=TRUE))
BiocManager::install("hgu133plus2.db")
library("hgu133plus2.db")
ids=toTable(hgu133plus2SYMBOL)#提取GPL芯片平台的探针ID和对应的基因
length(unique(ids$symbol))#查看一共多少个基因
tail(sort(table(ids$symbol)))#查看每个基因对应多少个探针
经过查询，对应的平台的GPL570，该平台对应的R包是hgu133plus2
able() 函数可以生成频数统计表，这里就是统计每个基因symbol出现的次数然后将其表格化；sort()函数将symbol出现的频率从小到大进行排序；tail()取最后6个即出现频率最大的6个
table一下我们可以看到，9786基因设计了1个探针；5326个基因设计了2个探针；2879个基因设计了3个探针…也就是说大部分的基因只设计了1个探针。
其实一般基因都会设计很多探针的，我们下载的表达矩阵是作者处理之后的，把许多不好的探针都过滤掉了，我们处理作者的数据要默认人家做的是对的，否则就要下载原始数据自己处理。
table(sort(table(ids$symbol)))
table(rownames(exp) %in% ids$probe_id)
发现有11574个探针没有对应的基因名；43101个探针有对应的基因名。
x %in% y表示 x 的元素在y中吗？然后返回逻辑值。rownames(exp)即表达矩阵exp的行名是文章数据中用到的所有探针ID（probe_id）；ids$probe_id是具有对应基因的所有探针。所以返回的TRUE就是文章数据中有对应基因的探针数。
现在我们对探针进行过滤，把没有对应基因名的探针过滤掉：
exp = exp[rownames(exp) %in% ids$probe_id,]
dim(exp)
过滤的本质就是矩阵取子集，如：matrix[2,]意思就是取矩阵matrix的第2行和所有的列。同理，我们这里exp[rownames(exp) %in% ids$probe_id,]就是取具有对应基因的所有探针（行），和所有的列。
对比一下过滤之前和过滤之后的探针数量：
table(rownames(exp) %in% ids$probe_id)
dim(exp)
exp = exp[rownames(exp) %in% ids$probe_id,]
dim(exp)
可以发现过滤前54675个探针，过滤后剩下43101个探针。
然后，我们使用match函数把ids里的探针顺序改一下，使ids里探针顺序和我们表达矩阵的顺序完全一样。
ids=ids[match(rownames(exp),ids$probe_id),]
然后使用merge函数将exp的探针与芯片平台探针id相匹配，合并到exp
exp<-as.data.frame(exp)
exp$probe_id<-rownames(exp)
exp<-merge(exp,ids,by="probe_id")
2、多个探针检测一个基因，合并在一起，取其平均值
library(limma)
exp<-avereps(exp,ID=exp$symbol)
exp<-as.data.frame(exp)
row.names(exp)<-exp$symbol#将exp行名改为基因名
save(exp,Group,file="Expreset1.Rdata")
七、基因的差异性分析
rm(list=ls())
options(stringsAsFactors = F)
load(file = "Expreset1.Rdata")
library(limma)
exp<-avereps(exp,ID=exp$symbol)
exp<-as.data.frame(exp)
exp<-exp[,-c(1,122)]#去掉第一列与最后一列
save(exp,Group,file="Expreset1.Rdata")
exp[1:4,1:4]#最好每次都检测数据
table(Group)#table函数，查看Group中的分组个数
exp<-data.frame(exp,stringsAsFactors=F)#转为数量型
exp<-data.matrix(exp)
boxplot(exp[2,] ~Group)#第一个基因以Group分组画箱线图
t.test(exp[2,] ~Group)#第一个基因做t检验
#也可以自定义一个函数画图
bp=function(g){
library(ggpubr)
df=data.frame(gene=g,stage=Group)
p<-ggboxplot(df,x="stage",y="gene",
color="stage",palette="jco",
add="jitter")
p+stat_compare_means()#增加P值
}
#调用定义好的函数，避免同样的绘图代码重复多次
library(ggplot2)
bp(exp[1,])
bp(exp[2,])
bp(exp[3,])
bp(exp[4,])
以第一个基因为例，结果：
差异分析走标准的limma流程：
#创建一个分组的矩阵
library(limma)
design=model.matrix(~0+factor(Group))
colnames(design)=levels(factor(Group))
rownames(design)=colnames(exp)
head(design)
#创建差异比较矩阵
contrast.matrix<-makeContrasts(paste0(unique(Group),collapse="-"),levels=design)
contrast.matrix#这个矩阵声明，我们要肿瘤组和非肿瘤组进行差异分析比较
#第1步lmFit,为每个基因给定一系列的阵列来拟合线性模型
fit<-lmFit(exp,design)
#第2步eBayes，给出了一个微阵列线性模型拟合，通过经验贝叶斯调整标准误差到一个共同的值来计算修正后的t统计量、修正后的f统计量和微分表达式的对数概率
fit2<-contrasts.fit(fit,contrast.matrix)
fit2<-eBayes(fit2)
#第3步topTabe，从线性模型拟合中提取出排名靠前的基因表
#topTable(fit2,coef=2,adjust='BH')
tempOutput<-topTable(fit2,coef=1,n=Inf)
tempOutput<-na.omit(tempOutput)#移除NA值
head(tempOutput)
tempOutput["YY1",]
#上调基因和下调基因
tempOutput$g=ifelse(tempOutput$P.Value>0.01,"stable",
ifelse(tempOutput$logFC>1.5,"up",
ifelse(tempOutput$logFC<1.5,"down","stable")))
table(tempOutput$g)
save(exp,Group,tempOutput,file="tempOutput.Rdata")
说明：if判断，如果这一基因P值大于0.01，则为stable接上句else否则：接下来开始判断哪些P小于0.01的基因，再if判断，如果logFC大于1.5，则为上调基因
降维结果为：
八、绘制火山图和热图
火山图
rm(list=ls())
options(stringsAsFactors = F)
load(file = "tempOutput.Rdata")
DEG=tempOutput
head(DEG)
attach(DEG)
plot(logFC,-log10(P.Value))#df新增加一列“v”，值为-log10(P.Vaule)
library(ggpubr)
df=DEG
df$v=-log10(P.Value)
ggscatter(df,x="logFC",y="v",size=0.5)
ggscatter(df,x="logFC",y="v",size=0.5,color="g")#根据基因上调下调添加颜色
均值差异图
ggscatter(df,x="AveExpr",y="logFC",size=0.2)
df$p_c=ifelse(df$P.Value<0.001,"p<0.001",
ifelse(df$P.Value<0.01,"0.001<p<0.01","p<0.01"))
table(df$p_c)
ggscatter(df,x="AveExpr",y="logFC",color="p_c",size=0.2,
palete=c("green","red","black"))
热图
rm(list=ls())
options(stringsAsFactors = F)
load(file = "tempOutput.Rdata")
exp[1:4,1:4]
table(Group)#挑选出差异最大100个基因和差异最小的100个基因
x=tempOutput$logFC#去这一列将其重新赋值给x
names(x)=row.names(tempOutput)#将基因名作为名字，命名给x
x[1:4]
cg=c(names(head(sort(x),100)),
names(tail(sort(x),100)))#x从小到大排列，取前100及后100，并取其对应的基因名，作为向量赋值给cg
library(pheatmap)
pheatmap(exp[cg,],show_colnames=F,show_rownames=F)#按照cg取行，所得到的矩阵来画热土
n=t(scale(t(exp[cg,])))#通过scale对log-ratio数值进行归一化，scale()函数需要行使样本名，列是基因名，所以需要转置，然后再转置回来
n[n>2]=2
n[n<-2]=-2
n[1:4,1:4]
pheatmap(n,show_colnames=F,show_rownames=F)
#以肿瘤与非肿瘤分组画热图
ac=data.frame(g=Group)
rownames(ac)=colnames(n)#将ac的行名也就是分组信息给到n的列明，即热图中位于上方的分组信息
pheatmap(n,show_colnames=F,
show_rownames=F,
cluster_cols=F,
annotation_col=ac)#列名注释信息为ac即分组信息
