I want to run a CNV call with breakdancer, so I tried two methods :
1) Using AWK to select only reads that have a mapping quality of 60 like that :
samtools view -h sample1-bam1.bam | awk 'substr($0,1,1) == "@" || $5 == 60 {print}'| samtools view -bS - > sample1-MQ60-bam1.bam
bam2cfg.pl -h sample1-MQ60-bam1.bam sample1-MQ60-bam2.bam > sample1.cfg
breakdancer_max sample1.cfg > sample1-MQ60.txt >> returns a file with found CNVs
2) In parallel on Breakdancer directly using -q flag
bam2cfg.pl -h sample1-1.bam sample1-2.bam > sample1-method2.cfg
breakdancer_max -q 60 sample1-method2.cfg > 1-MQ60.txt >> returns an empty file
breakdancer_max -q 59 sample1-method2.cfg > 1-MQ59.txt >> returns an empty file
Technically the two methods should return the same amount of CNV found .. why when I use -q flag on breakdancer this returns an empty file ??
•
written
by
madkitty
In the first code block, you're selecting reads with a mapping quality of exactly 60. In the second example, you're using reads that have a mapping quality of at least 60.
I'm not sure that completely explains your problem, but it's certainly somewhere to start. May I also ask why you're requiring a MQ of 60? That seems awfully high. 35 or 40 is a more typical cutoff for breakdancer, IIRC.
I think by default it's 35 on breakdancer. But I had horrible long and multiple CNVs called whenever I use the default option on breakdancer, you can see it here : http://i46.tinypic.com/iwjdar.jpg
Even with mapping quality of 60 of the first method the result is not that different than the image you can see. If I use method 2) I have an empty file ..
Log in to answer this question.
