TVP

# 关于批次效应矫正后出现负值

• 首先要了解一下什么叫批次效应
• 那么如何解决批次效应呢?
• limma 包中 removeBatchEffect 函数中出现负值问题
• 异常值的处理方法
• 结尾

### 首先要了解一下什么叫批次效应

a batch effect occurs when non-biological factors in an experiment cause changes in the data produced by the experiment. Such effects can lead to inaccurate conclusions when their causes are correlated with one or more outcomes of interest in an experiment. They are common in many types of high-throughput sequencing experiments, including those using microarrays, mass spectrometers

简单的总结一下:批次效应的来源有很多,例如,实验人员,温度,湿度,加入的药剂等等…由于批次效应的存在,世界上不可能存在两次一模一样的实验结果,如果有那就是造假(这句话是我瞎编的)

### 那么如何解决批次效应呢?

目前主流的方法有：ComBat方法(这个方法出现负值比较多)、替代变量分析法、距离加权判别法和基于比值的方法等

1. 构造一个合理的去批次模型
2. 带上批次信息和数据进行模型拟合
3. 将数据放入批次模型中进行校正

### limma 包中 removeBatchEffect 函数中出现负值问题

``````library(limma)
# 首先随机生成一个具有批次的表达谱
# a1,a2 为第一批次,a3,a4为第二批次
n = 20
a1 = rnorm(n,mean = 5,sd =1)
a2 = rnorm(n,mean = 5,sd =1)
a3 = rnorm(n,mean = 500,sd =1)
a4 = rnorm(n,mean = 500,sd =1)
dat = data.frame(a1,a2,a3,a4)
dat``````
``````##          a1       a2       a3       a4
## 1  4.325009 4.242443 500.7656 500.1623
## 2  3.036909 4.583820 499.0214 500.1241
## 3  6.467585 5.699004 501.6757 500.4596
## 4  6.952332 5.140634 500.2092 500.2292
## 5  5.842731 4.214526 499.4631 498.9063
## 6  6.367415 3.780075 499.9727 500.2209
## 7  5.716041 4.402080 499.6904 500.8299
## 8  3.473262 6.027185 501.0130 499.3186
## 9  4.373247 4.032372 499.7803 500.7890
## 10 5.008814 6.666476 499.6079 499.8026
## 11 5.370628 4.513518 500.6244 498.5864
## 12 3.733842 3.406336 499.9219 498.1257
## 13 6.012681 6.185878 501.4360 499.9811
## 14 3.848765 5.061008 499.6790 498.9002
## 15 5.424255 5.716941 499.5216 502.1032
## 16 4.400094 5.194406 500.9268 500.4789
## 17 4.424244 4.366474 501.5144 499.1714
## 18 5.843304 6.357230 499.5256 500.8641
## 19 4.956007 4.261948 500.3811 499.7551
## 20 6.384053 7.214116 501.0150 499.5818``````
``````library(pheatmap)
pheatmap(dat)``````

``````batch = c(1,1,2,2)
df = removeBatchEffect(dat,
batch = batch)
df``````
``````##             a1       a2       a3       a4
##  [1,] 252.4151 252.3326 252.6755 252.0722
##  [2,] 250.9181 252.4650 251.1402 252.2429
##  [3,] 253.9598 253.1912 254.1835 252.9674
##  [4,] 254.0387 252.2270 253.1228 253.1428
##  [5,] 252.9208 251.2926 252.3850 251.8283
##  [6,] 253.8789 251.2916 252.4612 252.7094
##  [7,] 253.3166 252.0026 252.0899 253.2294
##  [8,] 251.1810 253.7350 253.3052 251.6108
##  [9,] 252.4142 252.0733 251.7394 252.7481
## [10,] 251.9426 253.6003 252.6741 252.8688
## [11,] 252.7023 251.8452 253.2927 251.2547
## [12,] 251.4607 251.1332 252.1950 250.3989
## [13,] 253.3173 253.4905 254.1314 252.6764
## [14,] 251.2661 252.4783 252.2616 251.4828
## [15,] 253.0452 253.3379 251.9007 254.4823
## [16,] 252.3529 253.1472 252.9740 252.5261
## [17,] 252.3980 252.3402 253.5406 251.1977
## [18,] 252.8906 253.4045 252.4783 253.8168
## [19,] 252.6856 251.9915 252.6515 252.0255
## [20,] 253.1337 253.9638 254.2653 252.8322``````
``pheatmap(df)``

``````n = 20
a1 = rnorm(n,mean = 500,sd =1)
a2 = rnorm(n,mean = 500,sd =1)
a3 = rnorm(n,mean = 500,sd =100)
a4 = rnorm(n,mean = 500,sd =100)
dat = data.frame(a1,a2,a3,a4)
batch = c(1,1,2,2)
df = removeBatchEffect(dat,
batch = batch)
df``````
``````##             a1       a2       a3       a4
##  [1,] 537.7696 538.5662 569.9667 506.3691
##  [2,] 479.7175 479.1112 439.6990 519.1297
##  [3,] 590.1917 588.5127 569.1306 609.5738
##  [4,] 517.1115 517.0660 471.6131 562.5644
##  [5,] 502.9291 503.5741 474.4474 532.0558
##  [6,] 518.9635 516.6697 512.0909 523.5422
##  [7,] 496.8190 497.6665 558.2569 436.2286
##  [8,] 444.2579 446.6837 374.1672 516.7743
##  [9,] 485.5867 486.9465 494.4157 478.1175
## [10,] 547.5067 546.2181 497.8157 595.9091
## [11,] 529.6868 529.5018 594.4324 464.7562
## [12,] 487.9258 487.0325 501.5869 473.3714
## [13,] 460.6049 462.5335 484.3022 438.8362
## [14,] 492.7862 495.1724 391.6916 596.2669
## [15,] 558.8889 557.1967 557.0981 558.9875
## [16,] 487.6311 490.0558 463.7952 513.8916
## [17,] 504.2347 502.6942 480.5099 526.4190
## [18,] 504.6615 505.4736 463.3255 546.8096
## [19,] 481.0901 481.7881 491.7063 471.1720
## [20,] 507.4263 508.3922 567.8398 447.9787``````
``pheatmap(df)``

``````n = 20
a1 = rnorm(n,mean = 5,sd =1)
a2 = rnorm(n,mean = 5,sd =1)
a3 = rnorm(n,mean = 500,sd =100)
a4 = rnorm(n,mean = 500,sd =100)
a1[1] = 200
a4[2] = 1
dat = data.frame(a1,a2,a3,a4)
batch = c(1,1,2,2)
df = removeBatchEffect(dat,
batch = batch)
df``````
``````##             a1       a2        a3        a4
##  [1,] 396.6261 201.3172 343.85994  254.0834
##  [2,] 126.4910 124.8079 371.08459 -119.7857
##  [3,] 278.4899 279.1699 370.04624  187.6136
##  [4,] 198.5896 200.0464 210.75117  187.8849
##  [5,] 244.2498 245.1607 161.17401  328.2365
##  [6,] 256.2880 257.3001 258.18580  255.4023
##  [7,] 296.4463 295.4160 367.83188  224.0304
##  [8,] 133.5816 133.9541  95.22865  172.3071
##  [9,] 290.7292 292.2697 372.60820  210.3908
## [10,] 200.8624 201.0808 151.28795  250.6552
## [11,] 238.2652 238.0220 247.02427  229.2629
## [12,] 244.2496 245.2831 270.46833  219.0644
## [13,] 292.9604 294.6897 322.89647  264.7536
## [14,] 162.8317 162.3553 205.78186  119.4051
## [15,] 230.4572 233.2674 199.10435  264.6203
## [16,] 264.6253 265.3615 346.24899  183.7378
## [17,] 319.4633 320.2481 246.33326  393.3781
## [18,] 277.3527 279.1256 213.02937  343.4489
## [19,] 232.2332 232.0170 284.16990  180.0802
## [20,] 243.6998 243.0443 162.91349  323.8306``````
``pheatmap(df)``

``````n = 20
a1 = rnorm(n,mean = 5,sd =1)
a2 = rnorm(n,mean = 5,sd =1)
a3 = rnorm(n,mean = 5,sd =1)
a4 = rnorm(n,mean = 5,sd =1)
a5 = rnorm(n,mean = 5,sd =1)
a6 = rnorm(n,mean = 5,sd =1)
a7 = rnorm(n,mean = 500,sd =100)
a8 = rnorm(n,mean = 500,sd =100)
a9 = rnorm(n,mean = 500,sd =100)
a10 = rnorm(n,mean = 500,sd =100)
a11 = rnorm(n,mean = 500,sd =100)
a12 = rnorm(n,mean = 500,sd =100)
a1[1] = 20000000
a4[2] = 1
dat = data.frame(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12)
batch = c(1,1,1,1,1,1,2,2,2,2,2,2)
df = removeBatchEffect(dat,
batch = batch)
df``````
``````##                 a1            a2            a3            a4            a5
##  [1,] 1.833359e+07 -1666407.1479 -1666405.8406 -1666409.3761 -1666407.6376
##  [2,] 2.401895e+02      239.5788      242.0466      235.7371      238.9493
##  [3,] 2.475231e+02      246.3076      246.3042      247.2449      245.4299
##  [4,] 2.501823e+02      250.7809      249.2230      249.8569      249.3185
##  [5,] 2.063781e+02      206.3622      207.4780      204.7642      206.5903
##  [6,] 2.659123e+02      265.1276      265.0505      263.5802      265.7264
##  [7,] 2.380901e+02      239.6339      240.2970      238.1826      238.4551
##  [8,] 2.609592e+02      260.2369      259.4295      259.8068      263.0927
##  [9,] 2.413162e+02      241.6367      241.3342      241.3404      240.4079
## [10,] 2.861402e+02      283.9601      285.1293      285.7786      285.8991
## [11,] 2.636705e+02      263.9683      263.1181      263.1354      262.8506
## [12,] 2.436528e+02      242.4592      243.2073      241.5365      242.5469
## [13,] 2.618892e+02      264.8314      262.9322      262.7577      264.3643
## [14,] 2.777165e+02      276.6748      278.4512      275.5434      277.5528
## [15,] 2.369637e+02      235.4952      234.9260      236.1859      236.0994
## [16,] 2.591900e+02      258.1998      257.9451      260.3612      259.1091
## [17,] 2.570514e+02      256.7701      255.6074      257.2384      256.0814
## [18,] 2.618720e+02      261.1358      259.7916      260.5563      260.7569
## [19,] 2.438197e+02      243.6779      242.8511      242.1618      243.4374
## [20,] 2.644406e+02      264.7627      263.8011      262.1363      262.8141
##                  a6            a7           a8           a9          a10
##  [1,] -1666407.8218 1666921.64300 1.666932e+06 1666930.3327 1.666835e+06
##  [2,]      238.4606     216.76381 3.002159e+02     340.2846 2.561326e+02
##  [3,]      245.7600     307.63501 2.787366e+02     124.2921 2.050753e+02
##  [4,]      249.4805     377.31221 2.590395e+02     346.4215 2.565909e+02
##  [5,]      207.9340     169.87062 2.305704e+02     158.0066 3.963109e+02
##  [6,]      266.7230     239.98411 4.093580e+02     123.2362 2.565042e+02
##  [7,]      237.6206     268.73043 2.607454e+02     203.6789 2.165893e+02
##  [8,]      261.1724     180.42772 2.224725e+02     332.8927 2.577614e+02
##  [9,]      240.8256     -13.46822 2.840026e+02     309.0982 3.447213e+02
## [10,]      285.6554     112.91614 1.903634e+02     350.3533 4.788515e+02
## [11,]      263.0271     123.48691 2.895869e+02     303.2690 3.284093e+02
## [12,]      242.9073     222.26538 2.271866e+02     300.5197 9.125782e+01
## [13,]      262.7608     271.87566 5.131568e+01      82.0737 4.327190e+02
## [14,]      276.7248     372.87311 2.851697e+02     261.9885 4.083923e+02
## [15,]      236.8317     274.17676 1.904748e+02     228.5594 3.598058e+02
## [16,]      260.9709     279.28600 2.277301e+02     154.6868 2.784211e+02
## [17,]      255.1001     233.21750 1.369220e+02     141.7044 3.515268e+02
## [18,]      259.2409     345.50541 1.024226e+02     209.4460 3.846927e+02
## [19,]      244.2695     228.23103 1.311428e+02     240.3178 4.180164e+02
## [20,]      263.9173     154.92524 2.888070e+02     348.9575 3.746935e+02
##                 a11          a12
##  [1,] 1666971.88237 1666958.6084
##  [2,]     -10.05017     331.6151
##  [3,]     313.74916     249.0816
##  [4,]      87.99697     171.4811
##  [5,]     159.37026     125.3780
##  [6,]     191.19347     371.8440
##  [7,]     245.93847     236.5968
##  [8,]     211.97698     359.1663
##  [9,]     119.08776     403.4193
## [10,]     261.01627     319.0621
## [11,]     413.56616     121.4516
## [12,]     235.87191     379.2085
## [13,]     300.48787     441.0637
## [14,]     210.79181     123.4482
## [15,]     253.26687     110.2183
## [16,]     280.45769     335.1944
## [17,]     332.70539     341.7728
## [18,]     271.06342     250.2234
## [19,]     286.95991     155.5494
## [20,]     225.95455     188.5342``````
``pheatmap(df)``

### 异常值的处理方法

1. 使用 3σ或者 1.5IQR原则过滤异常值
2. log转换(这个方法可以把偏态数据进行拉回来)
3. sigmoid函数对数据进行压缩(这个方法适用于除了异常值后方差较小的数据)
4. 如果这个基因你压根就不关心直接删掉

0 条评论

• ### OSCA单细胞数据分析笔记12—Intergrating Datasets

无论是scRNA-seq，还是Bulk RNA-seq，批次效应都是一个很头疼的问题，如何有效地校正、并且正确地使用校正后的数据是很值得讨论的分析点。

• ### 你要的rmarkdown文献图表复现全套代码来了（单细胞）

强烈要求我们推荐纯粹的R语言的文献图表复现全套代码，其实很容易检索到，2020奶牛7月仅仅是单细胞高分（IF>9）文章就有一百多篇，全部的单细胞相关文章有六七百...

• ### 单细胞数据清洗的这5个步骤你会做吗？

写教程的话，我的优点仅仅是量大，坚持了七年多写了超1万篇教程。但实际上绝大部分都浮于表面，深度不够。

• ### 一文搞定高通量数据整合分析中批次效应的鉴定和处理

批次效应表示样品在不同的批次处理和测量时引入的与生物状态不相关的系统性的技术偏差。很多因素都可能导致批次效应的产生，如不同实验条件、不同操作者、不同公司的试剂、...

• ### SMNN：对单细胞数据进行批次校正

随着单细胞测序技术的成熟和测序成本的不断下降，产生了越来越多的单细胞数据。在整合来自多个批次的单细胞数据时，批次效应校正至关重要。

• ### 总被审稿人提起的多重假设检验校正是什么？

NGS系列文章包括NGS基础、在线绘图、转录组分析 （Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析 （ChIP-seq基本分...

• ### BatchBench比较scRNA批次矫正方法

当你的才华还撑不起你的野心时，请潜下心来，脚踏实地，跟着我们慢慢进步。不知不觉在单细胞转录组领域做知识分析也快两年了，通过文献速递这个栏目很幸运聚集了一些小伙伴...

• ### Cell 深度| 一套普遍适用于各类单细胞测序数据集的锚定整合方案

自北京大学汤富酬教授（当时为英国剑桥大学格登研究所(Gurdon Institute) Azim Surani实验室博士后）等人于2009年在Nature Me...

• ### 【Briefings in Bioinformatics】四篇好文简读-专题17

A network-based method for brain disease gene prediction by integrating brain co...

• ### 高通量数据中批次效应的鉴定和处理（六）- 直接校正表达矩阵

处理批次因素最好的方式还是如前面所述将其整合到差异基因鉴定模型中，降低批次因素带来的模型残差的自由度。但一些下游分析，比如数据可视化，也需要直接移除效应影响的数...

• ### 批次效应到底是个什么东东？

简单翻译一下的话，就是： 批次效应是在进行实验的时候附带产生了和实验结果没有关系的数据偏差。例如， 1. 一组实验在星期一进行一次而另一组在星期二进行，...

• ### 重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程 （原理、代码和评述）

NGS系列文章包括NGS基础、转录组分析 （Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析 （ChIP-seq基本分析流程）、...

• ### 单细胞转录组测序中的批次效应知多少？ （上）

当然，这样做的第一件问题就是数据来自不同样本，实验室和实验的数据并不能“很好地混合”。数据不能很好的联合分析主要是基于细胞注释和使用UMAP 和/或 tSNE ...

• ### 异质脑：自闭症谱系障碍患者自发连接模式畸变

来自以色列魏茨曼科学研究学院的Avital Hahamy等人在Nature neuroscience上发表文章，发现自闭症谱系障碍（Autism spectru...

• ### NC |SCALE准确鉴定单细胞ATAC-seq数据中染色质开放特征

SCALE全称是Single-Cell ATAC-seq analysis vie Latent feature Extraction, 从名字中就能知道这个软...

• ### 复制粘贴就能运行的100套R实战演练代码也有错误

前面整理了100多套R代码，因为时间跨度有点长，而且公众号写作后没办法修改，所以安排实习生进行代码审查，看看是不是确实复制粘贴就可以运行。

• ### 综述：高维单细胞RNA测序数据分析工具（上）

当你的才华还撑不起你的野心时，请潜下心来，脚踏实地，跟着我们慢慢进步。不知不觉在单细胞转录组领域做知识分析也快两年了，通过文献速递这个栏目很幸运聚集了一些小伙伴...

• ### 超详细 | 生物医学研究和临床应用中scRNA-seq的数据分析指南

随着高通量scRNA-seq（包括临床样本）能力的扩大，对这些海量数据的分析能力已成为进入该领域研究人员的必备技能。近日，《Military Medical R...

10元无门槛代金券