<p>Honestly I'd just use <code>multiBigwigSummary</code> and then <code>plotCorrelation</code> from deepTools for this, but I'm a bit biased. There, the idea would be to consider each gene as a unit (you could instead use bins, but I don't think that would as nicely do what you want), namely by giving the tools a BED or GTF file input. It would then calculate the average signal in each gene/transcript and you could do your spearman's correlation. Features with 0 in all samples could optionally be removed (<code>plotCorrelation --skipZeros</code>).</p>

<p>While you certainly could go the whole 9 yards and use per-base comparisons, that seems a bit overkill and I suspect that it won't really yield appreciably more information (especially when one considers the additional time overhead).</p>

CC BY-SA 3.0
2017-05-18T08:44:44.850Z
See my clarification. - I want to look at the spatial patterns within genes, not between them.
Ian Sudbery
CC BY-SA 3.0
2017-05-18T09:02:34.753Z
@IanSudbery Ah, good to know, my answer isn't applicable then.
Devon Ryan
CC BY-SA 3.0
2017-05-17T18:04:16.737Z
2017-05-17T18:04:16.737Z
Devon Ryan
<p>Rather than working on the base level, you could probably work on say gene level counts. <a href="https://en.wikipedia.org/wiki/Kendall_rank_correlation_coefficient" rel="nofollow noreferrer">Kendall's tau</a>, an ordinal association metric, can then be used as an appropriate correlation measure.</p>

<p>If $X$ and $Y$ are your iCLIP replicates, $i$ represents gene index and $(x_i, y_i)$ represents the number of RBP binding sites in $X$ and $Y$ respectively for the $i^{th}$ gene, Kendall's tau is defined as :</p>

<p>$$
\frac{\text{#(concordant pairs)} - \text{#(discordant pairs)}}{n(n-1)/2}
$$</p>

<p>Where any two pairs $(x_i, y_i)$ and $(x_j, y_j)$ are concordant if:</p>

<ul>
<li>$x_i &gt; x_j$ AND $y_i &gt; y_j$</li>
</ul>

<p>OR</p>

<ul>
<li>$x_i &lt; x_j$ AND $y_i &lt; y_j$</li>
</ul>

<p>Correspondingly they are discordant if:</p>

<ul>
<li>$x_i &lt; x_j$ AND $y_i &gt; y_j$</li>
</ul>

<p>OR</p>

<ul>
<li>$x_i &gt; x_j$ AND $y_i &lt; y_j$</li>
</ul>

CC BY-SA 3.0
2017-07-26T16:05:04.010Z
Mathjax support hadn't been activated. It now has been so your answer renders correctly.
terdon
CC BY-SA 3.0
2017-05-17T18:13:02.973Z
Not sure why tex rendering is off.
Saket Choudhary
CC BY-SA 3.0
2017-05-18T08:46:13.987Z
See my clarification: I'm interested in the patterns within genes, rather than between them. Is tau better at dealing with data where many observations are zero?
Ian Sudbery
CC BY-SA 3.0
2017-05-19T10:34:29.080Z
If your interest is to study the pattern within genes, does your second sample act as a replicate? I am not sure I completely understand what "within genes" would mean here.
Saket Choudhary
CC BY-SA 3.0
2017-05-19T11:57:07.337Z
No, I'm trying to determine if two proteins show similar within gene binding patterns.
Ian Sudbery
CC BY-SA 3.0
2017-05-17T18:08:38.830Z
2017-07-26T16:30:47.080Z
2017-07-26T16:30:47.080Z
Saket Choudhary
<p>It depends whether you want to treat the peak intensities as binary (comparing presence/absence of peaks in the sets) or continuous (comparing the relative magnitudes of the peaks).</p>

<p><strong>Binary</strong></p>

<p>For starting out, a simple binary comparison may be appropriate. You can use a <a href="http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0096303" rel="nofollow noreferrer">peak caller</a> of your choosing to identify peaks in each sample according to your desired criteria. Then you can use a similarity metric such as the <a href="https://en.wikipedia.org/wiki/Jaccard_index" rel="nofollow noreferrer">Jaccard index</a> to quantify the level of agreement among the peaks in the two samples.</p>

<p>One potential obstacle is that defining the boundaries of your peaks won't be totally straightforward. For example, a peak in one sample might have 2 overlapping peaks in the other sample, one on each end. A rough solution for this is to divide the genome into bins (maybe around 100-1000 bp, depending on your desired resolution). You can treat a peak as present in a bin if more than half of the peak lies in the bin. That way, bins in one sample can be directly compared to the corresponding bins in the other sample. Obviously, this isn't the only way to do this; other appropriate methods exist too.</p>

<p><strong>Continuous</strong></p>

<p>If you want to treat the peak intensities as continuous, you could apply a similar binning method, taking the "score" of a bin to be the average peak intensity at positions within that bin. You could then throw away all bins with no peaks or only low-intensity peaks throughout the genome. Then you could compute the Spearman's correlation for the remaining bins. I'm guessing it will be harder to find a strong correlation for continuous intensities, because of the amount of experimental variability that is inherently present.</p>

<p>If, after following these steps, the Spearman's correlation is still "artificially low" as you suggested, then this is likely a problem with the underlying data, not the overall analysis; maybe your two datasets actually don't agree that well.</p>

CC BY-SA 3.0
2017-05-18T08:47:55.113Z
I was worried about artificially low values when most observations were 0. Excluding cases where observations are 0 in both signals just didn't "feel" right.
Ian Sudbery
CC BY-SA 3.0
2017-05-18T16:16:26.250Z
@IanSudbery I see. So you're also interested in the extent to which the "non-peaks" agree with one another. One possibility then is not to throw away the bins where both values are zero, and to report the Spearman correlation. However, since you said the peaks are quite sparse, your zero values will make it hard to see the extent of correlation between the true positives. I agree that it's not "right" to exclude the zero values in Spearman correlation scores, but if you instead use a metric for set intersections (like the Jaccard index) excluding zeroes would be okay.
CloudyGloudy
CC BY-SA 3.0
2017-05-22T12:25:48.907Z
Looking at the intersection solves the sparsity problem, but not the zero-inflated one: signals from similar looking tracks tend to be close to each other, but not right on top of each other. Honestly i'm starting to think your second suggest is the best there is.
Ian Sudbery
CC BY-SA 3.0
2017-05-22T15:37:29.640Z
That's true, but if you use the binning approach to define intersections instead of looking for strict overlap, this could help. (i.e. it counts as an intersection if the same bin contains peaks in both tracks, even if the peaks themselves don't directly overlap)
CloudyGloudy
CC BY-SA 3.0
2017-05-17T20:18:09.070Z
2017-05-17T20:18:09.070Z
CloudyGloudy
<p>It is one of my favorite stories. </p>

<p>Drop a glance to <a href="http://stereogene.bioinf.fbb.msu.ru/" rel="nofollow noreferrer">StereoGene software</a>, it for genomic track correlation, it described in a <a href="http://biorxiv.org/content/early/2017/05/25/059584" rel="nofollow noreferrer">preprint</a>.</p>

<p>You also can run MACS or another peak caller and estimate the correlation of two interval sets using the <a href="http://genometricorr.sourceforge.net/" rel="nofollow noreferrer">GenomtriCorr</a> package.</p>

CC BY-SA 3.0
2017-05-30T19:18:22.740Z
2017-06-15T18:16:11.983Z
2017-06-15T18:16:11.983Z
Sasha Favorov
<p>This question pertains to iCLIP, but it could just as easily be ChIP-seq or ATAC-seq or mutation frequencies. </p>

<p>I have iCLIP read counts across the transcriptome and I wish to know if the signals are correlated - that is, where one of them is high, is the other likely to be high.</p>

<p>Often when dealing with such data (e.g. iCLIP data) we know that the data is generally sparse - that is at most positions both signals are zero and this is correct, and also zero-inflated - that is many bases that "should" have a signal are missing that data. So just calculating the Spearman's correlation is likely to give an artificially low value.</p>

<p><strong>What might be a way to asses the association? I should add that the aim is the assess association of binding patterns within genes, rather than (or as well as) between genes.</strong></p>

<p>Things I have thought of:</p>

<ul>
<li>Apply some sort of smoothing to the data (eg a rolling mean). Remove any bases with 0 in both samples. Compute the spearmans. </li>
<li>Calculate the average pairwise distance between every read in sample one and every read in sample two. Compare this to data where the reads have been randomised within genes. </li>
</ul>

<p>In the first case removing all bases with 0 in both samples seems wrong. But if 99.99% of all bases have zero in both samples, then this seems like its necessary for Spearman. </p>

<p>In the second case, the result seems like it would be non-intuitive to interpret. And also calculating this would be massively computationally intensive. </p>

<r>
How to correlate two zero inflated bedgraph-like signals?
# How to correlate two zero inflated bedgraph-like signals?

(Asked by: username_0 on 2017-05-17)

This question pertains to iCLIP, but it could just as easily be ChIP-seq or ATAC-seq or mutation frequencies.

I have iCLIP read counts across the transcriptome and I wish to know if the signals are correlated - that is, where one of them is high, is the other likely to be high.

Often when dealing with such data (e.g. iCLIP data) we know that the data is generally sparse - that is at most positions both signals are zero and this is correct, and also zero-inflated - that is many bases that "should" have a signal are missing that data. So just calculating the Spearman's correlation is likely to give an artificially low value.

**What might be a way to asses the association? I should add that the aim is the assess association of binding patterns within genes, rather than (or as well as) between genes.**

Things I have thought of:

* Apply some sort of smoothing to the data (eg a rolling mean). Remove any bases with 0 in both samples. Compute the spearmans.
* Calculate the average pairwise distance between every read in sample one and every read in sample two. Compare this to data where the reads have been randomised within genes.

In the first case removing all bases with 0 in both samples seems wrong. But if 99.99% of all bases have zero in both samples, then this seems like its necessary for Spearman.

In the second case, the result seems like it would be non-intuitive to interpret. And also calculating this would be massively computationally intensive.

**Tags:** r

---

## 4 answers

### Answer by username_1 on 2017-05-17. Score: 4

Honestly I'd just use `multiBigwigSummary` and then `plotCorrelation` from deepTools for this, but I'm a bit biased. There, the idea would be to consider each gene as a unit (you could instead use bins, but I don't think that would as nicely do what you want), namely by giving the tools a BED or GTF file input. It would then calculate the average signal in each gene/transcript and you could do your spearman's correlation. Features with 0 in all samples could optionally be removed (`plotCorrelation --skipZeros`).

While you certainly could go the whole 9 yards and use per-base comparisons, that seems a bit overkill and I suspect that it won't really yield appreciably more information (especially when one considers the additional time overhead).

> **Comments:**
> **username_0**: See my clarification. - I want to look at the spatial patterns within genes, not between them.  
> **username_1**: @IanSudbery Ah, good to know, my answer isn't applicable then.  

---

### Answer by username_2 on 2017-05-17. Score: 3

Rather than working on the base level, you could probably work on say gene level counts. [Kendall's tau](https://en.wikipedia.org/wiki/Kendall_rank_correlation_coefficient), an ordinal association metric, can then be used as an appropriate correlation measure.

If $X$ and $Y$ are your iCLIP replicates, $i$ represents gene index and $(x\_i, y\_i)$ represents the number of RBP binding sites in $X$ and $Y$ respectively for the $i^{th}$ gene, Kendall's tau is defined as :

$$
\frac{\text{#(concordant pairs)} - \text{#(discordant pairs)}}{n(n-1)/2}
$$

Where any two pairs $(x\_i, y\_i)$ and $(x\_j, y\_j)$ are concordant if:

* $x\_i > x\_j$ AND $y\_i > y\_j$

OR

* $x\_i < x\_j$ AND $y\_i < y\_j$

Correspondingly they are discordant if:

* $x\_i < x\_j$ AND $y\_i > y\_j$

OR

* $x\_i > x\_j$ AND $y\_i < y\_j$

> **Comments:**
> **username_2**: Not sure why tex rendering is off.  
> **username_0**: See my clarification: I'm interested in the patterns within genes, rather than between them. Is tau better at dealing with data where many observations are zero?  
> **username_2**: If your interest is to study the pattern within genes, does your second sample act as a replicate? I am not sure I completely understand what "within genes" would mean here.  
> **username_0**: No, I'm trying to determine if two proteins show similar within gene binding patterns.  
> **username_3**: Mathjax support hadn't been activated. It now has been so your answer renders correctly.  

---

### Answer by username_4 on 2017-05-17. Score: 3

It depends whether you want to treat the peak intensities as binary (comparing presence/absence of peaks in the sets) or continuous (comparing the relative magnitudes of the peaks).

**Binary**

For starting out, a simple binary comparison may be appropriate. You can use a [peak caller](http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0096303) of your choosing to identify peaks in each sample according to your desired criteria. Then you can use a similarity metric such as the [Jaccard index](https://en.wikipedia.org/wiki/Jaccard_index) to quantify the level of agreement among the peaks in the two samples.

One potential obstacle is that defining the boundaries of your peaks won't be totally straightforward. For example, a peak in one sample might have 2 overlapping peaks in the other sample, one on each end. A rough solution for this is to divide the genome into bins (maybe around 100-1000 bp, depending on your desired resolution). You can treat a peak as present in a bin if more than half of the peak lies in the bin. That way, bins in one sample can be directly compared to the corresponding bins in the other sample. Obviously, this isn't the only way to do this; other appropriate methods exist too.

**Continuous**

If you want to treat the peak intensities as continuous, you could apply a similar binning method, taking the "score" of a bin to be the average peak intensity at positions within that bin. You could then throw away all bins with no peaks or only low-intensity peaks throughout the genome. Then you could compute the Spearman's correlation for the remaining bins. I'm guessing it will be harder to find a strong correlation for continuous intensities, because of the amount of experimental variability that is inherently present.

If, after following these steps, the Spearman's correlation is still "artificially low" as you suggested, then this is likely a problem with the underlying data, not the overall analysis; maybe your two datasets actually don't agree that well.

> **Comments:**
> **username_0**: I was worried about artificially low values when most observations were 0. Excluding cases where observations are 0 in both signals just didn't "feel" right.  
> **username_4**: @IanSudbery I see. So you're also interested in the extent to which the "non-peaks" agree with one another. One possibility then is not to throw away the bins where both values are zero, and to report the Spearman correlation. However, since you said the peaks are quite sparse, your zero values will make it hard to see the extent of correlation between the true positives. I agree that it's not "right" to exclude the zero values in Spearman correlation scores, but if you instead use a metric for set intersections (like the Jaccard index) excluding zeroes would be okay.  
> **username_0**: Looking at the intersection solves the sparsity problem, but not the zero-inflated one: signals from similar looking tracks tend to be close to each other, but not right on top of each other. Honestly i'm starting to think your second suggest is the best there is.  
> **username_4**: That's true, but if you use the binning approach to define intersections instead of looking for strict overlap, this could help. (i.e. it counts as an intersection if the same bin contains peaks in both tracks, even if the peaks themselves don't directly overlap)  

---

### Answer by username_5 on 2017-05-30. Score: 2

It is one of my favorite stories.

Drop a glance to [StereoGene software](http://stereogene.bioinf.fbb.msu.ru/), it for genomic track correlation, it described in a [preprint](http://biorxiv.org/content/early/2017/05/25/059584).

You also can run MACS or another peak caller and estimate the correlation of two interval sets using the [GenomtriCorr](http://genometricorr.sourceforge.net/) package.


---

