#step 1--------下载GSE4107的芯片数据---------------
library(Biobase);library(BiocGenerics);library(parallel)#GEOquery的依赖包
library(GEOquery)
gse="GSE4107"#将GSE4107赋值给gse
eSet <- getGEO(gse,destdir = ".", getGPL = F)#使用getGEO函数下载表达数据(此时并不是表达矩阵)
#(1)使用exprs()函数提取表达矩阵
class(eSet)
length(eSet)
eSet[[1]]
length(eSet[[1]])
exp <- exprs(eSet[[1]])
class(exp)
colnames(exp)
rownames(exp)
exp[1:4,1:4]
#使用log2将芯片表达数据归一化exp=log2(exp+1)
exp=log2(exp+1)
exp[1:4,1:4]
dim(exp)
#(2)使用pData()函数从第一步中下载的eSet提取临床信息
pd <- pData(eSet[[1]])
class(pd)
dim(pd)
pd[1:4,1:4]
colnames(pd)
rownames(pd)
colnames(exp)
#(3)调整pd的行名和exp的列名完全一致
identical(rownames(pd),colnames(exp))
p <- identical(rownames(pd),colnames(exp))#检验pd行名是否和exp列名一致
match(rownames(pd),colnames(exp))
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]#如不一致,需运行此行命令
#(4)提取芯片平台编号,平台编号列在annotation列
gpl <- eSet[[1]]@annotation
#保存提取的表达矩阵、临床数据(主要是分组信息)、芯片平台编号供下一步使用
save(gse,exp,pd,gpl,file = "step1out.Rdata")
#。。。。。。至此,第一步获取表达数据和分组信息已完成。。。。。。。。。。。。。。。
#。。。。。。。。。。。。。。。。。。。。。。。。。
#step 2------分组-------------------
#(1)group_list(实验分组)和ids(芯片注释),主要是字符串操作,每次都需要改
rm(list = ls())
load("step1out.Rdata")
library(stringr)
pd$title#使用此命令查看pd(临床信息)的分组信息
group_list=ifelse(str_detect(pd$title,"control"),"control","treat")
dim(pd)
pd[1:4,1:4]
#使用factor()函数设置参考水平,对照在前,处理在后,谁在前谁是参考
group_list <- factor(group_list,levels = c("control","treat"))
#(2)--------ids---------
#方法1 BioconductorR包
gpl
#http://www.bio-info-trainee.com/1399.html
#此时需要网页搜索GPL对应的包,找到后就可以找到芯片探针ID和基因ID
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
#用toTable函数将带有symbol的gene生成数据框
ids <- toTable(hgu133plus2SYMBOL)
class(ids)
head(ids)
#至此获得芯片探针ID和基因symbol
#保存exp/ids/group_list/至step2.Rdata,下一步就是PCA/heatmap差异分析了
save(exp,ids,group_list,file = "step2out.Rdata")
#--------------美妙的分割线---------------
#载入第一步和第二步的数据
rm(list = ls())
load("step1out.Rdata")
load("step2out.Rdata")
#PCA差异分析只需要exp(表达矩阵)和group_list(分组信息)
#Principal Component Analysis
#输入数据:exp和group_list
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials
#画图需要数据库,而exp是矩阵,需将其转化成数据库
#先看一下exp数据,确定横纵坐标,PCA横坐标是GSM(样本),纵坐标是group_list
class(exp)
colnames(exp)
rownames(exp)
dat=as.data.frame(t(exp))#先将exp表达矩阵转置,后将其转化成数据框
library(FactoMineR)#画主成分分析图需要加载这两个包
library(factoextra)
# pca的统一操作走起
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_list, # color by groups
#palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
pca_plot
#使用ggplot 中的ggsave()函数保存PCA
ggsave(plot = pca_plot,filename = paste0(gse,"PCA.png"))
#将图保存为R.data,以供后续拼图用
save(pca_plot,file = "pca_plot.Rdata")
#热图 ------
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,]
#绘制热图
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n)
library(pheatmap)
heatmap=pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
scale = "row")
ggsave(plot = heatmap,filename = paste0(gse,"heatmap.png"))
save(heatmap,file = "heatmap.Rdata")
save(heatmap,cg,n,file = "heatmap.Rdata")
#至此初步分析芯片分组情况结束,下一步是挑选表达差异基因
#-----------------华丽丽的分割线--------------
https://blog.sciencenet.cn/blog-3414136-1229004.html
上一篇:
学习R语言的第一步:安装软件+安装包下一篇:
我的第一个转录组分析