<p>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?</p>

<p>My current pseudo script for normalising to RPM:</p>

<ul>
	<li>mappedReads = samtools view -c -F 4 input.bam</li>
	<li>scalingFactor = 1000000 / mappedReads</li>
	<li>bedtools genomecov -ibam input.bam -bg - scale scalingFactor -g chrom.sizes &gt; output.bedGraph</li>
	<li>bedGraphToBigWig input.bedGraph chrom.sizes output.bigWig</li>
</ul>

You can use [deepTools][1] 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.

 [1]: http://deeptools.github.io/