Reading in coverage data

To do this, we must first read in the coverage information from the .bam file. This can be done with the getCoverage function:

chr4_coverage <- getCoverage(
  regions_of_interest = "chr4", 
  bam_file = pasillaBamSubset::untreated3_chr4(),
  window_size = 350, smoothing_window_size = 400)

Here, we get the smoothed and windowed coverage for chr4.³ As we wanted the coverage for the entire chr4, we could simply make regions_of_interest = "chr4". However, we could also have supplied a GRanges describing that area instead. Whichever input is used, it is really important that the sequence names match exactly. For example, the following will t hrow an error, because there is no sequence named “4” or “Chr4” in the .bam file:

getCoverage(regions_of_interest = "4", 
            bam_file = pasillaBamSubset::untreated3_chr4(),
            window_size = 300, smoothing_window_size = 400)
#> Error in getCoverage(regions_of_interest = "4", bam_file = pasillaBamSubset::untreated3_chr4(), : Make sure all of the chromsomes in regions_of_interest are in
#>                 the bam file and spelled exactly the same as in the bam
getCoverage(regions_of_interest = "Chr4", 
            bam_file = pasillaBamSubset::untreated3_chr4(),
            window_size = 300, smoothing_window_size = 400)
#> Error in getCoverage(regions_of_interest = "Chr4", bam_file = pasillaBamSubset::untreated3_chr4(), : Make sure all of the chromsomes in regions_of_interest are in
#>                 the bam file and spelled exactly the same as in the bam

Plotting coverage

Now that we have the coverage data, we can plot the ideogram again using this information. To draw a ‘coverage rectangle’ we need to firstly specifiy the coverage_data to be used (as either a GRanges or a data frame) and then also supply to coverage_rectangle a vector of the sector names to plot the coverage data for.⁴

gmovizInitialise(ideogram_data = fly_ideogram_chr4_only, 
                 coverage_rectangle = "chr4", 
                 coverage_data = chr4_coverage,
                 xaxis_spacing = 30) 

As you can see, the chr4 ideogram rectangle is replaced with a line graph showing the coverage over the entire chromosome. The coloured area represents the coverage, allowing easy identification of high and low coverage areas.

Smoothing and windowing coverage data

When reading in the coverage data, there are two additional parameters window_size and smoothing_window_size that modify the values.

• window_size controls the window size over which coverage is calculated (where a window size of 1 is per base coverage. A larger window size will reduce the time taken to read in, smooth and plot the coverage. It will also remove some of the variation in the coverage, although this is not its primary aim. If you have more than 10-15,000 points, it is highly recommended to use a larger window size, as this will take a long time to plot.

• smoothing_window_size controls the window used for moving average smoothing, as carried out by the pracma package. It does not reduce the number of points and so offers no speed improvement (in fact, it increases the time taken to read in the coverage data). It does, however, reduce the variation to produce a smoother, more attractive plot.

For example, try running the following:

# default window size (per base coverage)
system.time({getCoverage(regions_of_interest = "chr4", 
                         bam_file = pasillaBamSubset::untreated3_chr4())})

# window size 100
system.time({getCoverage(regions_of_interest = "chr4", 
                         bam_file = pasillaBamSubset::untreated3_chr4(),
                         window_size = 100)})

# window size 500
system.time({getCoverage(regions_of_interest = "chr4", 
                         bam_file = pasillaBamSubset::untreated3_chr4(),
                         window_size = 500)})

Notice how going from the default window size of 1 (per base coverage) to a relatively modest window size of 100 dramatically reduces the time needed to read in the coverage data.

In terms of the appearance of the plot: (note: for speed, we will plot only a subset of the chromosome: from 70000-72000bp)

# without smoothing
chr4_region <- GRanges("chr4", IRanges(70000, 72000))
chr4_region_coverage <- getCoverage(regions_of_interest = chr4_region,
                          bam_file = pasillaBamSubset::untreated3_chr4())
gmovizInitialise(ideogram_data = chr4_region, coverage_rectangle = "chr4", 
                 coverage_data = chr4_region_coverage, custom_ylim = c(0,4))

# with moderate smoothing 
chr4_region_coverage <- getCoverage(regions_of_interest = chr4_region,
                          bam_file = pasillaBamSubset::untreated3_chr4(),
                          smoothing_window_size = 10)
gmovizInitialise(ideogram_data = chr4_region, coverage_rectangle = "chr4", 
                 coverage_data = chr4_region_coverage, custom_ylim = c(0,4))

# with strong smoothing
chr4_region_coverage <- getCoverage(regions_of_interest = chr4_region,
                          bam_file = pasillaBamSubset::untreated3_chr4(),
                          smoothing_window_size = 75)
gmovizInitialise(ideogram_data = chr4_region, coverage_rectangle = "chr4", 
                 coverage_data = chr4_region_coverage, custom_ylim = c(0,4))

Notice how adding smoothing dramatically improves the appearance of the plot. It also slightly reduces the time taken, because there are less extreme points. However, it does result in the loss of the finer detail of the coverage data. Thus, it is recommended that you play around with the values of smoothing_window_size and window_size and choose a value that is best suited to your own data.

Setting custom sector widths

One way to manipulate the width/size of the sectors is to specify a custom_sector_width (custom sector width) vector. This vector should be the same length as the number of sectors. For example:

gmovizInitialise(fly_ideogram, 
                 custom_sector_width = c(0.2, 0.2, 0.2, 0.2, 0.2, 0.1, 0.1))

Notice that the custom_sector_width vector had length 7, because this is how many sectors there are. custom_sector_width can also be used for the insertionDiagram and featureDiagram functions in the same way.

Zooming

Whilst it is quite easy to set custom sector widths when there are only a few sectors, it can be quite troublesome for entire genomes. Also, using this method loses the relative sizing of all sectors, potentially leading to misinterpretation.

We can solve this problem by using the zooming functionality of gmovizInitialise. Doing this is relatively easy, all we need to do is supply the names of sector(s) to zoom to the zoom_sectors argument:

gmovizInitialise(fly_ideogram, zoom_sectors = c("chr4", "chrYHet"),
                 zoom_prefix = "z_")

Now, chr4 and chrYHet are clearly visible alongside the rest of the sectors. Notice that chrYHet is still around 1/4 of the size of chr4, as is expected from their relative sizes (347038bp and 1351857bp, respectively). Also, all of the other chromosomes are still proportional. Another advantage of using the zooming is that the zoom_prefix applied to the start of the zoomed sector label makes it clear which sectors have been zoomed and which have not.