Pearson Correlation test between two Chip-Seq data sets

Question

[dally]

I was wondering how to run a pearson correlation test between two chip-seq data sets. I have narrowpeak files for both which are essentially bed files.

Is there any bioconductor packages in R that can do this? I've tried a variety of bedtools, and programs such as CHANCE, but they haven't given me what I required.

I used MACS2 to call peaks.

I used Bowtie to align to hg19.

Thanks!

[karl.stamm]

Do you know what Pearson correlation is? It measures the collinearity between two one dimensional measures.

[dally]

I am only vaguely familiar with the test. I've seen it on a large amount of papers recently while reading up on a new project. Essentially what I am trying to do as quoted from a paper is "find the number of uniquely aligned reads per Pol II gene for two complete biological replicates generated from * insert cell type * cells (Pearson's correlation, R = 0.##)

[karl.stamm]

There's no test. There's no hypothesis, and there's no p-value. The correlation is just a measurement.

[Madelaine Gogol]

You can window the data along the genome and calculate a correlation for the values in each window, but you definitely want to filter low windows. This post had some additional ideas.

Answer 1

[Constantine]

Try the bamCorrelate function from the deeptools software

[Fidel]

We have updated to deepTools 2.0 and now bamCorrelate is called multiBamCoverage. The output of that program can be used to compute a correlation (using plotCorrelation) of a PCA (using plotPCA). The new documentation is here:

http://deeptools.readthedocs.org/en/latest/content/tools/multiBamCoverage.html

[dally]

Thank you, this was exactly what I needed.

Answer 2

[karl.stamm]

The "number of uniquely aligned reads per PolII gene" is really just a row of a table of integers.

Sample    Gene1    Gene2    Gene3
Rep1      10       20       30
Rep2      12       18       40

For that table, you can compute the correlation in R like

> cor(c(10,20,30),c(12,18,40))
[1] 0.9496528

and say they are 95% correlated.

So your first step is to generate the table, specifying gene column headers to help you count reads at each site.

[dally]

This makes much more sense. Is there any software that will generate this table for me? Or is this a sort of "gotta build your own custom script" thing? I do mostly wet lab work, so something with custom scripts is probably out of my ballpark.

[karl.stamm]

There's a lot of different Biostars questions asked and answered for counting reads. Look for "how to count reads". The most popular tool is probably "htseq-count"

I like coverageBed from bedtools. You give it the BAM and a BED annotating your regions of interest, and it spits out how many reads were in each site.