Coverage Comparison Of Different Mappers

Question

[rob234king]

Is there a tool or maybe a method in R to produce a line graph of the coverage of mapped resequencing, for each chromosome of the mapped bam file? I'm looking to compare the coverage across the genome of mapped reads using two different mapping software.

Answer 1

[Pierre Lindenbaum]

Do a simple loop with bash & samtools:

cut -d '       ' -f 1 reference.fa.fai |\
 while read CHROM
 do
     for BAM in *.bam
     do
         echo -n -e "${CHROM}\t${BAM}\t"
         samtools  view -F 4  -c ${BAM} ${CHROM}
     done
 done

[Nicolas Rosewick]

one-line command rocks !

[rob234king]

I use your command from above (shown below) but remove some of the spaces in -d ' ' part to get it to work but I don't think my output looks correct? Seems to give me the total reads mapped for each chromosome (note chromosome 0 are a series of contigs which may explain why repeats it) but I'm thinking now that it would be more useful if I could do something more detailed like average coverage every 100 bases and print figure and average position so I could evaluate the coverage between the two mappers over different regions. Couldn't see an option in samtools to refine this command, any suggestions would be welcomed.

cut -d ' ' -f 1 /home/rob/Desktop/E9/samtoolsE9/S_lycopersicum.fa |\
 while read CHROM
 do
     for BAM in *.bam
     do
         echo -n -e "${CHROM}\t${BAM}\t"
         samtools view -F 4  -c ${BAM} ${CHROM}
     done
 done > /home/rob/Desktop/coverage1.vcf

Below is the output I'm getting

SL2.40ch00    A_BWA.marked.bam    17419246
SL2.40ch00    A_novoalign.marked.bam    13106767
AATAATAATAATAATAATAATAATAAATAAATAAATAAATAATAATAATAATAATAATAATAAATAAATAAATAAATAAA    A_BWA.marked.bam    17419246
AATAATAATAATAATAATAATAATAAATAAATAAATAAATAATAATAATAATAATAATAATAAATAAATAAATAAATAAA    A_novoalign.marked.bam    13106767
TAAATAAATAAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATA    A_BWA.marked.bam    17419246
TAAATAAATAAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATA    A_novoalign.marked.bam    13106767
ATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAAAAATAATAATAATAATAATAATAATAATAAT    A_BWA.marked.bam    17419246
ATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAATAAAAATAATAATAATAATAATAATAATAATAAT    A_novoalign.marked.bam    13106767
AATCATCATCATCATCATCATCATCATCATCATCATCATAATAATAATAATCATAATAAGGATAGATGATCTTTTCAACT    A_BWA.marked.bam    17419246
AATCATCATCATCATCATCATCATCATCATCATCATCATAATAATAATAATCATAATAAGGATAGATGATCTTTTCAACT    A_novoalign.marked.bam

Thanks