专栏首页生信技能树多种批次效应去除的方法比较

多种批次效应去除的方法比较

前面我在生信技能树推文:你确定你的差异基因找对了吗? 提出了文章的转录组数据的60个样品并没有按照毒品上瘾与否这个表型来区分,而是不同人之间的异质性非常高,这个时候我提出来了一个解决方案,就是理论上就可以把人当做是一个批次效应,使用sva包的combat函数,把这样的效应去除一下,接着再找差异。当然了,去除批次效应的方法,肯定不止这一个,现在让我们列举并且比较一下吧!

昨天我在生信菜鸟团也再次强调了批次效应的问题:批次效应不得不防

首先载入数据并且包装一个PCA可视化函数

rm(list = ls())  ## 魔幻操作,一键清空~
options(stringsAsFactors = F)

pca_plot = function(dddd,ggggg){
  library("FactoMineR")
  library("factoextra")
  df.pca <- PCA(t(dddd), graph = FALSE)
  fviz_pca_ind(df.pca,
               #axes = c(2,3),
               geom.ind = "point",
               col.ind = ggggg ,
               addEllipses = TRUE,
               legend.title = "Groups"
  )
}
# 下面的 step1-output.Rdata 文件,大家可以去学习我的GEO课程# 就知道如何制作啦。
load(file = 'step1-output.Rdata')
# 每次都要检测数据
dat[1:4,1:4]
table(group_list)
library(limma)
g=factor( group_list )
g
g=relevel(g,'con')
design=model.matrix(~g) 
fit=lmFit(dat,design) 
fit=eBayes(fit) 
options(digits = 4) 
topTable(fit,coef=2,adjust='BH') 
# 首先是瘾君子与正常人的差异分析
deg1=topTable(fit,coef=2,adjust='BH',number = Inf)
pca_plot(dat,g)

数据如下,是一个表达矩阵和分组信息,我在B站的GEO课程多次讲解了,大家读懂:

使用 limma 的 removeBatchEffect 函数

需要注意的是removeBatchEffect 函数这里表达矩阵和需要被去除的批次效应是必须参数,然后本来的分组也是需要添加进入,这样与真实分组相关的差异就会被保留下来。

load('trait.Rdata')
batch=trait$patient

## 使用 limma 的 removeBatchEffect 函数
dat[1:4,1:4]
ex_b_limma <- removeBatchEffect(dat,
                                batch = batch,
                                design = design)
dim(ex_b_limma) 
ex_b_limma[1:4,1:4]
fit=lmFit(ex_b_limma,design)  
fit=eBayes(fit) 
options(digits = 4) 
topTable(fit,coef=2,adjust='BH') 
deg2=topTable(fit,coef=2,adjust='BH',number = Inf)
pca_plot(ex_b_limma,g)

然后使用 sva 的 ComBat 函数

可能是我没有理解 sva 的 ComBat 函数用法,下面的代码发现并不能把个体差异抹去,无法凸显出毒品成瘾与否两个分组的差异。

## 使用 sva 的 ComBat 函数
library(sva)
table(batch)
table(g)
ex_b_sva = ComBat(dat=as.matrix(dat), 
                       batch=batch 
                  )
ex_b_sva[1:4,1:4]
fit=lmFit(ex_b_sva,design)  
fit=eBayes(fit) 
options(digits = 4) 
topTable(fit,coef=2,adjust='BH') 
# 首先是瘾君子与正常人的差异分析
deg3=topTable(fit,coef=2,adjust='BH',number = Inf)
pca_plot(ex_b_sva,g)

比较原始矩阵和去除批次效应后

可以看到,只有 limmaremoveBatchEffect 函数做到了把矩阵区分成为毒品上瘾与否的截然不同的两个部分。

毫无疑问,使用这样的去除了人的效应的表达矩阵后再做差异分析肯定是能找到非常多的有统计学显著效果的基因列表。单独PCA可视化 limmaremoveBatchEffect 函数后的表达矩阵,如下:

把3次差异分析结果结合起来

ids=rownames(deg1)
combDEG=cbind(deg1[ids,c(1,5)],
              deg2[ids,c(1,5)],
              deg3[ids,c(1,5)])

cor(combDEG[,c(1,3,5)])

那些之前被人群效应淹没的毒品上瘾相关的基因有:

然后我们去可视化这些基因,在两个表达矩阵:

代码很简单:

combDEG$probe_id=rownames(combDEG)
library(illuminaHumanv3.db)
ids=toTable(illuminaHumanv3SYMBOL)
deg=merge(combDEG,ids,by='probe_id')#通过merge函数,由于deg和ids都有probe_id这一列,因此通过'probe_id'合并为新的deg
deg=deg[order(deg[,4]),]
cg=c(head(deg[,1],100),
     tail(deg[,1],100))

# ex_b_limma, ex_b_sva, dat 
library(pheatmap)
n=t(scale(t(ex_b_limma[cg,])))  
n[n>2]=2 #限定上限,使表达量大于2的等于2
n[n< -2]= -2 #限定下限,使表达量小于-2的等于-2
n[1:4,1:4] 
load('trait.Rdata')
ac=data.frame(group_list=group_list)
rownames(ac)=colnames(n) 
pheatmap(n,show_colnames =F,show_rownames = F,annotation_col = ac)

现在的问题就是,我们定位的这些差异基因,是否在真正的两个组别的差异呢,还是仅仅是因为我们使用了算法抹去个体差异后的产物。

本文分享自微信公众号 - 生信技能树(biotrainee),作者:生信技能树

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-10-05

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 眼馋Y叔的可视化函数却还想使用DAVID结果肿么办

    实际上,我很能理解粉丝的心情, 确实就这么一个数据集,使用Y叔的clusterProfiler包的结果不满意,就不得不使用在线网页工具:Functional A...

    生信技能树
  • 论怎么黏在一起-R(paste|paste0)|Excel(concatenate)

    Excel里的一系列已经写好的函数,用起来之后,会不禁感叹,那么长一段时间,我们真的是暴殄天物了!!!

    生信技能树
  • TCGA数据库临床资料官方大全

    因为TCGA计划跨时太长,纳入研究的病人数量太多, 或多或少有点资料继续错误或者不完整,所以TCGA团队下功夫在计划结束后(April 2018)完整的系统性的...

    生信技能树
  • 浅谈集合数据结构之栈与队列的区别

    吾爱乐享
  • 2019第三季度网络威胁分析

    在前几个季度新安装包的数量明显下降,而第三季度的数据比上个季度增加了117067个包。

    FB客服
  • 如何判断一个数据分布是否是正太分布?

    在数据分析过程中,得到一组数据,在分析之前,通常需要判断数据是否符合正态分布与否,再决定下一步分析方法。那么,如何判断数据是否属于正太分布呢?

    小末快跑
  • Ubuntu20.04修改ip地址的方法示例

    今天在继续上次的办公协作平台实验的时候发现虚拟机因为之前是DHCP的方式获取的ip地址,而这次虚拟机开机后DHCP获取的地址变更了,因为网站上已经配置了固定ip...

    砸漏
  • 爬虫框架Scrapy(一)

    Absorb what is useful. Discard what is not. Add what is uniquely your own.

    小闫同学啊
  • 微信小程序接入百度OCR(身份证识别)

    支持对二代居民身份证正反面所有8个字段进行结构化识别,包括姓名、性别、民族、出生日期、住址、身份证号、签发机关、有效期限,识别准确率超过99%;同时支持身份证正...

    Kindear
  • Android开发(54) AIDL示例

    最近在考虑项目重构的时候,考虑将项目拆分成两个APK,一个用于数据服务,一个用于UI展示。 数据服务APK向自己编写APK提供数据,同时也可以向第三方提供数据...

    zhangyunfeiVir

扫码关注云+社区

领取腾讯云代金券