venn | 多样本间peak重叠韦恩图的解决方案

日常瞎掰

  ChIP-seq、ATAC-seq等会获取到基因组范围内蛋白结合的peak区间，后续为了研究目的可能需要对不同条件下的样本做一个overlap，看看样本间peak重叠的情况。并不像做基因集间的overlap那样简单直接，也有很多工具可以用来绘制韦恩图。Peak集间的重叠情况从数据处理的角度来说就稍显复杂了。不过，好在还是有人造好了轮子，让咱们可以直接使用，下面来展示一下收集的两款好用的工具。

Intervene

  基于python的命令行工具，可用于peak间的重叠情况，主要包含两种功能，涉及三种可视化方式：

1. venn to compute Venn diagrams of up-to 6 sets
2. upset to compute UpSet plots of multiple sets
3. pairwise to compute and visualize intersections of genomic sets as clustered heatmap

  当然了，咱们今天主要展示一下如何使用Intervene做韦恩图。下面使用两个公共数据数据做演示一下：

head -n3 sample1_peaks.narrowPeak sample2_peaks.bed
==> sample1_peaks.narrowPeak <==
GL000009.2      200882  201351  GSM2661795_peak_1       93      .       6.46021 12.82362        9.31556 313
GL000194.1      183861  184100  GSM2661795_peak_2       33      .       3.66392 6.49670 3.34866 121
GL000205.2      60999   61212   GSM2661795_peak_3       37      .       3.26887 6.95261 3.77481 82

==> sample2_peaks.bed <==
chr1    629902  630062  peak1   12.19972
chr1    633920  634054  peak2   7.78683
chr1    906807  906941  peak3   7.78683

# 安装
pip install intervene

# 运行
intervene venn -i sample1_peaks.narrowPeak sample2_peaks.bed --save-overlaps --names sample1,sample2 --project sample1_sample2 --output sample_ovlap_intervene

# 结果
tree sample_ovlap_intervene
sample_ovlap_intervene
├── sample1_sample2_venn.pdf
└── sets
    ├── 01_sample2.bed
    ├── 10_sample1.bed
    └── 11_sample1_sample2.bed

结果如下：

  Intervene venn也可用于基因间的韦恩图，通过参数--type list控制，还有一些其他可以调整的参数，大家根据软件的帮忙选项查看吧，这里不再介绍了。
  生成的bed文件分别为各样本特有的区间和共有的区间，如10_sample1.bed表示sample1特有的peak，11_sample1_sample2.bed表示sample1、sample2共有的peak。

ChIPpeakAnno

  该软件是一个R包，功能也比较多，overlpap及韦恩图展示只是其中一个功能，也可用于注释peak、功能富集等。当然，今天咱们也仅展示一下如何获取多样本间peak重叠情况及韦恩图：

BiocManager::install("ChIPpeakAnno")
library(ChIPpeakAnno)

gr1 <- toGRanges("sample1_peaks.narrowPeak", format="narrowPeak") 
gr2 <- toGRanges("sample2_peaks.bed", format="BED")
ol <- findOverlapsOfPeaks(gr1, gr2)
makeVennDiagram(ol, NameOfPeaks=c('sample1','sample2'),
                fill=c("#009E73", "#F0E442"), 
                col=c("#D55E00", "#0072B2"), 
                cat.col=c("#D55E00", "#0072B2"))

# 获取特有和重叠的peak
names(ol$peaklist)
[1] "gr2"       "gr1"       "gr1///gr2"

ol$peaklist$gr1
GRanges object with 2365 ranges and 1 metadata column:
           seqnames               ranges strand |                  peakNames
              <Rle>            <IRanges>  <Rle> |            <CharacterList>
     [1] GL000009.2     [200882, 201351]      * |     gr1__GSM2661795_peak_1
     [2] GL000194.1     [183861, 184100]      * |     gr1__GSM2661795_peak_2
     [3] GL000205.2     [ 60999,  61212]      * |     gr1__GSM2661795_peak_3
     [4] GL000208.1     [ 29557,  29770]      * |     gr1__GSM2661795_peak_4
     [5] GL000213.1     [ 57524,  57786]      * |     gr1__GSM2661795_peak_5
     ...        ...                  ...    ... .                        ...
  [2361]       chrY [25381875, 25382099]      * | gr1__GSM2661795_peak_10046
  [2362]       chrY [56687821, 56688084]      * | gr1__GSM2661795_peak_10047
  [2363]       chrY [56723433, 56723672]      * | gr1__GSM2661795_peak_10048
  [2364]       chrY [56767239, 56767439]      * | gr1__GSM2661795_peak_10049
  [2365]       chrY [56857035, 56857278]      * | gr1__GSM2661795_peak_10050

结果如下：

  overlap和韦恩图是相互独立的操作，如果不想要特有和共有的peak，其实可以直接将gr1、gr2以列表的形式传入makeVennDiagram函数直接画韦恩图。findOverlapsOfPeaks函数返回结果是一个列表，其中peaklist为peak的重叠情况，如这里的gr1、gr2、gr1///gr2分别为两个样本特有和共有的结果，以GRanges形式存储。

结束语

  可以看出Intervene、ChIPpeakAnno使用起来都挺容易的，但后者的韦恩图效果还是更好看一些，尤其是根据数值大小调整重叠区域的面积。不过，虽然两个软件的结果大体一致，如特有的peak数目完全一样，但细究的话就会发现重叠的peak数目略微有些差别。当然，这点差异是可以忽略的，也是可以解释的，那么，这个现象的根本原因是什么呢？答案留给朋友们自己思考吧。

往期回顾

【海外招聘】美国威尔康奈尔医学院招聘博士后
R编程技巧 | 学习高手实现函数多功能化的两种方法
linux | while + read 实现本地 or 集群批处理，实用且优雅
pyscenic | 单细胞转录因子分析，原理图文详解
一网打尽scRNA矩阵格式读取和转化(h5 h5ad loom)