Coverage Comparison Of Different Mappers

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.

coverage

• 2,961 views • link

updated 12.9 years ago by Pierre Lindenbaum

• written 12.9 years ago by rob234king

1 answer

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

• 312 views • link

updated 6.8 years ago by Ram

• written 12.9 years ago by Pierre Lindenbaum

one-line command rocks !

• 0 views • link

written 12.9 years ago by Nicolas Rosewick

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

• 0 views • link

updated 6.8 years ago by Ram

• written 12.9 years ago by rob234king