This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Mapq filtered bigwigs have higher peaks than unfiltered bigwigs?

I have some ChIP-seq data that I am aligning to a region. To view the alignments in bigwig format on UCSC browser, I'm using a simple pipeline like:

bwa mem -t 2 -k 50 -c 1000000 ref.fasta reads1.fq reads2.fq >  align.sam
samtools view -b -h  -F 3852  align.sam |samtools sort -o align.sort.bam
samtools index align.sort.bam
bamCompare -b1 align.sort.bam -b2 input.sort.bam --operation ratio -o align.bw

Works great. I am also interested in viewing the alignments after filtering to remove any reads that could possibly align to another location (there are repetitive elements and I want to remove reads that are assigned randomly among candidate locations). I've done this by filtering based on mapq score and then generating bigwigs with the same pipeline:

samtools view -bq 20 align.sort.bam > align_mapq.bam
samtools index align_mapq.bam
bamCompare -b1 align_mapq.bam -b2 input.sort.bam --operation ratio -o align_mapq.bw

This also works. The weird thing comes in the two outputs. In a simple way, I understand that mapq filtering is removing a large part of the data. However, when I compare the two resulting bigwigs on the same scale in UCSC, there are some peaks where my protein is enriched that are higher in the mapq dataset than in the unfiltered dataset.

How can that be true if there is less data in the mapq-filtered dataset? There are no additional reads being mapped to any location--it's the same dataset, only smaller. Is there something happening in the normalization process? I'm wondering if this is due to the reduced number of reads being normalized to the same input with deeptools, which scales the samples to each other using read count. Maybe I need to specify to deeptools not to scale the mapq-filtered bigwig (--scaleFactorsMethod None), but I'm wondering if this will also impart scaling differences to the two bigwigs (normal and mapq_filtered)?

Any help would be appreciated! I'm at a loss to explain why this is happening at some peaks. I would like the mapq_alignment to be presented on the same scale as the original data, so the loss of data in the mapq_alignment can be visually interpreted (hope this makes sense).

bigwig normalization deeptools mapq

0 answers

No answers yet.

Log in to answer this question.