8
$\begingroup$

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.

Improve this question
$\endgroup$

4 Answers 4

Reset to default
4
$\begingroup$

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).

Improve this answer
$\endgroup$
2
  • $\begingroup$ See my clarification. - I want to look at the spatial patterns within genes, not between them. $\endgroup$
    – Ian Sudbery
    May 18, 2017 at 8:44
  • $\begingroup$ @IanSudbery Ah, good to know, my answer isn't applicable then. $\endgroup$
    – Devon Ryan
    May 18, 2017 at 9:02
3
$\begingroup$

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 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 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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Ian Sudbery
    May 18, 2017 at 8:47
  • $\begingroup$ @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. $\endgroup$
    – CloudyGloudy
    May 18, 2017 at 16:16
  • $\begingroup$ 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. $\endgroup$
    – Ian Sudbery
    May 22, 2017 at 12:25
  • $\begingroup$ 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) $\endgroup$
    – CloudyGloudy
    May 22, 2017 at 15:37
3
$\begingroup$

Rather than working on the base level, you could probably work on say gene level counts. Kendall's tau, 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$
Improve this answer
$\endgroup$
5
  • $\begingroup$ Not sure why tex rendering is off. $\endgroup$
    – rightskewed
    May 17, 2017 at 18:13
  • $\begingroup$ 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? $\endgroup$
    – Ian Sudbery
    May 18, 2017 at 8:46
  • $\begingroup$ 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. $\endgroup$
    – rightskewed
    May 19, 2017 at 10:34
  • $\begingroup$ No, I'm trying to determine if two proteins show similar within gene binding patterns. $\endgroup$
    – Ian Sudbery
    May 19, 2017 at 11:57
  • $\begingroup$ Mathjax support hadn't been activated. It now has been so your answer renders correctly. $\endgroup$
    – terdon ♦
    Jul 26, 2017 at 16:05
2
$\begingroup$

It is one of my favorite stories.

Drop a glance to StereoGene software, it for genomic track correlation, it described in a preprint.

You also can run MACS or another peak caller and estimate the correlation of two interval sets using the GenomtriCorr package.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.