4
$\begingroup$

Given a aligned bam file (wgs/bwa-mem) what steps do i need to perform to generate a plot as seen below (from this paper):

enter image description here

We are trying to see how much of genome was covered using pacbio and illumina for cow.


You mean remove this left plot? enter image description here

$\endgroup$

2 Answers 2

5
$\begingroup$

The general steps are:

  1. Install deepTools
  2. Run plotCoverage on your BAM files
  3. Remove the left-most plot in photoshop/gimp/imagemagick
$\endgroup$
2
  • $\begingroup$ updated my question $\endgroup$ Jan 11, 2018 at 0:36
  • $\begingroup$ Yes, I mean the left plot in the example, since it looks like you want % >= the coverage value on the X axis. $\endgroup$
    – Devon Ryan
    Jan 11, 2018 at 1:07
1
$\begingroup$

Devon gives a practical answer using the great deepTools package. These should be the theoretical, lower level steps behind the plots:

  • Compute read depth at each genomic position, using for example samtools mpileup with appropriate flags to accept/reject reads. You get a table like:

    Pos       Depth
    chr1:1    0
    chr1:2    10
    chr1:3    11
    ...
    chrX:100  10
    
  • For each unique value x in the Depth column count how many genomic positions have depth >= x. You can use Unix tools like sort ... | uniq -c | awk ... for this. You get a table like:

    Depth   N bases with depth >= x   % bases with depth >= x
    0       (all genomic positions)   100 %
    1       1000                      90 %
    2       995                       85 %
    .
    .
    .
    N       2                         0.1 %                         
    
  • Now just plot column Depth (x) vs % bases with depth >= x (y) with R or you favourite plotting tool.

Did I get it right?

$\endgroup$
1
  • $\begingroup$ yes, not sure bust does picar collectwgsmetrics do something similar? $\endgroup$ Jan 11, 2018 at 13:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.