表观调控13张图之二相关性热图看不同样本相关性

们已经公布了：6个小时的表观调控13张图视频课程免费大放送哦 其实很多朋友并没有留意到我们不仅仅是有视频，还有配套的学徒解读：

• 表观调控13张图之一证明基因干扰有效性

现在我们再解读一下第二张图，如果你对视频感兴趣，还是可以继续留邮箱，我们在圣诞节统一发邮件给大家全部的视频云盘链接和配套代码哈！

关于视频审查员

我把表观调控数据分析，拆分成为了13张图，分别录制为13个视频，即将免费发布在B站，这个期间我们的视频编辑师还在兢兢业业的奋斗，希望这13张图能带领大家学会表观调控数据分析的一般流程， 然后应用到自己的课题哈!同时我也招募了4位视频审查员和细节知识补充员。接下来我们就连载第一位视频审查员的13个笔记：

第一位视频审查员大家也许并不陌生了，早在2018-08-29我发布给学徒的ATAC-seq数据实战（附上收费视频） 他就学完了全部课程内容，还写了笔记在简书，大家搜索二货潜，就可以看到。

以下是正文

参禅之初，看山是山，看水是水；禅有悟时，看山不是山，看水不是水；禅中彻悟，看山仍然山，看水仍然是水。

——青原行思

当我们拿到数据时候，除了前面的质控等分析外，我们一般需要查看样品内的重复性怎么样，一般目前市面上的 RNA-seq、ChIP-seq 测序样品内的相关性都能高达 0.9 以上。

我们可以通过两种策略来计算样品内的相关性

1、根据基因的表达量信息来计算样品之间的相关性，比如 RNA-seq 。

2、将全基因组等分 bin 的方法，然后计算每个 bin 里面的 reads 数, 然后通过均一化等过程，再对数据进行计算相关性, 比如 ChIP-seq 等 DNA 类型测序数据。

热图一 通过基因的表达量来计算样品相关性

rm(list = ls())
options(stringsAsFactors = F)
a = read.table('../figure-01-check-gene-expression/all.counts.id.txt', header = T)
dim(a)
dat = a[,7:16]
library(stringr)
ac = data.frame(group=str_split(colnames(dat),'_',simplify = T)[,1])
rownames(ac) = colnames(dat)
M=cor(log(dat+1))
pheatmap::pheatmap(M,
                   annotation_col = ac,
                   # breaks = seq(0, 100, length.out = 100)
                   ) 

这里需要注意几点 cov()函数计算相关性有三种方法

参考来源：

• Pearson，Kendall和Spearman三种相关分析方法的异同
• 百度百科大家不会的百度就好，统计学概念有很多经典的解释，比如 StatQuest 的视频真的精彩，这里只做略微提示，就是告诉大家在使用一个函数之前一定一定一定要阅读它的 manual，不要一味的只顾 copy，甚至大家可以进一步去理解背后实现的原理。
• 《白话统计》

1

第1点

• pearson: 即我们所说的 皮尔逊相关系数，更加强调的是是否具有线性关系，如果样本数据点精确的落在直线上（计算样本皮尔逊系数的情况），或者双变量分布完全在直线上（计算总体皮尔逊系数的情况），则相关系数等于 1 或 -1。
• kendall: 肯德尔相关系数，接触的少。
• Spearman: 即我们所说的 斯皮尔曼相关系数, 又称 秩相关系数，是秩的排序后所处的位置的相关，往往侧重两者是正相关还是负相关。如果当 X 增加时，Y 趋向于增加，斯皮尔曼相关系数则为正。如果当 X 增加时，Y 趋向于减少，斯皮尔曼相关系数则为负。斯皮尔曼相关系数为零表明当 X 增加时 Y 没有任何趋向性。当 X 和 Y 越来越接近完全的单调相关时，斯皮尔曼相关系数会在绝对值上增加。

2

第2点

• 用《白话统计》中的话来说：线性相关系数小不等于没有相关性。有些线性相关关系系数小，但是其曲线（比如二次项）相关性较大。所以当我们发现相关系数较小的视化，最好通过散点图确定这是直线相关，否则相关系数小未必表示没有线性相关。得出的结论是能是 没有线性相关，但不能轻易说 没有相关。
• 同样是《白话统计》中的话：存在异常值的时候要谨慎对待相关性，具体这里就不赘述了。

热图二 分析deeptools软件的multiBigwigSummary和plotCorrelation得到的相关性结果

• linux 中运行：

multiBigwigSummary  bins -b  *.bw -o results.npz -p 4
plotCorrelation -in results.npz \
  --corMethod spearman --skipZeros \
  --plotTitle "Spearman Correlation of Read Counts" \
  --whatToPlot heatmap --colorMap RdYlBu --plotNumbers \
  --plotFileFormat pdf \
  -o heatmap_SpearmanCorr_readCounts.pdf   \
  --outFileCorMatrix SpearmanCorr_readCounts.tab

• 将结果导入 R

b = read.table('SpearmanCorr_readCounts.tab', header = T)
colnames(b) <- str_split(colnames(b),'\\.',simplify = T)[, 1]
rownames(b) <- str_split(colnames(b),'\\.',simplify = T)[, 1]
bc = data.frame(group = str_split(colnames(b),'_',simplify = T)[, 1])
rownames(bc) <- colnames(b)
pheatmap::pheatmap(b,
                   annotation_col = bc) 

我们可以很清楚的看到，样品内的重复性都是极高的。同一样品都聚类在一起。样品内的相关性显著高于样品间的相关性。说明数据重复性很好，可以进行进下一步。

最后想说点自己的感想：

在冯国双老师的《白话统计》开篇就有这样一句话：参禅之初，看山是山，看水是水；禅有悟时，看山不是山，看水不是水；禅中彻悟，看山仍然山，看水仍然是水。 当时我看这本书的时候，并没有意识到这句话的深刻含义，直到某一天我在怀疑自己的时候，一个朋友把这句话送给了我。再次强调不论我们是使用一个软件还是其他，我们一定要了解其背后的原理；软件是死的，但是人是活的。统计学是一定要学的，同一个东西，同一种方法，如果不注意背后的原理，有可能你得到的结果就是反的。希望我们都能从我以为是我以为→我以为的不是我以为→我以为的就是我以为 这个过程中走出来 。