Peak call summary

[NoemieL]

Hello,
I am interested in doing a vendiagramm showing the peaks overlap between my two groups, however I did not understand how you get the metadata(proj_all_4@peakSet)$PeakCallSummar ( and associated pdf) from the getPeakSet(proj_all_4) data. I expected to obtained 47168 peaks in the "No" group and I get 79705 peaks from getPeakSet(proj_all_4).

>metadata(proj_all_4@peakSet)$PeakCallSummar
                 Group     Var1   Freq
1           UnionPeaks   Distal 57.556
2           UnionPeaks   Exonic 16.252
3           UnionPeaks Intronic 91.620
4           UnionPeaks Promoter 23.621
5        No(n = 47168)   Distal 43.646
6        No(n = 47168)   Exonic 12.226
7        No(n = 47168) Intronic 66.538
8        No(n = 47168) Promoter 19.314
9  T(11;14)(n = 30860)   Distal 40.928
10 T(11;14)(n = 30860)   Exonic 14.111
11 T(11;14)(n = 30860) Intronic 72.367
12 T(11;14)(n = 30860) Promoter 21.776


>PeakSet=getPeakSet(proj_all_4)
>head(PeakSet)
GRanges object with 6 ranges and 13 metadata columns:
           seqnames        ranges strand |     score replicateScoreQuantile groupScoreQuantile Reproducibility
              <Rle>     <IRanges>  <Rle> | <numeric>              <numeric>          <numeric>       <numeric>
  T(11;14)     chr1 804685-805185      * |   23.6164                  0.585              0.466               2
  T(11;14)     chr1 817103-817603      * |  49.68177                  0.635              0.527               3
  T(11;14)     chr1 826645-827145      * |  16.05508                  0.493              0.348               2
  T(11;14)     chr1 827281-827781      * | 215.87407                   0.93              0.905               5
  T(11;14)     chr1 838249-838749      * |   23.6164                  0.585              0.466               3
  T(11;14)     chr1 857900-858400      * |  29.98427                   0.62              0.507               3
                      GroupReplicate distToGeneStart nearestGene    peakType distToTSS  nearestTSS        GC
                         <character>       <integer> <character> <character> <integer> <character> <numeric>
  T(11;14) T(11;14)._.                         12435      FAM87B      Distal         2  uc057aum.1    0.4012
  T(11;14)  T(11;14)._.                             17      FAM87B    Promoter        17  uc057aum.1     0.479
  T(11;14) T(11;14)._.                         626   LINC00115      Exonic        62  uc031tlo.2    0.5788
  T(11;14) T(11;14)._.                             8   LINC00115    Promoter         8  uc031tlo.2    0.6846
  T(11;14) T(11;14)._.                       10976   LINC00115    Intronic      6893  uc031tlo.2    0.4551
  T(11;14) T(11;14)._.                        18752      FAM41C      Exonic      6801  uc057aux.1    0.3812
                 idx         N
           <integer> <numeric>
  T(11;14)         1         0
  T(11;14)         2         0
  T(11;14)         3         0
  T(11;14)         4         0
  T(11;14)         5         0
  T(11;14)         6         0
  -------
>nrow(PeakSet2)
[1] 189049

>nrow(PeakSet2[which(grepl("No", PeakSet2$GroupReplicate)),])
[1] 79705

What type of filters did you performed? and how can I get the getPeakSet associated to the PeakCallSummary?
Thanks for you help

[rcorces]

I think there is some confusion about what the PeakSetSummary depicts.
No(n = 47168) means there were 47148 in the cells in the No cluster/group. That is not the number of peaks.
The peak sets generated by addReproduciblePeakSet() are not going to be useful for what you want. If you really want to do this, you would use the individual peak calls from each cluster/cell group which are in the ./PeakCalls folder.

However, its worth noting that peak calls are extremely sensitive to sequencing depth and data quality and comparing raw peak calls is often dubious.

[rcorces]

If you take a step back and think about how you would determine if an ATAC-seq peak is an enhancer without some orthogonal data (for ex. ChIP-seq), I think you will find that this is not possible. This is why we classify as distal, intronic, promoter, exonic instead.

[KuroToaster]

@rcorces I have a follow up question about this. I have a hard time of finding the number of peaks by types (e.g. distal, intronic, promoter, exonic). I really wants the number by peak types to produce an image by jupyter notebook.

[rcorces]

@KuroToaster - this information is stored in the peak set and obtained via getPeakSet().

[KuroToaster]

Thank you for the reply!
But how to do you get the number of peaks for each cluster by kind? For example, I would like to get the number of peak in the promoter region of cluster 1.

[rcorces]

look in the PeakCalls directory