admin管理员组文章数量:1794759
文章复现—bulkRNA转录组结合机器学习等进行相关疾病研究01—多数据集去除批次效应后联合分析以及火山图标准绘制
专题1—文章复现—bulkRNA转录组结合机器学习等进行相关疾病研究01—多数据集去除批次效应后联合分析以及火山图标准绘制
今天给大家复现一篇bulkRNA转录组结合机器学习等进行相关疾病研究的文章
《Identification of cuproptosisassociated subtypes and signature genes for diagnosis and risk prediction of Ulcerative colitis based on machine learning》基于机器学习的溃疡性结肠炎诊断和风险预测中铜死亡相关亚型和特征基因的鉴定（IF：5.7） Date：2023.04
1 文章思路
文章的大体思路如下：
可以看出这是一篇纯生信的文章，前期还是常规的GEO数据挖掘，取了三个与UC相关的数据集，处理后合并在一起，去除批次后，进行常规的差异基因和富集分析，之后与铜死亡相关基因取交集，取交集基因进行机器学习建模和亚群分类等。这篇文档主要介绍数据收集和预处理，差异表达基因鉴定，以及给大家复现一下原文章的Fig 2.A，后续的复现等之后再给大家一一复现。
2 数据收集与预处理
文章挑选了三个GSE38713 , GSE87473 , GSE92415，基于芯片的数据集，联合起来分析，共298个实验组，55个对照组，数据集的芯片平台并不相同，我们要先单独处理每个数据集后，拿到相应的表达矩阵（行名基因名，列名样本名）和分组信息后，才能根据基因名取交集，cbind后再去除批次效应。
1 GSE38713
代码语言：r复制rm(list = ls())
options(stringsAsFactors = F)
library(AnnoProbe)
library(GEOquery)
library(ggplot2)
library(ggstatsplot)
library(patchwork)
library(reshape2)
library(stringr)
getOption('timeout')
options(timeout=10000)
#获取并且检查表达量矩阵
##～～～gse编号需修改～～～
gse_number <- 'GSE38713'
#dir.create(gse_number)
#setwd( gse_number )
#list.files()
if(T){gset <- geoChina(gse_number)}
gset[[1]]
a=gset[[1]]
dat=exprs(a) #提取表达矩阵
dim(dat)#看一下dat这个矩阵的维度
dat[1:4,1:4]
pd = pData(a)
head(pd)
##～～～查看数据是否需要log～～～
boxplot(dat[,1:4],las=2)
range(dat)
#dat=log2(dat+1)
boxplot(dat,las=2)
#去批次
library(limma)
dat=normalizeBetweenArrays(dat)
#根据生物学背景及研究目的人为分组
##通过查看说明书知道取对象a里的临床信息用pData
colnames(pd)
library(stringr)
##～～～分组信息编号需修改～～～
group_list=ifelse(grepl("control", pd$description,ignore.case = T ),
"control","case")
table(group_list)
group_list <- factor(group_list, levels = c("control","case"))
table(group_list)
#查看芯片编号
a
a@annotation
# #获取芯片注释信息
# ##～～～芯片编号需修改～～～
ids=idmap( a@annotation ,'soft')
#
head(ids)
colnames(ids)=c('probe_id','symbol')
dim(ids)
ids=ids[ids$symbol != '',]
ids=ids[ids$probe_id %in% rownames(dat),]
ids$probe_id=as.character(ids$probe_id)
dat[1:4,1:4]
dat=dat[ids$probe_id,]
if(T){
#探针基因ID对应以及去冗余
dat[1:4,1:4]
ids=ids[ids$probe_id %in% rownames(dat),]
dat[1:4,1:4]
dat=dat[ids$probe_id,]
probe_exp = dat
probe_info = ids
#ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids$median=apply(dat,1,median)
#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
ids=ids[!duplicated(ids$symbol),]
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4] #保留每个基因ID第一次出现的信息
}
save(gse_number,dat,group_list,
#probe_exp,probe_info,pd,
file = 'step1_output.Rdata')2 GSE87473
注意 GSE87473 , GSE92415的GPL平台是GPL13158，不能通过ids=idmap( a@annotation ,'soft')函数直接获取，应该去GEO网站上找到对应的文件下载导入到R中。
代码语言：r复制rm(list = ls())
options(stringsAsFactors = F)
library(AnnoProbe)
library(GEOquery)
library(ggplot2)
library(ggstatsplot)
library(patchwork)
library(reshape2)
library(stringr)
getOption('timeout')
options(timeout=10000)
#获取并且检查表达量矩阵
gse_number <- 'GSE87473'
#dir.create(gse_number)
#setwd( gse_number )
#list.files()
if(T){gset <- geoChina(gse_number)}
#gset = getGEO("GSE12366", destdir = '.', getGPL = F)
gset[[1]]
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd = pData(a)
head(pd)
##～～～查看数据是否需要log～～～
boxplot(dat[,1:4],las=2)
range(dat)
#dat=log2(dat+1)
boxplot(dat,las=2)
library(limma)
dat=normalizeBetweenArrays(dat)
a=gset[[1]]
#根据生物学背景及研究目的人为分组
##通过查看说明书知道取对象a里的临床信息用pData
pd=pData(a)
##挑选一些感兴趣的临床表型。
colnames(pd)
pd$characteristics_ch1.2
library(stringr)
##～～～分组信息编号需修改～～～
group_list=ifelse(grepl("Normal", pd$characteristics_ch1.2,ignore.case = T ),
"control","case")
table(group_list)
group_list <- factor(group_list, levels = c("control","case"))
table(group_list)
a
a@annotation
# #获取芯片注释信息
# ##～～～芯片编号需修改～～～
ids=idmap( a@annotation ,'soft')
head(ids)
#上面的函数用不了，自己下载导入GPL文件
library(data.table)
tmp <- fread("GPL13158.txt", sep = "\t", skip = 18, data.table = FALSE)
tmp <- tmp[-nrow(tmp), ]
dim(tmp)
ids <- tmp[,c("ID","Gene Symbol")]
colnames(ids) <- c("probe_id","symbol")
dim(ids)
ids=ids[ids$symbol != '',]
dim(ids)
ids=ids[ids$probe_id %in% rownames(dat),]
dim(ids)
ids$probe_id=as.character(ids$probe_id)
dat[1:4,1:4]
dat=dat[ids$probe_id,]
dim(dat)
#
if(T){
#探针基因ID对应以及去冗余
dat[1:4,1:4]
ids=ids[ids$probe_id %in% rownames(dat),]
dat[1:4,1:4]
dat=dat[ids$probe_id,]
probe_exp = dat
probe_info = ids
#ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids$median=apply(dat,1,median)
#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
ids=ids[!duplicated(ids$symbol),]
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4] #保留每个基因ID第一次出现的信息
}
save(gse_number,dat,group_list,
#probe_exp,probe_info,pd,
file = 'step1_output.Rdata')3 GSE92415
代码语言：r复制rm(list = ls()) ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
library(AnnoProbe)
library(GEOquery)
library(ggplot2)
library(ggstatsplot)
library(patchwork)
library(reshape2)
library(stringr)
getOption('timeout')
options(timeout=10000)
#获取并且检查表达量矩阵
##～～～gse编号需修改～～～
gse_number <- 'GSE92415'
#dir.create(gse_number)
#setwd( gse_number )
#list.files()
if(T){gset <- geoChina(gse_number)}
#gset = getGEO("GSE12366", destdir = '.', getGPL = F)
gset[[1]]
a=gset[[1]]
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
pd = pData(a)
head(pd)
##～～～查看数据是否需要log～～～
boxplot(dat[,1:4],las=2)
range(dat)
#dat=log2(dat+1)
boxplot(dat,las=2)
library(limma)
dat=normalizeBetweenArrays(dat)
a=gset[[1]]
#根据生物学背景及研究目的人为分组
##通过查看说明书知道取对象a里的临床信息用pData
pd=pData(a)
##挑选一些感兴趣的临床表型。
colnames(pd)
#
library(stringr)
##～～～分组信息编号需修改～～～
group_list=ifelse(grepl("Healthy", pd$characteristics_ch1.2,ignore.case = T ),
"control","case")
table(group_list)
group_list <- factor(group_list, levels = c("control","case"))
table(group_list)
a
a@annotation
# #获取芯片注释信息
# ##～～～芯片编号需修改～～～
ids=idmap( a@annotation ,'soft')
head(ids)
#上述函数用不了
library(data.table)
tmp <- fread("../GSE87473/GPL13158.txt", sep = "\t", skip = 18, data.table = FALSE)
tmp <- tmp[-nrow(tmp), ]
dim(tmp)
ids <- tmp[,c("ID","Gene Symbol")]
colnames(ids) <- c("probe_id","symbol")
dim(ids)
ids=ids[ids$symbol != '',]
dim(ids)
ids=ids[ids$probe_id %in% rownames(dat),]
dim(ids)
ids$probe_id=as.character(ids$probe_id)
dat[1:4,1:4]
dat=dat[ids$probe_id,]
dim(dat)
#
if(T){
#探针基因ID对应以及去冗余
dat[1:4,1:4]
ids=ids[ids$probe_id %in% rownames(dat),]
dat[1:4,1:4]
dat=dat[ids$probe_id,]
probe_exp = dat
probe_info = ids
#ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids$median=apply(dat,1,median)
#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
ids=ids[!duplicated(ids$symbol),]
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4] #保留每个基因ID第一次出现的信息
}
save(gse_number,dat,group_list,
#probe_exp,probe_info,pd,
file = 'step1_output.Rdata')4 Remove Batch Effects
代码语言：r复制rm(list = ls())
options(stringsAsFactors = F)
#数据导入
load('GSE38713/step1_output.Rdata')
GSE38713_dat <- dat
GSE38713_g <- group_list
table(GSE38713_g)#control 13 /case 30
load('GSE87473/step1_output.Rdata')
GSE87473_dat = dat
GSE87473_g = group_list
table(GSE87473_g)#control 21 /case 106
load('GSE92415/step1_output.Rdata')
GSE92415_dat = dat
GSE92415_g = group_list
table(GSE92415_g)#control 21 /case 162
#数据合并
tmp <- intersect(rownames(GSE38713_dat),
rownames(GSE87473_dat))
ids <- intersect(tmp,
rownames(GSE92415_dat))
merge_dat = cbind(GSE38713_dat[ids,],
GSE87473_dat[ids,],
GSE92415_dat[ids,])
dim(merge_dat)
meta=data.frame(
gse_number = c(rep('GSE38713',length(GSE38713_g)),
rep('GSE87473',length(GSE87473_g)),
rep('GSE92415',length(GSE92415_g))),
group = c(GSE38713_g,GSE87473_g,GSE92415_g)
)
table(meta)
#数据整合去除批次效应
library(sva)
table(meta$gse_number)
combat_edata1 = ComBat(dat=as.matrix(merge_dat),
batch=meta$gse_number, mod=NULL,
par.prior=T, prior.plots=F)
combat_edata1[1:4,1:4]
boxplot(combat_edata1)
dat=combat_edata1
group_list=meta$group
boxplot(combat_edata1,col=as.numeric(as.factor(meta$gse_number)))
save(dat,group_list,
#probe_exp,probe_info,pd,
file = 'step1_output.Rdata')可以看到整合的control组和case组是与文章完全吻合的
整合后的表达矩阵
3 差异分析和火山图可视化
3.1 差异基因
参照文章的阈值 logFC = 0.3，adjust.p.value = 0.05，采用limma包进行分析。
rm(list = ls())
load(file = "step1_output.Rdata")
#差异分析
library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
deg = topTable(fit,coef = 2,number = Inf)
#3.加change列,标记上下调基因
logFC_t = 0.3
p_t = 0.05
k1 = (deg$adj.P.Val < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$adj.P.Val < p_t)&(deg$logFC > logFC_t)
deg = mutate(deg,change = ifelse(k1,"Down",ifelse(k2,"Up","Not")))
table(deg$change)可以看到分析出来的差异基因的数目和原文章有些许差异，这是极为正常的现象，因为你的分析流程和作者的不可能完全一样，但是这无伤大雅，不同的分析流程只会影响到那些介于是与不是的差异基因，不会影响到本身效果比较显著的差异基因。
3.2 火山图可视化
这里力求和原图一模一样的效果，所以这里调试了许久。特别是用EnhancedVolcano的包，想让上调和下调基因显示不同的颜色还挺困难的。
代码语言：r复制library(EnhancedVolcano)
#原文章图中显示的基因
genes <- c("AQP8", "ABCG2", "HMGCS2", "PCK1", "CLDN8",#下调
"SLC6A14", "DUOX2", "CXCL1", "MMP3", "DUOXA2", "MMP10", "S100A8",
"TNIP3", "SAA1", "REG1A", "DEFA5", "REG1B")
# 给不同状态的颜色加标签
keyvals <- rep('#727271', nrow(deg))
# set the base name/label as 'Mid'
names(keyvals) <- rep('Not', nrow(deg))
# fold change > 1.5 & p-value < 0.0001 为高表达
keyvals[which(deg$logFC > 0.3 & deg$adj.P.Val<0.05)] <- '#38559b'
names(keyvals)[which(deg$logFC > 0.3 & deg$adj.P.Val<0.05)] <- 'Up'
# fold change < -1.5 & p-value < 0.0001为低表达
keyvals[which(deg$logFC < -0.3 & deg$adj.P.Val<0.05)] <- '#D8AE45'
names(keyvals)[which(deg$logFC < -0.3 & deg$adj.P.Val<0.05)] <- 'Down'
volcano_plot<- EnhancedVolcano(deg,
lab = rownames(deg),
x = 'logFC',
y = 'adj.P.Val',
pCutoff = 0.05,
FCcutoff = 0.3,
xlim = c(-5, 5),
ylim = c(0,80),
title = NULL,
subtitle = NULL,
caption = NULL,
xlab = 'logFC',
ylab = '-log10(adj.P.Val)',
pointSize = 2,#点大小
colAlpha = 1,#透明度
selectLab = genes,
labSize = 3,#标签的大小 之前的版本transcriptLabSize
colCustom = keyvals,
cutoffLineType = 'twodash',
cutoffLineWidth = 0.8,
legendPosition = 'right',
drawConnectors = TRUE,
boxedLabels = T,#是否加框
)
ggsave("volcano_plot.pdf", plot = volcano_plot, width = 8, height = 6)最后做出来的效果图如下
对比下两张图，可以看到两张图基本一致，但是原图（右图）还是有点小问题的，原文中logFC的阈值设置为0.3，但是从火山图中来看，明显logFC的阈值不是0.3，可能是1左右，还有就是原文中是采用adjust.p.value才call差异基因的，但是火山图的纵坐标是pValue。这里应该是作者画图的时候的错误。
本文标签： 文章复现bulkRNA转录组结合机器学习等进行相关疾病研究01多数据集去除批次效应后联合分析以及火山图标准绘制
版权声明：本文标题：文章复现—bulkRNA转录组结合机器学习等进行相关疾病研究01—多数据集去除批次效应后联合分析以及火山图标准绘制 内容由林淑君副主任自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://www.xiehuijuan.com/baike/1754624905a1704512.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论