Test Data

In this example we’ll be using the WGS example report used in the MultiQC documentation. You can find view the report in your browser here, and find the associated multiqc_data.json file here. Feel free to download this file and use it to follow along with this vignette.

In the rest of the vignette we will use the multiqc_data_path variable to indicate the path to this multiqc_data.json file. Feel free to set this variable to the path to this file on your system.

Setup

First, install the package using:

install.packages("TidyMultiqc")

Then if you want, you can load the package:

library(TidyMultiqc)

However, for the sake of this tutorial we will explicitly use namespaced functions, so we won’t be using library a great deal.

Basic Usage

The main entry point to the TidyMultiqc package is the function load_multiqc. A basic invocation of this function looks like this:

(note the arrow for scrolling through columns)

df = TidyMultiqc::load_multiqc(multiqc_data_path)
df

We’ve now generated a tibble (a kind of data frame), whose rows are samples in the QC report, and whose columns are QC data and metadata about these samples.

By default this function only returns the “general” statistics, which are the ones in the “General Statistics” table at the top of the MultiQC report. In TidyMultiqc, these statistics are all prefixed by general. We can also extract the “raw” statistics, which includes some fields normally hidden from the report. These statistics will have the prefix raw.<toolname>. where <toolname> is the QC tool used to calculate it.

TidyMultiqc::load_multiqc(multiqc_data_path, sections = 'raw')

Often you won’t care about fields like raw.qualimap_bamqc_genome_results.bam_file, the path to the original BAM file, but ‘raw’ at least provides this option.

You can also combine both general and raw sections by passing in a longer vector:

df_both = TidyMultiqc::load_multiqc(multiqc_data_path, sections = c('raw', 'general'))
ncol(df_both)
#> [1] 388

That’s a lot of columns!

Uses

This section will briefly talk about some downstream use-cases for this package.

Plotting

One use for this data frame is creating QC plots. For example, to visualise the duplication rate per sample:

library(magrittr)

df %>%
  ggplot2::ggplot(ggplot2::aes(x=metadata.sample_id, y=general.percent_duplication)) +
  ggplot2::geom_col()

Of course, this is basically just replicating a plot already in the MultiQC report, but now we can customise it how we like!

Hypothesis Testing

With all this data, we might also want to test a specific hypothesis! For example, we might want to test the hypothesis that the mean GC content is the same as the mean GC content in the human genome (41%). If we assume that GC content is normally distributed, we can do the following:

t.test(df$general.percent_gc, mu=41)
#> 
#>  One Sample t-test
#> 
#> data:  df$general.percent_gc
#> t = -1, df = 5, p-value = 0.3632
#> alternative hypothesis: true mean is not equal to 41
#> 95 percent confidence interval:
#>  40.40490 41.26176
#> sample estimates:
#> mean of x 
#>  40.83333

It seems that we cannot reject this hypothesis, so these may well be human samples!

Extracting Metadata

It may be the case that your samples have important metadata that you want in your data frame. For example, it is common for the sample names or file names to be composed of a number of metadata fields, and indeed the report we are working with has IDs such as P4107_1003, which is composed of two identifiers.

To include this metadata in our output, we need to provide the find_metadata argument to load_multiqc, which is a function that is called for each sample, and which returns a named vector of metadata fields for that sample. It also gets passed the entire parsed MultiQC JSON report, so the function can traverse the structure as it wants to extract metadata.

From the File Name

Here is an example that parses the input file names to annotate additional metadata. Notice that the first argument our function is passed is a string which is the sample identifier for a sample, and how in this example we ignore the parsed argument. Also notice that the names we give to our return value (“batch” and “sample”) are prefixed by “metadata” to become the final names in the data frame.

TidyMultiqc::load_multiqc(
  multiqc_data_path, 
  find_metadata = function(sample, parsed) {
    # Split the sample ID to obtain some metadata
    segments <- stringr::str_split(sample, "_")[[1]]
    c(
      batch = segments[[1]],
      sample = segments[[2]]
    )
  }
)

From the File Path

We can extend this approach, but this time actually look up the file paths within the report_data_sources section of the MultiQC report:

TidyMultiqc::load_multiqc(
  multiqc_data_path, 
  find_metadata = function(sample, parsed) {
    # This gives us the path to the fastqc output file
    filepath = parsed$report_data_sources$FastQC$all_sections[[sample]]
    # Split into path segments
    path_segments = stringr::str_split(filepath, "/")[[1]]
    # The filename is the last path segment
    filename = dplyr::last(path_segments)
    # Split the filename using dots and underscores
    name_segments = stringr::str_split(filename, "[_\\.]")[[1]]
    # Arbitrarily assign names for the outputs
    name_segments %>% purrr::set_names(LETTERS[1:length(name_segments)])
  }
)

Of course in a real application we would choose specific names for each field.

From MultiQC Configuration

Finally, we might want to include metadata that doesn’t relate to the sample at all. For example, MultiQC has a number of report fields prefixed by config_ that we might want to store:

TidyMultiqc::load_multiqc(
  multiqc_data_path, 
  find_metadata = function(sample, parsed) {
    parsed[c(
      "config_creation_date",
      "config_version"
    )]
  }
)

Extracting Plot Data

Motivation

It is occasionally useful to extract QC data from the MultiQC plots. For example, let’s say we want to calculate the median quality score of every base in each sample. Unfortunately, MultiQC provides no numerical summary statistic for the read quality, it only has mapping quality and pass/fails for the per-base sequence quality:

df_both %>% dplyr::select(dplyr::contains('quality'))

However, our MultiQC report does have plots that contain this information. In particular, let’s look at the “Per Sequence Quality Scores” plot.

Listing Plots

Firstly, we need the ID of the plot we want. This isn’t necessarily obvious from just looking at the report, so we can use a utility function here:

TidyMultiqc::list_plots(multiqc_data_path)

Now, we know we want the “Per Sequence Quality Scores” plot, and by looking at the data frame above we can tell that the corresponding ID is fastqc_per_sequence_quality_scores_plot.

Loading Plots

Now that we have the plot ID, we can load the plot data. First, we need to tell TidyMultiqc to load include some plots by using sections = "plot" (you can load other sections at the same time, as explained above). Also, we need to pass the plot ID from the previous step into the plots argument:

df = TidyMultiqc::load_multiqc(
  multiqc_data_path, 
  sections = 'plot',
  plots = "fastqc_per_sequence_quality_scores_plot"
)
df

We now have the plot data, but it’s not in a very usable form! This is because each sample has an entire data frame of plot data. At this point if you’re comfortable using dplyr and tidyr to deal with nested data frames, you probably know what to do. Otherwise, read on.

Converting Plot Data

Recall that we are after the median quality score of each sample. First, we should look at the plot data for a single sample to know what we’re dealing with:

df$plot.fastqc_per_sequence_quality_scores_plot[[1]]

So each data frame is a set of x, y pairs. As it’s a histogram plot, we know that the x value is the quality score, and y is the number of times that score has been counted.

Unnesting

One possible way to process this nested data is to use tidyr:

df %>%
  tidyr::unnest(cols = plot.fastqc_per_sequence_quality_scores_plot)

As you can see, if we unnest in this way, we now have multiple rows for the same sample, which is a bit confusing (and not Tidy). However, if we use group_by and then summarise, this can be a useful way to calculate summary statistics. For example, if we want the total number of reads, we could do the following:

df %>%
  tidyr::unnest(cols = plot.fastqc_per_sequence_quality_scores_plot) %>%
  dplyr::group_by(metadata.sample_id) %>%
  dplyr::summarise(total_reads = sum(y))

Mapping with Purrr

Although unnesting worked well in this example, it can get a bit awkward for more complex operations. In these cases we can use purrr. Refer to the next section for an example that compares both approaches. The below example sums the number of reads for each sample, as we have done above, but this time it uses purrr::map_dbl to map over the list of data frames:

df %>%
  dplyr::mutate(
    total_reads = purrr::map_dbl(plot.fastqc_per_sequence_quality_scores_plot, ~sum(.$y)),
    plot.fastqc_per_sequence_quality_scores_plot = NULL
  )

Properly Handling Histogram Data

Of course, we actually want to find the median here, which is a bit harder. Luckily there exists a package called HistDat for generating summary statistics from histogram-type data. You can check out the package’s manual and vignettes here, but in brief, we want to convert each of these plot data frames into a HistDat object, which we can do using the same strategies as before. Then, using HistDat, we can calculate our summary statistics in one of the two ways mentioned above.

Using the tidyr approach, we can unnest the plot data, group it, create a HistDat object for each group, and then produce new columns using the new hist column:

df %>%
  tidyr::unnest(cols = plot.fastqc_per_sequence_quality_scores_plot) %>%
  dplyr::group_by(metadata.sample_id) %>%
  dplyr::mutate(hist = list(HistDat::HistDat(vals = x, counts = y)), .keep = "unused") %>%
  dplyr::mutate(
    mean_coverage = hist %>% dplyr::first() %>% mean(),
    median_coverage = hist %>% dplyr::first() %>% median(),
    max_coverage = hist %>% dplyr::first() %>% max(),
    hist= NULL
  ) %>%
  dplyr::slice(1)

Alternatively, using the purrr method, we can just map each plot data frame into a row of summary statistics. Much neater!

df %>%
  dplyr::mutate(
    purrr::map_dfr(plot.fastqc_per_sequence_quality_scores_plot, function(plot_df){
      hist = HistDat::HistDat(vals=plot_df$x, counts = plot_df$y)
      list(
        mean_coverage = mean(hist),
        median_coverage = median(hist),
        max_coverage = max(hist)
      )
    }),
    plot.fastqc_per_sequence_quality_scores_plot = NULL
  )

Custom Plot Parsers

So far, we have used TidyMultiqc’s built-in parsers. Each parser is a function that converts the plot JSON into a list of data frames. However, it is possible that we might encounter a new type of plot that is not yet implemented in TidyMultiqc. If this happens, the first thing you should do is file an issue against this package.

Then, if you’re daring, you can try to implement a parser. For your reference, you can refer to the existing parsers in the source code. Then, you can pass in your custom parser using the plot_parsers argument to load_multiqc:

TidyMultiqc::load_multiqc(
  multiqc_data_path, 
  sections = 'plot',
  plots = "fastqc_per_sequence_quality_scores_plot",
  plot_parsers = list(
    # This fake parser function takes a plot and just returns the iris dataset
    xy_line = function(plot_data, name){
      list(
        sample_1 = list(
          plot_name = list(iris)
        )
      )
    }
  )
)

Finally, if your parser works, please submit a pull request to TidyMultiqc to share this with everyone!

Using the strategies and patterns explained in this section, you should be in good stead to handle whatever plots MultiQC throws at you.