This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Normalising ChIP-seq data to RPM using a 1 or 2 count for paired-end data?

I want to normalise my paired-end ChIP-seq data (4 samples treated, 4 samples control) to reads per million (RPM). So far my calculation takes the number of reads mapped to a position, divides it by the total number of mapped reads in the sample, and then multiples this by one million. My question is, should the number of mapped paired-end reads be counted as 1 or 2 in the calculation? I believe Samtools reports the total number of mapped reads as half the actual number of paired-end reads which map?

My current pseudo script for normalising to RPM:

  • mappedReads = samtools view -c -F 4 input.bam
  • scalingFactor = 1000000 / mappedReads
  • bedtools genomecov -ibam input.bam -bg - scale scalingFactor -g chrom.sizes > output.bedGraph
  • bedGraphToBigWig input.bedGraph chrom.sizes output.bigWig
chip-seq normalisation paired-end
• 12,407 views • link
updated by Ram Ram • written by James Ashmore James Ashmore

For paired end data, rpm=fpm (fragments per million). So, it should be counted as one.

• 0 views • link
updated by Ram Ram • written by poisonAlien poisonAlien

Just to clarify, you are saying that each pair of reads should be counted as just one? So in a perfect world where 100% of my reads align, the number of total alignments would be half the total number of reads in the data?

• 0 views • link
updated by Ram Ram • written by James Ashmore James Ashmore

Yup. Count a pair of reads as one. For single end reads each read represents a fragment it is sequenced from whereas for paired-end, a pair of sequence represents the fragment it is sequenced from. As simple as that.

• 1 views • link
updated by Ram Ram • written by poisonAlien poisonAlien

3 answers

You can use deepTools to do the normalisation that you want. The output is a bigwig file. It is quite fast if you have multiple cores.

bamCoverage -b file.bam --normalizeUsingRPKM -o file.bw

The program first extends the read to match the paired-end length before computing the coverage. I opted to count paired end reads as 2 and not as 1 to avoid a bias when a read is not properly paired which could be a significant fraction.

• 0 views • link
updated by Ram Ram • written by Fidel Fidel

Thank you for the suggestion, I used Bowtie2 to map my reads but indicated to remove discordant and mixed alignments, in this case counting each pair of reads as 1 would be preferred

• 0 views • link
updated by Ram Ram • written by James Ashmore James Ashmore

One option is to filter your bam file to include only the first mate, then use it with bamCompare. If all your reads are properly paired this should work.

• 0 views • link
updated by Ram Ram • written by Fidel Fidel

I just pushed a commit to github enabling bamCoverage to select reads based on the SAM flag. Now you can do something like this to get what you want:

bamCoverage -b file.bam --normalizeUsingRPKM -o file.bw --samFlag 64

--samFlag 64 means to select only the first mate.

• 0 views • link
updated by Ram Ram • written by Fidel Fidel

I think each read should be counted as 1

or

you can try MAnorm for normalization if it suits for your analysis or comparisons

• 0 views • link
updated by Ram Ram • written by Manvendra Singh Manvendra Singh

Initially I just want to create a coverage track for the normalised data.

• 0 views • link
written by James Ashmore James Ashmore

In case you do not already know MACS2 will generate coverage plots for the non-redundant fragments used in the analysis (by default 1 unique fragment per strand) normalised to millions of reads using the --SPMR flag.

• 0 views • link
updated by Ram Ram • written by Ian Ian

I think MACS does not use the paired-end data to extend the reads. Does it extend the reads?

• 0 views • link
written by Fidel Fidel

MACS2 does not need to extend reads as it knows the start and end of the fragment by using the paired end reads information in the BAM file.

• 0 views • link
written by Ian Ian

For this you need to use the BAMPE mode of MACS2 which I am not sure it properly works. We tried it and didn't work, also I see numerous complains in the mailing list. Moreover, the bam file has to be sorted by name. I always use MACS2 for peak calling but not coverage tracks.

• 0 views • link
written by Fidel Fidel

Log in to answer this question.

More posts like this