Inspecting the PASS/WARN/FAIL Status of each module

The simplest of the plots is to summarise the PASS/WARN/FAIL flags as produced by FastQC for each module. This plot can be simply generated using plotSummary()

plotSummary(fdl)

Default summary of FastQC flags.

Visualising Read Totals

The next most informative plot may be to summarise the total numbers of reads in each associated Fastq file. By default, the number of duplicated sequences from the Total_Duplicated_Percentage module are shown, but this can be disabled by setting duplicated = FALSE.

plotReadTotals(fdl)

As these are ggplot2 objects, the output can be modified easily using conventional ggplot2 syntax. Here we’ll move the legend to the top right as an example.

plotReadTotals(fdl) +
    theme(
        legend.position = c(1, 1), 
        legend.justification = c(1, 1),
        legend.background = element_rect(colour = "black")
    )

Per Base Sequence Qualities

Turning to the Per base sequence quality scores is the next most common step for most researchers, and these can be obtained for an individual file by selecting this as an element (i.e. FastqcData object ) of the main FastqcDataList object. This plot replicates the default plots from a FastQC report.

plotBaseQuals(fdl[[1]])

Example showing the Per_base_sequence_quality plot for a single FastqcData object.

When working with multiple FastQC reports, these are summarised as a heatmap using the mean quality score at each position.

plotBaseQuals(fdl)

Example showing the Mean Per Base Squence Qualities for a set of FastQC reports.

Boxplots of any combinations can also be drawn from a FastqcDataList by setting the argument plotType = "boxplot". However, this may be not suitable for datasets with a large number of libraries.

plotBaseQuals(fdl[1:4], plotType = "boxplot")

Mean Sequence Quality Per Read

Similarly, the Mean Sequence Quality Per Read plot can be generated to replicate plots from a FastQC report by selecting the individual file from the FastqcDataList object.

plotSeqQuals(fdl[[1]])

Example plot showing Per_sequence_quality_scores for an individual file.

A heatmap of mean sequence qualities can be generated when inspecting multiple reports.

plotSeqQuals(fdl)

Example heatmaps showing Per_sequence_quality_scores for a set of files.

An alternative view may be to plot these as overlaid lines, which can be simply done by setting plotType = "line". Again, discretion should be shown when choosing this option for many samples.

r2 <- grepl("R2", names(fdl))
plotSeqQuals(fdl[r2], plotType = "line")

Per Base Sequence Content

The Per_base_sequence_content module can also be plotted for an individual report with the layout being identical to that from FastQC.

plotSeqContent(fdl[[1]])

Individual Per_base_sequence_content plot

These are then combined across multiple files as a heatmap showing a composite colour for each position. Colours are combined using rgb() with A,C, G and T being represented by green, blue, black and red respectively.

plotSeqContent(fdl)

Combined Per_base_sequence_content plot

NB These plots can be very informative setting the argument usePlotly = TRUE, however they can be slow to render given the nature of how the package plotly renders graphics.

Again, supplying multiple files and setting plotType = "line" will replicate multiple individual plots from the original FastQC reports. The argument nc is passed to facet_wrap() from the package ggplot2 to determine the number of columns in the final plot.

plotSeqContent(fdl[1:2], plotType = "line", nc = 1)

Adapter Content

Adapter content as identified by FastQC is also able to be plotted for an individual file.

plotAdapterContent(fdl[[1]]) 

Adapter Content plot for a single FastQC report

When producing a heatmap across a set of FastQC reports, this will default to Total Adapter Content, instead of showing the individual adapter types.

plotAdapterContent(fdl)

Heatmap showing Total Adapter Content by position across a set of FastQC reports

Sequence Duplication Levels

As with all modules, the Sequence Duplication Levels plot is able to be replicated for an individual file.

plotDupLevels(fdl[[1]])

Example Sequence Duplication Levels plot for an individual file.

When plotting across multiple FastQC reports, duplication levels are shown as a heatmap based on each default bin included in the initial FastQC reports. By default, the plotted values are the “Pre” de-duplication values. Note that values are converted to percentages instead of read numbers to ensure comparability across files. In the plot below, CCGC_R2 shows very low duplication levels, whilst ATTG_R1 shows high levels of duplication. The commonly observed ‘spikes’ around >10 are also evident as the larger red blocks.

plotDupLevels(fdl)

Sequence Duplication Levels for multiple files

The class TheoreticalGC

A selection of Theoretical GC content is supplied with the package in the object gcTheoretical, which has been defined with the additional S4 class TheoreticalGC. GC content was calculated using scripts obtained from
https://github.com/mikelove/fastqcTheoreticalGC. Available genomes and transcriptomes can be obtained using the function gcAvail() on the object gcTheoretical and specifying the type.

gcAvail(gcTheoretical, "Genome")

## # A tibble: 24 × 4
##    Name          Group    Source Version   
##    <chr>         <chr>    <chr>  <chr>     
##  1 Alyrata       Plants   JGI    V1.0      
##  2 Amellifera    Animals  UCSC   apiMel2   
##  3 Athaliana     Plants   TAIR   Arapot11  
##  4 Btaurus       Animals  UCSC   bosTau8   
##  5 Celegans      Animals  UCSC   ce11      
##  6 Cfamiliaris   Animals  UCSC   canFam3   
##  7 Dmelanogaster Animals  UCSC   dm6       
##  8 Drerio        Animals  UCSC   danRer10  
##  9 Ecoli         Bacteria NCBI   2008/08/05
## 10 Gaculeatus    Other    UCSC   gasAcu1   
## # ℹ 14 more rows

Inspecting GC Content

As with all modules, data for an individual file replicates the default plot from a FastQC report, but with one key difference. This is that the Theoretical GC content has been provided in the object gcTheoretical based on 100bp reads. This empirically determined content is shown as the Theoretical GC content line.

plotGcContent(fdl[[1]], species = "Hsapiens", gcType = "Transcriptome")

Example GC Content plot using the Hsapiens Transcriptome for the theoretical distribution.

Again, data is summarised as a heatmap when plotting across multiple reports, with the default value being the difference between the observed and the theoretical GC content.

plotGcContent(fdl)

Example GC content showing the difference between observed and theoretical GC content across multiple files.

Line plots can also be produce an alternative viewpoint, with read totals displayed as percentages instead of raw counts.

plotGcContent(fdl, plotType = "line",  gcType = "Transcriptome")

Example GC content plot represented as a line plot instead of a heatmap.

Customized theoretical GC content can be generated using input DNA sequences from a supplied fasta file.

faFile <- system.file(
    "extdata", "Athaliana.TAIR10.tRNA.fasta", 
    package = "ngsReports")
plotGcContent(fdl, Fastafile = faFile, n = 1000)

Log file import

Bowtie log files can be parsed and imported

fls <- c("bowtiePE.txt", "bowtieSE.txt")
bowtieLogs <- system.file("extdata", fls, package = "ngsReports")
df <- importNgsLogs(bowtieLogs, type = "bowtie")

df %>%
    dplyr::select("Filename", starts_with("Reads")) %>%
    pander(
        split.tables = Inf,
        style = "rmarkdown", 
        big.mark = ",",
        caption = "Select columns as an example of output from bowtie."
    )

STAR log files can be parsed and imported

starLog <- system.file("extdata", "log.final.out", package = "ngsReports")
df <- importNgsLogs(starLog, type = "star")

The output of the samtools flagstat module can be parsed and imported

flagstatLog <- system.file("extdata", "flagstat.txt", package = "ngsReports")
df <- importNgsLogs(flagstatLog, type = "flagstat")

In addition to the files produced by the above alignment tools, the output from Duplication Metrics (picard) can also be imported. This is imported as a list with a tibble containing the detailed output in the list element $metrics and the histogram data included as the second elements.

sysDir <- system.file("extdata", package = "ngsReports")
fl <- list.files(sysDir, "Dedup_metrics.txt", full.names = TRUE)
dupMetrics <- importNgsLogs(fl, type = "duplicationMetrics", which = "metrics")
str(dupMetrics)

## tibble [1 × 10] (S3: tbl_df/tbl/data.frame)
##  $ Library                       : chr "~/myExpt/2_alignedData/bams/Sample1.bam"
##  $ Unpaired Reads Examined       : int 93
##  $ Read Pairs Examined           : int 43490914
##  $ Secondary Or Supplementary Rds: int 5366030
##  $ Unmapped Reads                : int 0
##  $ Unpaired Read Duplicates      : int 79
##  $ Read Pair Duplicates          : int 11494386
##  $ Read Pair Optical Duplicates  : int 642648
##  $ Percent Duplication           : num 0.264
##  $ Estimated Library Size        : int 69609522

Plot log file imports

Summaries of log files from select mapping and alignment tools can be plot using the function plotAlignemntSummary().

plotAlignmentSummary(bowtieLogs, type = "bowtie")

Example Bowtie logs for PE and SE sequencing

plotAlignmentSummary(starLog, type = "star")

Example STAR aligner logs