Pearson correlation between replicates

[afreund]

Hi all,

I'm analyzing some ChIP-seq data, and I'd like to see how different my replicate samples are (pairwise comparison only is fine). To be clear, I want to compare the read density across the whole genome, not just the peaks. I know the process of comparing peaks has been covered in other threads, but I couldn't find a description of how to generate a simple scatterplot of read density at all genomic positions (after normalization), with the corresponding pearson correlation coefficient. It seems to be a common form of analysis, and it looks like people generally use ~200 bp windows. It also seems important to eliminate regions with no reads in either sample, to avoid artificially increasing the correlation. Any suggestions on how to tackle this would be much appreciated!

I'm dealing with FASTQ files, already mapped to hg19.

[dpryan]

Sounds like you'll be needing R and one of the packages from Bioconductor that can handle SAM/BAM files (unless you prefer to use the C or whatever API).

[afreund]

Thanks, that's definitely a good start, although a little more detail would be helpful. I have used R in the past, and I've looked around Bioconductor's website, but I'm still not sure how to go forward.

Does anyone know if Galaxy's pre-loaded SAMtools offer a way to do this?

[bioinfosm]

thats interesting. Why would you need that correlation? why not just use a statistical model, like the ones in useq or macs, and when you input the replicates as case and control, obtain not many significantly different regions.

[afreund]

Hi bioinfosm,

I'm fairly new to this, so correct me if I'm wrong, but although the two approaches are superficially similar (i.e. they both are a measure of the similarity between two samples), the Pearson correlation coefficient (PCC) is a much more robust measurement because it doesn't rely on an arbitrary significance cutoff, as well as being more quantitative, allowing for simpler comparisons. I'd like to assess the reproducibility of my samples between experiments, conditions, etc, and I believe it's easier to do that with the PCC. From what I've read, it's a fairly standard approach.

[kopi-o]

I think it's very easy to go wrong with simple Pearson correlations on genome-wide ChIP-seq tag count profiles. I would advise to go for something like this recent method: http://genomebiology.com/2012/13/3/R16/abstract

Note, I haven't tried it myself yet

[tir_al]

If you want to do genome wide correlation in R, you will need 2 bioconductor packages: RSamtools and GenomicRanges, and also the lengths for the chromosomes from the corresponding genome.

The pipeline should go something like this:

readBam -> convert bam to GRanges -> extend the reads (resize function) -> make the coverage (coverage function). Now you need to smooth the coverage a bit because, otherwise the resulting vector will be to big. This can be done easily using the first function in the Rle tips and tricks manual.
Repeat that for each sample, put all of the data in a data frame, and use the cor function on the resulting object.

I would advise you to go chromosome by chromosome (if you have an indexed bam file), and aggregate the data only before doing the correlation analysis.

If you want to do a comparison of peak regions, a very cool method was published in Bioinformatics recently: An effective statistical evaluation of ChIPseq dataset similarity.

Cheers!

[afreund]

Thanks very much, tir_al, I'll give that a try!

[archi]

Hi all, I'm also doing analysis of some nextgen data. right now,I'm handling fastq files of three experimental replicates and want to calculate their correlation(i.i correlation of the three replicates).
Can someone please guide on how to go about this in R using or without using Bioconductor.

I'm from a math background and know how to calculate pearson correlation coefficient between two set of numbers,but how to do that for fastq data??


thanks in advance.

[dpryan]

The short answer is that you don't.

Firstly, map the reads to the appropriate genome. Then perform whatever type of quantification you need (for ChIP-seq, this would be peak calling; for RNA-seq this would be counting reads per gene/transcript/whatever; etc.) and then calculate the correlation from the resulting metrics.