Workflows that generate, aggregate, and visualize quality control files for MSK-ACCESS.

• Free software: Apache Software License 2.0

• Documentation: https://msk-access.gitbook.io/access-quality-control-v2/

Installation

Clone the repository:

git clone --depth 50 https://github.com/msk-access/access_qc_generation

Credits

• CMO cfDNA Informatics Team

• Cookiecutter: https://github.com/audreyr/cookiecutter

• audreyr/cookiecutter-pypackage: https://github.com/audreyr/cookiecutter-pypackage

Given the output files from Nucleo, there are workflows to generate the quality control files, aggregate them files across many samples, and visualize them using MultQC. You can choose to run these workflows whether you have just one or hundreds of samples. Depending on your use case, there are two main options:

(1) Run qc_generator.cwl followed by aggregate_visualize.cwl. This approach first generates the QC files for one or more samples, and you use the second CWL script to aggregate the QC files and visualize them with MultiQC. This option can be useful for when you want to generate the QC files for some samples just once and then reuse those samples in multiple MultiQC reports.

(2) Run just access_qc.cwl. This option just combines the two steps from the first option into one workflow. This workflow can

Step 1: Create a virtual environment.

Option (A) - if using cwltool

If you are using cwltool only, please proceed using python 3.6 as done below:

Option (B) - recommended for Juno HPC cluster

If you are using toil, python 3 is required. Please install using Python 3.6 as done below:

Step 2: Clone the repository

Step 3: Install requirements using pip

We have already specified the version of cwltool and other packages in the requirements.txt file. Please use this to install.

Step 4: Generate an inputs file

For details on how to create this file, please follow this example (there is a minimal example of what needs to be filled in at the end of the page):

It's also possible to create and fill in a "template" inputs file using this command:

Once we have successfully installed the requirements we can now run the workflow using cwltool/toil .

Step 5: Run the workflow

Requirements

Following are the requirements for running the workflow:

• • Python Packages (will be installed as part of pipeline installation):

    • toil[cwl]==5.1.0

    • pytz==2021.1

    • typing==3.7.4.3

    • ruamel.yaml==0.16.5

    • pip==20.2.3

    • bumpversion==0.6.0

    • wheel==0.35.1

    • watchdog==0.10.3

    • flake8==3.8.4

    • tox==3.20.0

    • coverage==5.3

    • twine==3.2.0

    • pytest==6.1.1

    • pytest-runner==5.2

    • coloredlogs==10.0

    • pytest-travis-fold==1.3.0

Parameter Used by Tools

Common Parameters Across Tools

Template Inputs File

biometrics_bed_file:
  class: File
  path: /path/to/MSK-ACCESS-v1_0-probe-B.sorted.bed
biometrics_json: true
biometrics_plot: true
biometrics_vcf_file:
  class: File
  path: /path/to/MSK-ACCESS-v1_0-TilingaAndFpSNPs.vcf
collapsed_bam:
- class: File
  path: /path/to/bam
collapsed_biometrics_coverage_threshold: null
collapsed_biometrics_major_threshold: null
collapsed_biometrics_min_base_quality: null
collapsed_biometrics_min_coverage: null
collapsed_biometrics_min_homozygous_thresh: null
collapsed_biometrics_min_mapping_quality: null
collapsed_biometrics_minor_threshold: null
duplex_bam:
- class: File
  path: /path/to/bam
duplex_biometrics_major_threshold: null
duplex_biometrics_min_base_quality: null
duplex_biometrics_min_coverage: null
duplex_biometrics_min_homozygous_thresh: null
duplex_biometrics_min_mapping_quality: null
duplex_biometrics_minor_threshold: null
group_reads_by_umi_bam:
- class: File
  path: /path/to/bam
hsmetrics_coverage_cap: 30000
hsmetrics_minimum_base_quality: 1
hsmetrics_minimum_mapping_quality: 1
noise_sites_bed:
  class: File
  path: /path/to/MSK-ACCESS-v1_0-probe-A_no_msi_sorted_deduped.bed
pool_a_bait_intervals:
  class: File
  path: /path/to/MSK-ACCESS-v1_0-probe-A_baits.sorted.interval_list
pool_a_target_intervals:
  class: File
  path: /path/to/MSK-ACCESS-v1_0_panelA_targets.interval_list
pool_b_bait_intervals:
  class: File
  path: /path/to/MSK-ACCESS-v1_0-probe-B_baits.sorted.interval_list
pool_b_target_intervals:
  class: File
  path: /path
reference:
  class: File
  path: /path
sample_group:
- patient_id
sample_name:
- sample_id
sample_sex:
- M
sample_type:
- tumor
sequence_qc_min_basq: 1
sequence_qc_min_mapq: 1
sequence_qc_threshold: null
sequence_qc_truncate: null
simplex_bam:
- class: File
  path: /path
uncollapsed_bam:
- class: File
  path: /path
uncollapsed_bam_base_recal:
- class: File
  path: /path

Introduction

At the top of the MultiQC report are one or two tables showing some per-sample information. One table is for plasma samples and another is for buffy-coat samples; so only one table may show up depending on your sample composition.

Interpretation

In the above figure you'll notice that most columns are highlighted as either red, yellow or green, which indicates if the metric fails, is borderline, or passes the thresholds set for each, respectively. This allows you to quickly glance at all samples to see where potential issues are. Below are the descriptions for each column and were the data was obtained from.

Introduction

This figure plots the normalized coverage against the % GC content from the ACCESS target regions. Each line is data from one sample.

Methods

Tool used: GATK-CollectHsMetrics BAM type: (1) collapsed BAM and (2) uncollapsed BAM. Regions: Pool A

The data used to produce this figure are the values under the normalized_coverage and %gc columns, which are in the *_per_target_coverage.txt output file from CollectHsMetrics. For each sample separately, the % GC content for each target region is calculated, followed by binning the target regions by their GC content (in 5% intervals). Then for each bin, the mean coverage is calculated and then normalized across all regions that fall into each GC bin.

Interpretation

Extreme base compositions, i.e., GC-poor or GC-rich sequences, lead to an uneven coverage or even no coverage of reads across the genome. This can affect downstream small variant and copy number calling. Both of which rely on consistent sequencing depth across all regions. Ideally this plot should be as flat as possible. The above example depicts a slight decrease in coverage at really high GC-rich regions, but is a good result for ACCESS.

Introduction

This figure shows the insert size distribution from the ACCESS target regions. Insert size is calculated from the start and stop positions of the reads after mapping to the reference genome.

Methods

Tool used: GATK-CollectInsertSizeMetrics BAM type: Collapsed BAM Regions: Pool A

The data used to produce this figure are the values under the MODE_INSERT_SIZE column contained in the output file from CollectInsertSizeMetrics.

Interpretation

Cell free DNA has distinctive features due to the natural processes behind its fragmentation. One such feature is the set of 10-11 bp fluctuations that indicate the preferential splicing of fragments due to the number of bases per turn of the DNA helix, which causes a unique pattern of binding to the surface of histones.

The more pronounced peak at 166 bp indicate complete wrapping of the DNA around the histones' circumference, and similarly the second more pronounced peak indicates two complete wraps.

Buffy coat samples are mechanically sheared and thus do not exhibit these distinctive features, hence the different shape for their distribution.

Introduction

This figure shows the density plot of coverage values from the ACCESS target regions. Each line is data from one sample. Each sample is normalized by the median coverage value of that sample to align all peaks with one another and correct for sample-level differences.

Methods

Tool used: GATK-CollectHsMetrics BAM type: Collapsed BAM Regions: Pool A

The data used to produce this figure are the values under the normalized_coverage column, which are in the *_per_target_coverage.txt output file from CollectHsMetrics. Then the gaussian_kde function from the python scipy package is used to produce the density plot.

Interpretation

Each distribution should be unimodal, apart from a second peak on the low end due to X chromosome mapping from male samples. Narrow peaks are indicative of evenly distributed coverage across all bait regions. Wider distributions indicate uneven read distribution, and may be correlated with a large GC bias. Note that the provided bed file lists start and stop coordinates of ACCESS design probes, not the actual genomic target regions.

Introduction

There are several sections displaying bait set capture efficiency. Each section corresponds to a separate BAM type and bait set combination. The tool used to produce the metrics is GATK-CollectHsMetrics. By default, only the mean bait coverage, mean target coverage, and % Usable bases on-target are displayed. However, there are many more metrics that can be toggled to display by clicking on the Configure Columns button.

Methods

Tool used: GATK-CollectHsMetrics BAM type: (1) Uncollapsed BAM, (1) Collapsed BAM, (1) Duplex BAM, and (4) Simplex BAM Regions: Pool A and Pool B

Interpretation

The aim is to have high coverage across Pool A and Pool B panels.

Introduction

This figure shows the mean base quality by cycle for before and after BaseQualityScoreRecalibration (BQSR). The sequencer uses the difference in intensity of the fluorescence of the bases to give an estimate of the quality of the base that has been read. The BQSR tool from GATK recalculates these values based on the empirical error rate of the reads themselves, which is a more accurate estimate of the original quality of the read.

Methods

Tool used: GATK-MeanQualityByCycle BAM type: Uncollapsed BAM. Regions: Pool A

Interpretation

It is normal to see a downwards trend in pre and post-recalibration base quality towards the ends of the reads. Average post-recalibration quality scores should be above 20. Spikes in quality may be indicative of a sequencer artifact.

Introduction

Two metrics are used to estimate sample contamination: minor contamination and major contamination. Moreover, minor contamination is calculated separately for collapsed and duplex BAMs. Both contamination metrics are produced by the fingerprinting SNP set. However, minor contamination is calculated using just the homozygous sites, whereas the major contamination is via the ratio of heterozygous to homozygous sites. For each contamination-BAM type combination there is a table showing per-sample contamination values and any associated metrics.

Methods

It is a two step process to produce the table: (1) extract SNP genotypes from each sample using biometrics extract command and (2) perform a pairwise comparison of all samples to determine sample relatedness using the biometrics minor and biometrics major commands. Please see the biometrics documentation for further documentation on the methods.

Interpretation

Minor contamination

Samples with minor contamination rates of >0.002 are considered contamination.

Major contamination

The fraction of heterozygous positions should be around 0.5. If the fraction is greater than 0.6, it is considered to have major contamination.

Introduction

This section contains a table showing the samples clustered into groups, where each row in the table corresponds to one sample. The table will show whether your samples are grouping together in unexpected ways, which would indicate sample mislabelling.

Methods

Tool used: biometrics BAM type: Collapsed BAM Regions: MSK-ACCESS-v1_0-curatedSNPs.vcf

It is a two step process to produce the table: (1) extract SNP genotypes from each sample using biometrics extract command and (2) perform a pairwise comparison of all samples to determine sample relatedness using the biometrics genotype command. Please see the biometrics documentation for further documentation on the methods.

Interpretation

Below is a description of all the columns.