•画图通用,仿制数据的思维通用,富集分析基本通用
•GEO数据库的背景知识
•GEO表达芯片的原理
•GEO表达芯片特有的下载方式
•表达芯片的差异分析(就几句代码)
•表达芯片的复杂分析
•主要学思维和方法,后面重点学习转录组的具体分析代码
·输入数据是一个连续型向量和一个有重复值的离散型向量—横坐标;
·上下五条线的意思 中间的又黑又粗的—中位数;上下两条线是最大值和最小值;方框的上下两条线是75%和25%(四分位数);在外面的点-离群点
>rm(list = ls())
>library(GEOquery)
#先去网页确定是否是表达芯片数据,不是的话不能用本流程。
>gse_number = "GSE56649"
>eSet <- getGEO(gse_number, destdir = '.', getGPL = F) #getGEO下载并读取文件
> class(eSet)
[1] "list"
> length(eSet)
[1] 1
> eSet = eSet[[1]]
> class(eSet)
[1] "ExpressionSet" #哪种数据类型
attr(,"package") # 属性之包的名字
[1] "Biobase" #Biobase包规定的格式"ExpressionSet"
> exp <- exprs(eSet)
> dim(exp)
[1] 54675 22
> exp[1:4,1:4] #检查矩阵是否正常,如果是空的就会报错,空的和有负值的、有异常值的矩阵需要处理原始数据。是否取过loglog 是否有负值
GSM1366348 GSM1366349 GSM1366350 GSM1366351
1007_s_at 279.156 202.866 406.818 232.668
1053_at 487.700 409.018 393.810 397.127
117_at 666.869 388.589 704.633 953.481
121_at 240.646 361.198 305.229 374.757
> #如果表达矩阵为空,大多数是转录组数据,不能用这个流程(后面另讲)。
> #自行判断是否需要log 可以画boxplot()
> exp = log2(exp+1) 正常log在0-20之间
> boxplot(exp)
### (2)提取临床信息
>pd <- pData(eSet)
#(3)让exp列名与pd的行名顺序完全一致
>p = identical(rownames(pd),colnames(exp));p#判断信息是否以一对应
>if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#分组信息来自临床信息,分组信息需要与表达矩阵列名一一对应
#临床信息需要与表达矩阵一一对应
>gpl_number <- eSet@annotation;gpl_number #提取子集 注意什么时候用@,什么时候用$,看图1
[1] "GPL570"
>save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")
#保存 gse_number(原本的编号),pd(临床信息),exp(表达矩阵),gpl_number(芯片编号)
# 从临床样本中获得实验分组(在表格中慢慢找,代码如何实现看下)
rm(list = ls())
load(file = "step1output.Rdata")
library(stringr)
# 标准流程代码是二分组,多分组数据的分析后面另讲
# 生成Group向量的三种常规方法,三选一,选谁就把第几个逻辑值写成T,另外两个为F。如果三种办法都不适用,可以继续往后写else if
# 1.Group----
# 第一种方法,有现成的可以用来分组的列
Group = pd$`disease state:ch1` #pd$`cell type:ch1`注意这个``不能删
}else if(F){
# 第二种方法,自己生成,但需要看清哪是RA哪是control
Group = c(rep("RA",times=13),
rep("control",times=9))
Group = rep(c("RA","control"),times = c(13,9)) #简写
}else if(T){
# 第三种方法,使用字符串处理的函数获取分组
Group=ifelse(str_detect(pd$source_name_ch1,"control"),
"control",
"RA")
#felse(str_detect()) 两个函数连用是用来分组的神器
#str_detect(pd$source_name_ch1,"control")-source_name_ch1这一列是否含有control
#ifelse(str_detect(pd$source_name_ch1,"control"),"control","RA")
是control就填control不是填RA
#三种代码完整方法
if(F){
# 1.Group----
# 第一种方法,有现成的可以用来分组的列
Group = pd$`disease state:ch1`
}else if(F){
# 第二种方法,自己生成
Group = c(rep("RA",times=13),
rep("control",times=9))
Group = rep(c("RA","control"),times = c(13,9))
}else if(T){
# 第三种方法,使用字符串处理的函数获取分组
Group=ifelse(str_detect(pd$source_name_ch1,"control"),
"control",
"RA")
}
# 需要把Group转换成因子,并设置参考水平,指定levels,对照组在前,处理组在后 ;因子正文与levels不对应时会产生NA
Group = factor(Group,levels = c("control","RA"))
Group
#Group是一个有重复值的向量 是分类型数据,适合用因子的形式
#factor直接转换并自动生成levels (control和RA),顺序以字母排序为准
#levels顺序有意义,在第一个位置的水平是参考水平
#参考水平将在做差异分析时,被设为对照组
#所以需要控制levels的顺序
#levels = c("control","RA") 写了按照写的顺序,control位参考水平
注释来源:不是所有的GPL都可以找到注释
>library(tinyarray) #tinyarry这个包是生信技能树写的,find_anno()是其中的一个函数
>find_anno(gpl_number) #打出找注释的代码
[1] "`library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)` (方法一,R包)and `ids <- AnnoProbe::idmap('GPL570')` are both avaliable"(方法二,idmap,idmap包含这个平台)
> ids <- AnnoProbe::idmap('GPL570')
trying URL 'http://49.235.27.111/GEOmirror/GPL/GPL570_bioc.rda'
Content type 'application/octet-stream' length 264292 bytes (258 KB)
==================================================
downloaded 258 KB
file downloaded in /Users/zhuo/learn /pipeline
>gpl_number
#http://WWW.bio-info-trainee.com/1399.html 网站(看图)上找GPL对应的R包编号,
>if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db") #下载R包
>library(hgu133plus2.db)
>ls("package:hgu133plus2.db") #只是展示包里有什么,但只要SYMBOL这个
>ids <- toTable(hgu133plus2SYMBOL) #ids <- toTable()提取 不同的包进行替换(看图)
>head(ids) #看到所需要的结果
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
#注:表格读取参数、文件列名不统一(如Gene Symbol 变成了Gene_Symbol),活学活用,有的表格里没有symbol列,也有的GPL平台没有提供注释表格
b = read.delim("GPL570-55999.txt",
check.names = F,
comment.char = "#")
colnames(b) #看图colnames(b)
ids2 = b[,c("ID","Gene Symbol")] #提取出来想要的两列
colnames(ids2) = c("probe_id","symbol") #统一重新命名,后面就不用改代码了 看图ids2 ;symbol中有的名称有/// 前后有两个名称,代表是两个基因-非特异性探针,直接去除
k1 = ids2$symbol!="";table(k1) #检测 空的 看图K1
k2 = !str_detect(ids2$symbol,"///");table(k2) #检测非特异性探针
ids2 = ids2[ k1 & k2,] #取既不是空的也不是非特异性的 & 同时满足两个条件
# ids = ids2 ids2是为了教课时好区分 后面用的时候把前面的#去掉,省改后面的代码
}
#若用第二种方法把最前面的if(F)变成T
仿制实例数据
列—两个部分(前四列是用于求PCA的值-探针/基因;最后一列为分组信息)
行—样本名称
需要对原始数据进行转换(如图a)
#仿制的前四列
dat=as.data.frame(t(exp)) #t() 转置 as.data.frame()作为数据框
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
save(pca_plot,file = "pca_plot.Rdata") #保存pca的数据
cg=names(tail(sort(apply(exp,1,sd)),1000)) # 挑选出1000个方差最大的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)
# 按行标准化
pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
scale = "row", #进行标准化,为了凸显行之间差别 缩小列之间的差别
breaks = seq(-3,3,length.out = 100) #breaks() -3,3(不同的结果设置的色带分配值不一样)是设置色带分布范围 分配颜色色带分配100种颜色
)
----注意 在代码没有错,但也不画图的情况下 用dev.off()-----
library(limma) #limma 包
design=model.matrix(~Group) #构建模型矩阵
fit=lmFit(exp,design) #做线性拟合
fit=eBayes(fit) #贝叶斯检验
deg=topTable(fit,coef=2,number = Inf) #topTable 提取结果
#简化思维 把上面的函数写成一个自己新的函数如egp()
egp=function(exp,Group){
design=model.matrix(~Group)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
return(deg)
}
a = egp(exp,Group)
## 鉴定两个函数结果是否一致
identical(a,deg)
#这一步结束得到的是deg(六列数据,还需4列,看图差异分析后的数据整理)
#为deg数据框添加几列
#1.加probe_id列,把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
#2.加上探针注释
ids = ids[!duplicated(ids$symbol),] # 这个代码是随机去重的方式
ids =distinct(ids,symbol,.keep_all = T)#这个代码也是随机去重的方式
###出现多个探针对应一个基因的情况,所以需对基因进行去重
####方法1:随机去重
####方法2:保留行和/行平均值最大的探针
####方法3:取多个探针的平均值
#其他去重方式在 “zz.去重方式.R”这个文件里
deg <- inner_join(deg,ids,by="probe_id") #inner_join 取交集
nrow(deg)
#3.加change列,标记上下调基因
logFC_t=1 #在这里更改阈值
p_t = 0.05 #在这里更改阈值
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable"))) #mutate加上updown...
table(deg$change)
#4.加ENTREZID列,用于富集分析(symbol转entrezid,然后inner_join)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, #bitr() 数据类型转换
fromType = "SYMBOL", #从..ID转换为
toType = "ENTREZID", # ...ID
OrgDb = org.Hs.eg.db) #转换一句的数据库-人类
#转换完结果看图b
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL")) #两个表的样本名称不一样,用inner_join() 后面要写by=c(""="")
save(Group,deg,logFC_t,p_t,gse_number,file = "step4output.Rdata") #保存
#到此有了10列 用于后续做图
library(dplyr)
library(ggplot2)
dat = distinct(deg,symbol,.keep_all = T)
p <- ggplot(data = dat,
aes(x = logFC, #横坐标
y = -log10(P.Value))) + #纵坐标
geom_point(alpha=0.4, size=3.5, #画点图 alpha透明度 size点的大小
aes(color=change)) + #这步之前是做了普通的点图
scale_color_manual(values=c("blue", "grey","red"))+ #设定颜色
geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +. #geom_vline()添加参考线(竖线)
geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) + #geom_hline()添加参考线(横线)
theme_bw() #去掉灰色的背景
p
#在火山图上标记感兴趣的基因
for_label <- dat%>%
filter(symbol %in% c("HADHA","LRRFIP1")) #挑选出想要的基因所在的行
volcano_plot <- p + #叠加涂层
geom_point(size = 3, shape = 1, data = for_label) + #画了空心的圆
ggrepel::geom_label_repel( #添加标签 (基因的名+边框)
aes(label = symbol),
data = for_label, #数据是来自新建的(挑选的基因所在的行,for_label)
color="black"
)
volcano_plot #结果图看图c
load(file = 'step2output.Rdata')
# 表达矩阵行名替换
exp = exp[dat$probe_id,] #转变为以基因为行名的
rownames(exp) = dat$symbol
if(T){
#取前10上调和前10下调 (可按logFC取也可按P value取)
library(dplyr)
dat2 = dat %>%
filter(change!="stable") %>% #筛选基因change不是stable的-找差异的
arrange(logFC) #按照logFC从小到大排序
cg = c(head(dat2$symbol,10), #取前十后十
tail(dat2$symbol,10))
}else{
#全部差异基因
cg = dat$symbol[dat$change !="stable"]
length(cg)
}
n=exp[cg,] #统计就会出现 20个差异表达基因和22个样本
dim(n)
#差异基因热图
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n)
heatmap_plot <- pheatmap(n,show_colnames =F,
scale = "row",
#cluster_cols = F, #如果基因数太多,就可以将基因名删掉
annotation_col=annotation_col,
breaks = seq(-3,3,length.out = 100)
)
heatmap_plot
# 拼图
library(patchwork)
volcano_plot +heatmap_plot$gtable #拼火山图和热图
library(corrplot)
g = sample(deg$symbol[1:500],10) # 这里是随机取样,注意换成自己感兴趣的基因
g
M = cor(t(exp[g,])) #cor()用于计算基因的相关性,提供矩阵数据,计算列于列之间的相关性,看图
pheatmap(M)
# 配色R包
library(paletteer)
my_color = rev(paletteer_d("RColorBrewer::RdYlBu")) #paletteer_d("RColorBrewer::RdYlBu") 在控制台上展示颜色,原本是蓝色-正相关,红色-负相关,rev将它转换一下 看图c
my_color = colorRampPalette(my_color)(10) #只取10种颜色
corrplot(M, type="upper",
method="pie", #以饼图代表相关系数大小
order="hclust",
col=my_color, #以颜色代表正负相关
tl.col="black",
tl.srt=45)
library(cowplot)
cor_plot <- recordPlot() #将画板上的图记录下来,为了后续可以拼图
# 拼图
load("pca_plot.Rdata")
plot_grid(pca_plot,cor_plot,
volcano_plot,heatmap_plot$gtable)
dev.off()
#保存
pdf("deg.pdf")
plot_grid(pca_plot,cor_plot,
volcano_plot,heatmap_plot$gtable)
dev.off()
a = deg$symbol[1]
boxplot(exp[a,]~Group)
deg$logFC[1]
-----来自生信技能树----
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。