Question: Coverage Comparison Of Different Mappers
1
gravatar for rob234king
6.7 years ago by
rob234king600
UK/Harpenden/Rothamsted Research
rob234king600 wrote:

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 • 1.7k views
ADD COMMENTlink modified 6.7 years ago by Pierre Lindenbaum128k • written 6.7 years ago by rob234king600
2
gravatar for Pierre Lindenbaum
6.7 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum128k wrote:

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
ADD COMMENTlink modified 7 months ago by RamRS27k • written 6.7 years ago by Pierre Lindenbaum128k

one-line command rocks !

ADD REPLYlink written 6.7 years ago by Nicolas Rosewick8.8k

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

ADD REPLYlink modified 7 months ago by RamRS27k • written 6.7 years ago by rob234king600
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1668 users visited in the last hour