芯片数据分析流程
2021-01-28 本文已影响0人
白米饭睡不醒
1
0.安装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就是成功！
#哪个报错，就回去安装哪个。如果你没有安装xx包，却提示你xx包不存在，这也正常，是因为复杂的依赖关系，缺啥补啥。
if(!require(AnnoProbe))devtools::install_local("./AnnoProbe-master.zip",upgrade = F)
library(AnnoProbe)
#上次遇到的报错解决：
BiocManager::install("graphlayouts")
library(clusterProfiler)
library(enrichplot)
1.下载数据，提取表达矩阵和临床信息
#数据下载
rm(list = ls())#清空环境
options(stringsAsFactors = F)#避免因子影响
library(GEOquery)#加载要用的包
gse_number = "GSE56649"#给gse编号一个变量，后面只需修改这里就好
eSet <- getGEO(gse_number,
destdir = '.',
getGPL = F)#下载且读取，destdir从哪读取文件"."代表当前目录;getGPL意思是下载matrix文件时要不要一起下载GPL注释
class(eSet)#看是什么数据类型
length(eSet)#看长度
eSet = eSet[[1]]#简化并重新赋值为简化版
#eSet@phenoData#提取临床信息
eSet@annotation#返回这个表达矩阵用哪个平台测序的
#(1)提取表达矩阵exp
exp <- exprs(eSet)
exp[1:4,1:4]
#range(exp)#可以查看矩阵取值范围，看有没有log过
exp = log2(exp+1)#确认好有没有取过log，没取过再取，避免矩阵里有零加的1（不改）甲基化要加更小的值
boxplot(exp)#这里也用来看范围
#(2)提取临床信息
pd <- pData(eSet)
#(3)调整pd的行名顺序与exp列名完全一致
p = identical(rownames(pd),colnames(exp));p#判断pd的行名和exp是否一致
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平台编号
gpl_number <- eSet@annotation
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")#输出结果
-
检查数据下载的完整性
看两个KB或MB大小是否一样
2.分组信息&探针注释
# Group(实验分组)和ids(探针注释)
rm(list = ls())
load(file = "step1output.Rdata")
library(stringr)
# 1.Group----
# 第一类，有现成的可以用来分组的列
if(F) Group = pd$`disease state:ch1`
#第二类，自己生成
if(F){
Group=c(rep("RA",times=13),
rep("control",times=9))
}#（1）
rep(c("RA","control"),times = c(13,9))#与(1)相同，简单写法
#第三类，**匹配关键词，自行分类
Group=ifelse(str_detect(pd$source_name_ch1,
"control"),"control","RA")
#设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,
levels = c("control","RA"))
Group
# 注意levels与因子内容必须对应一致
# Group = pd$`disease state:ch1`
# Group = factor(Group,
# levels = c("healthy control","rheumatoid arthritis"))
#2.ids-----------------
#方法1 BioconductorR包(最常用)
gpl_number #用这个GPL编号在以下网址搜索，第三列包名字的前缀，后加.db
#http://www.bio-info-trainee.com/1399.html
#安装并加载这个包
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")#看包里有什么
ids <- toTable(hgu133plus2SYMBOL)#toTable是变成数据框
#SYMBOL代表探针和symbol之间的关系
head(ids)
# 方法2 读取GPL平台的soft文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
#注：soft文件列名不统一，活学活用，有的GPL平台没有提供注释，如GPL16956
a = getGEO(gpl_number,destdir = ".")
b = a@dataTable@table
colnames(b)
ids2 = b[,c("ID","Gene Symbol")]
colnames(ids2) = c("probe_id","symbol")
ids2 = ids2[ids2$symbol!="" & #空字符串的去掉
!str_detect(ids2$symbol,"///"),]#一个探针对应多个基因属于非特异性探针，要去掉
}
# 方法3 官网下载,文件读取
##http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus
# 方法4 自主注释
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,gse_number,file = "step2output.Rdata")
- 当一列中一部分相同，一部分不相同，取不同部分时
pd$characteristics_ch1
#两种方法
str_remove(pd$characteristics_ch1,"disease state:")
str_split(pd$characteristics_ch1,": ",simplify = T)[,2]
3.数据探索(pca)和热图
rm(list = ls())
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")
#输入数据：exp和Group
#Principal Component Analysis
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials
# 1.PCA 图----展示两组之间的差异
dat=as.data.frame(t(exp))#把数据转置并变成数据框
library(FactoMineR)
library(factoextra)
dat.pca <- PCA(dat, graph = FALSE)#逗号后不用管
pca_plot <- fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = Group, # color by groups
palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
pca_plot
ggsave(plot = pca_plot,
filename = paste0(gse_number,"_PCA.png"))#保存图片
save(pca_plot,file = "pca_plot.Rdata")#保存Rdata，为了后面合图
# 2.top 1000 sd 热图---- 从总基因中挑变化比较大的一部分基因
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,]
# 直接画热图，对比不鲜明
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n)
pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col
)
# 用标准化的数据画热图，两种方法的比较：https://mp.weixin.qq.com/s/jW59ujbmsKcZ2_CM5qRuAg
## 1.使用热图参数
pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
scale = "row",#按行标准化，只比较列与列的区别
breaks = seq(-3,3,length.out = 100)
) #breaks 参数解读在上面链接，设置颜色分配范围
#breaks：小于-3的值显示和-3一个颜色，大于3显示和3一个颜色
# 常用-3到3 或 -2到2
dev.off()
## 2.自行标准化再画热图，与上面1画出的图没有区别
n2 = t(scale(t(n)))#scale只能按列标准化，所以先转置
pheatmap(n2,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
breaks = seq(-3,3,length.out = 100)
)
dev.off()
# 关于scale的进一步探索：zz.scale.R
# 3.相关性热图----
pheatmap::pheatmap(cor(exp),#看表达矩阵组成新矩阵列与列之间的相关性
annotation_col = annotation_col)
#一般选第二个，或者做完差异分析后拿差异基因做相关性热图比较明显
pheatmap::pheatmap(cor(n),#top1000差异基因的相关性，比表达矩阵有意义一些
annotation_col = annotation_col)
pheatmap::pheatmap(cor(n2),#标准化后的top1000的相关性
annotation_col = annotation_col
)
dev.off()
# 关于相关性背后的故事：https://mp.weixin.qq.com/s/IqMW6Qjf64dn30F4RQg5kQ
2
2
数据框、矩阵转置后都会变成矩阵，as.data.frame把矩阵变成数据框
4.差异分析
rm(list = ls())
load(file = "step2output.Rdata")
#差异分析，用limma包来做
#需要表达矩阵和Group，不需要改
library(limma)
design=model.matrix(~Group)#group要求对照组在前，实验组在后
fit=lmFit(exp,design)#线性拟合，数据：表达矩阵和根据分组信息生成的比较矩阵design
fit=eBayes(fit)#贝叶斯检验
deg=topTable(fit,coef=2,number = Inf)
#为deg数据框添加几列
#1.加probe_id列（探针id），把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)
#2.加上探针注释，每个探针对应的基因
table(!duplicated(ids$probe_id))
table(!duplicated(ids$symbol))
#按symbol列去重，常见标准有3个：最大值/平均值/随机去重
#随机去重，另两个见zz.去重方式.R
ids = ids[!duplicated(ids$symbol),]
deg <- inner_join(deg,ids,by="probe_id")#对deg和ids按照基因取交集
head(deg)
nrow(deg)#基因的数量
#3.加change列,标记上下调基因
logFC_t=1#设置logFC的阈值
P.Value_t = 0.01#设置P.Value的阈值
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)#下调基因
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)#上调基因
change = ifelse(k1,
"down",
ifelse(k2,"up","stable"))#标记上下调基因
deg <- mutate(deg,change)#添加到deg表格中
#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人类
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
dim(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
dim(deg)
length(unique(deg$symbol))
save(Group,deg,logFC_t,P.Value_t,gse_number,file = "step4output.Rdata")
5.可视化：火山图和热图
rm(list = ls())
load(file = "step1output.Rdata")
load(file = "step4output.Rdata")
#1.火山图----
library(dplyr)
library(ggplot2)
dat = deg#改名了
p <- ggplot(data = dat,
aes(x = logFC,
y = -log10(P.Value))) +
geom_point(alpha=0.4, size=3.5,
aes(color=change)) +
ylab("-log10(Pvalue)")+
scale_color_manual(values=c("blue", "grey","red"))+#改颜色
geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +#竖线
geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +#横线
theme_bw()
p
#在图上标记关心的基因
if(T){
#自选基因
for_label <- dat%>%
filter(symbol %in% c("HADHA","LRRFIP1"))
}
if(F){
#p值最小的10个
for_label <- dat %>% head(10)
}
if(F) {
#p值最小的前3下调和前3上调
x1 = dat %>%
filter(change == "up") %>%
head(3)
x2 = dat %>%
filter(change == "down") %>%
head(3)
for_label = rbind(x1,x2)
}
volcano_plot <- p +
geom_point(size = 3, shape = 1, data = for_label) +
ggrepel::geom_label_repel(#添加标记基因的图层
aes(label = symbol),
data = for_label,
color="black"
)
volcano_plot
ggsave(plot = volcano_plot,filename = paste0(gse_number,"_volcano.png"))
#2.差异基因热图----
load(file = 'step2output.Rdata')
if(T){
#全部差异基因
cg = deg$probe_id[deg$change !="stable"]
length(cg)
}else{
#取前30上调和前30下调
x=deg$logFC[deg$change !="stable"]
names(x)=deg$probe_id[deg$change !="stable"]
cg=names(c(head(sort(x),30),tail(sort(x),30)))
length(cg)
}
n=exp[cg,]
dim(n)
#差异基因热图
library(pheatmap)
annotation_col=data.frame(group=Group)#变成输入数据要求的样子
rownames(annotation_col)=colnames(n)
heatmap_plot <- pheatmap(n,show_colnames =F,
show_rownames = F,#要显示行名就注释掉，再将探针名替换为基因名
scale = "row",
#cluster_cols = F, #按列聚类，显示就不按列聚类也可
annotation_col=annotation_col,
breaks = seq(-3,3,length.out = 100)
)
heatmap_plot
ggsave(heatmap_plot,filename = paste0(gse_number,"_heatmap.png"))
load("pca_plot.Rdata")
library(patchwork)
library(ggplotify)
(pca_plot + volcano_plot +as.ggplot(heatmap_plot))
6.富集分析
rm(list = ls())
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)
# 1.GO 富集分析----
#(1)输入数据
gene_up = deg[deg$change == 'up',
'ENTREZID'] #上调基因对应的ENTREZID
gene_down = deg[deg$change == 'down','ENTREZID']
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']
#(2)富集
#以下步骤耗时很长，设置了存在即跳过
if(!file.exists(paste0(gse_number,"_GO.Rdata"))){
ego <- enrichGO(gene = gene_diff,#输入数据
OrgDb= org.Hs.eg.db,
ont = "ALL",
readable = TRUE)
#ont参数：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
save(ego,file = paste0(gse_number,"_GO.Rdata"))
}
load(paste0(gse_number,"_GO.Rdata"))
#(3)可视化
#条带图
barplot(ego)
#气泡图
dotplot(ego)
dotplot(ego, split = "ONTOLOGY", font.size = 10,
showCategory = 5) + facet_grid(ONTOLOGY ~ ., scale = "free") +
scale_y_discrete(labels = function(x) str_wrap(x, width = 45))
#geneList 用于设置下面图的颜色
geneList = deg$logFC
names(geneList)=deg$ENTREZID
geneList = sort(geneList,decreasing = T)
#(3)展示top通路的共同基因，要放大看。
#Gene-Concept Network
cnetplot(ego,categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(ego, showCategory = 3,foldChange=geneList, circular = TRUE, colorEdge = TRUE)
#Enrichment Map,这个函数最近更新过，版本不同代码会不同
#showCategory = 3，展示条目的数量，默认是五
Biobase::package.version("enrichplot")
if(F){
emapplot(pairwise_termsim(ego)) #新版本
}else{
emapplot(ego)#老版本
}
#(4)展示通路关系 https://zhuanlan.zhihu.com/p/99789859
#goplot(ego)
#(5)Heatmap-like functional classification
heatplot(ego,foldChange = geneList,showCategory = 8)
# 2.KEGG pathway analysis----
#上调、下调、差异、所有基因
#（1）输入数据
gene_up = deg[deg$change == 'up','ENTREZID']
gene_down = deg[deg$change == 'down','ENTREZID']
gene_diff = c(gene_up,gene_down)
gene_all = deg[,'ENTREZID']
#（2）对上调/下调/所有差异基因进行富集分析
if(!file.exists(paste0(gse_number,"_KEGG.Rdata"))){
kk.up <- enrichKEGG(gene = gene_up,
organism = 'hsa')#物种的缩写
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa')
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa')
save(kk.diff,kk.down,kk.up,file = paste0(gse_number,"_KEGG.Rdata"))
}
load(paste0(gse_number,"_KEGG.Rdata"))
#(3)看看富集到了吗？https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
table(kk.diff@result$p.adjust<0.05)
table(kk.up@result$p.adjust<0.05)
table(kk.down@result$p.adjust<0.05)
#(4)按照pvalue筛选通路
down_kegg <- kk.down@result %>%
filter(pvalue<0.05) %>% #筛选行
mutate(group=-1) #新增列
up_kegg <- kk.up@result %>%
filter(pvalue<0.05) %>%
mutate(group=1)
#(5)可视化
source("kegg_plot_function.R")
g_kegg <- kegg_plot(up_kegg,down_kegg)
g_kegg
#g_kegg +scale_y_continuous(labels = c(4,2,0,2,4))#改横坐标用的，横坐标不应该为负
ggsave(g_kegg,filename = 'kegg_up_down.png')
# 3.gsea作kegg和GO富集分析----
## https://www.jianshu.com/p/c5b7b7dbf29b
#(1)查看示例数据
data(geneList, package="DOSE")
#(2)将我们的数据转换成示例数据的格式
geneList=deg$logFC
names(geneList)=deg$ENTREZID
geneList=sort(geneList,decreasing = T)
#(3)富集分析
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
verbose = FALSE)
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
#(4)可视化
g2 = kegg_plot(up_kegg,down_kegg)
g2
# 4.能看懂的资料越来越多----
# GSEA学习更多：https://www.jianshu.com/p/baf85b51752e
# 富集分析学习更多：http://yulab-smu.top/clusterProfiler-book/index.html
# 弦图：https://www.jianshu.com/p/e4bb41865b7f
# GOplot：https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew
# 网上的资料和宝藏无穷无尽，学好R语言慢慢发掘~
* Annoprobe包三个函数的使用
#[1]
library(AnnoProbe)
geoChina("gse1009")#下载gse1009文件
load("GSE1009_eSet.Rdata")
eSet=gset[[1]]#即可对接标准流程
#[2]
ids=idmap("GPL570")#获取ids的一种方法
head(ids)
#[3]
annoGene()#可以知道你提供的基因ID的具体信息是什么
*小洁老师自己写的R包（简化）
#devtools::install_github("xjsun1221/tinyarray")
library(tinyarray)
geo = geo_download("GSE56649")
View(geo$pd)
pd = geo$pd
library(stringr)
Group=ifelse(str_detect(pd$source_name_ch1,"control"),
"control",
"RA")
#设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,
levels = c("control","RA"))
Group
find_anno(geo$gpl,install = T)
ids <- toTable(hgu133plus2SYMBOL)
geo$exp = log2(geo$exp+1)
deg = get_deg_all(geo$exp,Group,ids)
deg$plots
7.常见错误
3
8.自行探索部分
(1）配对数据分析
[1]
[2]
(2）多分组数据
[3]
[4]
(3)多数据联合分析
[5]
[6]
[7]
[8]
(4)标准流程的后续
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
