4 Analysis steps
A typical metagene analysis will consist steps:
- Extraction the read count of every BAM files in selected regions.
- Conversion in coverage.
- Noise removal
- Normalization of the coverage values.
- Matrices production.
- Data frame production.
- Generation of the metagene plot.
4.1 Minimal analysis
A minimal metagene analysis can be performed in 2 steps:
- Initialization (the
newfunction). plot
regions <- get_demo_regions()
bam_files <- get_demo_bam_files()
# Initialization
mg <- metagene$new(regions = regions, bam_files = bam_files)## Warning in parent.env(env): closing unused connection 6 (/tmp/RtmpWNPp3o/
## Rinst24b91f77cb75/metagene/extdata/list2.bed)## Warning in parent.env(env): closing unused connection 5 (/tmp/RtmpWNPp3o/
## Rinst24b91f77cb75/metagene/extdata/list1.bed)# Plotting
mg$plot(title = "Demo metagene plot")## align1_rep1_list1## align1_rep1_list2## align1_rep2_list1## align1_rep2_list2## align2_rep1_list1## align2_rep1_list2## align2_rep2_list1## align2_rep2_list2## ctrl_list1## ctrl_list2
As you can see, it is not mandatory to explicitly call each step of the metagene analysis. For instance, in the previous example, the plot function call the other steps automatically with default values (the next section will describe the steps in more details).
In this specific case, the plot is messy since by default metagene will produce a curve for each possible combinations of BAM file and regions. Since we have 5 BAM files and 2 regions, this gives us 10 curves.
If we want more control on how every step of the analysis are performed, we have to call each functions directly.
4.2 Complete analysis
In order to fully control every step of a metagene analysis, it is important to understand how a complete analysis is performed. If we are satisfied with the default values, it is not mandatory to explicitly call every step (as was shown in the previous section).
4.2.1 Initialization
During this step, the coverages for every regions specified are extracted from every BAM files. More specifically, a new GRanges is created by combining all the regions specified with the regions param of the new function.
regions <- get_demo_regions()
bam_files <- get_demo_bam_files()
mg <- metagene$new(regions = regions, bam_files = bam_files)## Warning in .Internal(as.vector(x, mode)): closing unused connection 6 (/
## tmp/RtmpWNPp3o/Rinst24b91f77cb75/metagene/extdata/list2.bed)## Warning in .Internal(as.vector(x, mode)): closing unused connection 5 (/
## tmp/RtmpWNPp3o/Rinst24b91f77cb75/metagene/extdata/list1.bed)4.2.2 Producing the matrices
To produce the metagene plot, the coverages must be converted in a matrix where the columns represent the positions and the rows the regions. Furthermore, to reduce the computation time during the following steps, the positions are also binned.
We can control the size of the bins with the bin_count argument. By default, a bin_count of 100 will be used during this step.
mg$produce_matrices()We can also use the design we produced earlier to remove background signal and combine replicates:
mg$produce_matrices(design = design)4.2.3 Producing the data.frame
The metagene plot are produced using the ggplot2 package, which require a data.frame as input. During this step, the values of the ribbon are calculated. The algorithm used for the estimation of the confidence interval can be specified with the stat parameter.
By default, metagene uses “bootstrap” to obtain a better estimation of the mean of enrichment for every positions in each group of regions/BAM files. Another approach, called “basic” will simply use the ribbon to represent 95% of the values for each positions in each groups.
mg$produce_data_frame(stat = "basic")## align1_list1## align1_list2## align2_list1## align2_list24.2.4 Plotting
During this step, metagene will use the data.frame to plot the calculated values using ggplot2. We show a subset of the regions by using the region_names and exp_names parameter. The region_names correspond to the names of the regions used during the initialization. The exp_name will vary depending if a design was added. If no design was added, this param correspond to the BAM name or BAM filenames. Otherwise, we have to use the names of the columns from the design.
mg$plot(region_names = "list1", title = "Demo plot subset")
