Longitudinal detection of circulating tumor DNA with bioinformatics analysis

Analysis of Roche KAPA Target Enrichment kit experimental data obtained on an Illumina sequencing system is most frequently performed using a variety of publicly available, open-source analysis tools.

The typical variant calling analysis workflow consists of sequencing read quality assessment, read filtering, mapping against the reference genome, duplicate removal, coverage statistic assessment, variant calling, and variant filtering.

For most of these steps, a variety of tools can be utilized. Beyond variant calling, the article provides guidance on using specific tools for data analysis in longitudinal mutation analysis applications. However, it also notes that alternative analysis workflows are applicable.

The article is designed to help readers who have a background in bioinformatics comprehend the fundamental analysis workflow that Roche employs for evaluating capture efficiency and conducting longitudinal mutation analysis. However, careful consideration of additional options is a must when deciding the most appropriate workflow for your research.

Solution

There are various free and open-source third-party tools accessible for processing raw sequencing data. These tools facilitate tasks such as converting data into suitable file formats, mapping reads to a reference sequence, assessing sequencing quality, and analyzing variant calls.

This article describes a number of steps and mini-workflows that use such third-party tools, which can be combined together into a variety of data analysis workflows.

It is recommended that users create a workflow tailored to their experimental data by using benchmark or control samples. These samples should contain somatic variants at various levels to ensure the workflow is optimized for their specific data.

Figure 1. Schematic of longitudinal detection of circulating tumor DNA (ctDNA). The plasma and germline samples are obtained at a baseline time point T₀. A new plasma sample is obtained at a follow-up time point T_N to be evaluated for ctDNA detection. The analysis modules, from FASTQ processing to variant calling, can be applied to all three samples. Afterward, the longitudinal analysis module will utilize specific output files from the 3 samples, including BAM files and variant calls, to perform an evaluation of ctDNA detection at the follow-up timepoint. Image Credit: Roche Sequencing and Life Science

Note that where the text “SAMPLE” appears throughout the examples shown here, the user should replace it with a unique sample name. Similarly, replace “DESIGN” with the name of the target enrichment design that matches the design files supplied by Roche. “NumProcessors” should be replaced with the number of CPU cores available.

Users need to substitute “/path/to/...” in the provided examples with an appropriate path on their system. It is assumed that the current directory contains all input files, and it will also serve as the location for output files and report files.

This article discusses tools that necessitate running a .jar file via Java. An exception is the Genome Analysis Toolkit (GATK), which, while requiring Java, operates through a wrapper. Users whose default system version is not Java 1.8 must run the GATK .jar file by specifying a direct path to Java 1.8.

Despite how it may appear in the printed format, the entire command for each step should be entered on a single line. Ensure that there are no spaces within a file path. However, it is important to include spaces before and after each option in the command.

The backslash “\” at the end of a line is used for line continuation, usually for a long command. There should be no spaces after the backslash.

Due to idiosyncrasies in most—if not all—PDF viewers, the underscores in command line examples in this article may not display properly at all zoom percentages. One way to confirm whether or not underscores are present is to print the page. Alternatively, users could try temporarily switching to a very high zoom percentage (e.g., 400%).

Examples included in this article show how to perform the analysis with paired end Illumina sequencing reads. Many of the tools work with single end reads if paired end reads are not available, though input and output formats may vary.

Note that using single end reads will artificially increase the duplicate rate due to the decreased ability to resolve a unique fragment from the library.

Table 1. Third-party data analysis tools were used in this white paper. The examples described in this document were tested using the software versions listed in parentheses, and different software versions may require different function calls and/or flags. See section Reference Links for installation instructions and explanations of command options. These tools were tested on a Redhat Linux system. Many GATK4 tools were originally developed as part of Picard, which maintains detailed documentation referenced throughout this white paper. Source: Roche Sequencing and Life Science

Reference genome and annotation

Before beginning the analysis, it is essential to select the appropriate human reference genome. The Longitudinal Mutation Analysis module employs the ctDNAtools package, ¹ which necessitates installing a compatible BSgenome package from Bioconductor.

The BAM alignment files inputted into ctDNAtools must align with a reference genome supported by BSgenome, like BSgenome.Hsapiens.UCSC.hg38, to ensure optimal compatibility.

The complete list of genomes that BSgenome supports can be accessed by using the “available.genomes() function” in R.

More information on BSgenome can be found at https://bioconductor.org/packages/BSgenome. Examples of available human genomes include BSgenome.Hsapiens.UCSC.hg38, BSgenome.Hsapiens.NCBI.GRCh38, BSgenome.Hsapiens.UCSC.hg19, and BSgenome.Hsapiens.1000genomes.hs37d5.

If there is a preference for using a version of the human reference genome that is not supported by BSgenome, necessary alterations must be made within the scripts of the ctDNAtools package to achieve compatibility. Alternatively, creating a custom BSgenome data package is also a viable option.

The corresponding Database of Single Nucleotide Polymorphisms (dbSNP0) annotation database (Variant Call Format) .VCF file for the human reference genome can be obtained from NCBI dbSNP Builds (https://www.ncbi.nlm.nih.gov/snp/).

The hg38 dbSNP build 138 .VCF file is also available in the GATK resource bundle (https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/ hg38/v0), which is described at https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle.

The dbSNP annotation is used to identify common Single Nucleotide Polymorphisms (SNPs) in variant calls. During reporter variant selection in the longitudinal analysis module, the common SNP variants are removed.

Index a reference genome

Before mapping can occur in most next-generation sequencing (NGS) applications, an indexed genome must be created. While various algorithms operate differently, the majority employ the Burrows-Wheeler (BWA) algorithm, which is effective for mapping millions of relatively short reads against the reference genome.

A genomic index is used to quickly identify the mapping location. Genomic indexing is required only once per genome version. The genomic index files that are created can then be used for all subsequent mapping jobs against that genome assembly version.

KDM Communications Ltd. recommends the FASTA formatted genome sequence be indexed with chromosomes in “karyotypic” sort order, i.e., chr1, chr2, ..., chr10, chr11, … chrX, chrY, chrM, etc. In these examples, reference genome files are referred to as “ref.fa”, which should be replaced by the actual file name (e.g., “hg38.fa”).

The use of an indexed reference genome in a later step is indicated by “ref.fa {indexed}” in the Input(s) section. An indexed reference genome includes the FASTA file of the genome along with all of its index files, which must be located in the same directory.

Note that multiple versions of the same reference genome can be available; thus, the user should select the appropriate genome for their analyses. It is essential for the user to confirm that the reference genome build used to create panel design files, such as primary target or capture target BED files, is identical to the one used in the analysis pipeline.

It is important to ensure that chromosome names are consistent between the two files in use. If the chosen version of the reference genome includes additional chromosomes or contigs (such as ALT contigs in the human reference), users must evaluate whether these elements are relevant and necessary for their specific analysis.

Certain areas of the human reference genome, like the pseudo-autosomal regions, can be depicted in various ways. These representations can significantly impact downstream analyses, including those related to variant analysis.

Decompress a FASTQ file

If the FASTQ files have been compressed (with a .GZ extension), some tools require them to be decompressed before use.

Examine sequence read quality

Before spending time evaluating mapping statistics, use FASTQ on raw reads and generate a per-base sequence quality plot and report to evaluate sequencing quality. The FASTQ tool can work on both compressed and uncompressed FASTQ files.

FastQC has a --threads option that allows users to specify the number of files that can be processed simultaneously. FastQC describes, “Each thread will be allocated 250MB of memory so the user shouldn’t run more threads than their available memory will cope with, and not more than 6 threads on a 32 bit machine”.

In the provided example, users have the option to use the --threads 2 command to accelerate the calculation process. For each SAMPLE input file, a .zip file is generated. Additionally, an HTML report named fastqc_report.html is created, which can be viewed using an internet browser.

The authors of FastQC have posted the following examples of the Quality Control (QC) report for a good and a bad sequencing run:

http://www.bioinformatics.babraham.ac.uk/projects/fastqc/good_sequence_short_fastqc.html

http://www.bioinformatics.babraham.ac.uk/projects/fastqc/bad_sequence_fastqc.html

Remove duplicates by utilizing unique molecular identifiers (UMIs)

The recommended methodology to obtain consensus reads from Polymerase Chain Reaction (PCR) and optical duplicates are described in this section. It is specific to sequencing data from libraries constructed using the KAPA Universal Unique Molecular Identifiers (UMIs) Adapter and the KAPA Unique Device Identification (UDI) Primer Mixes.

The methodology illustrated below consists of three steps:

1) Extraction of the UMI from insert reads;

2) Grouping of these UMIs into families/groups based on alignment coordinates and UMI sequence composition; and

3) Consensus calling of all reads within a particular UMI group.

The Roche UMI adapters described above utilize a mix of “fixed” sequence adapters. In case the exact UMI sequences are needed for additional analysis, refer to Table 2 on the following page.

Figure 2. Recommended pipeline for UMI extraction, grouping, and consensus read calling. Image Credit: Roche Sequencing and Life Science

Table 2. Sequence mix of the Roche UMI adapters (KAPA Universal UMI Adapter) The bolded T is the 3' T-overhang of the adapter. The sequence in italics is not part of the UMI sequence; it is added to increase the sequence diversity in order to ensure optimal sequencing performance. Source: Roche Sequencing and Life Science

Convert FASTQ to binary alignment/map (BAM)

The first step is to convert the demultiplexed, raw-sequencing FASTQ files to BAM files using the FastqToSam tool in GATK.

Extract UMIs from BAM

The read structure is specified as 3M3S+T. This involves extracting the first three bases, which are then stored as the UMI in the RX tag of the BAM file, as indicated by the 3M component of the structure.

Remove the next three bases from the beginning of the read (indicated as 3S). These bases form a punctuation sequence, which is used to enhance the sequence diversity and thereby ensure optimal sequencing performance.

Maintain the remaining sequence as part of the insert read (+T). The UMIs extracted from read 1 and read 2 are stored in the RX tag of the unmapped BAM file as UMI1-UMI2 (hereafter referred to as “the UMI” and considered as a single sequence).

Perform adapter trimming and quality filtering

The BAM file containing reads with extracted UMIs must be converted to a FASTQ format to facilitate adapter trimming and quality filtering. It is important to perform adapter trimming and quality filtering only after UMI extraction. This approach prevents potential bias and guarantees that only the template or insert is trimmed.

In this workflow, the unmapped BAM file is initially converted to a FASTQ file using GATK. Following this, adapter trimming and quality filtering are executed with fastp. The parameters are configured to allow automatic detection of adapter sequences, or users can manually set adapter sequences, as detailed in the Illumina Adapter Sequences Document #1000000002694 v.11 or later versions.

Note: Converting BAM files to FASTQ format does not preserve the extracted UMI information. Therefore, it is crucial to keep the output file from the UMI extraction process, named SAMPLE_unmapped_umi_extracted.bam. This file retains the UMI data stored in the RX tag, which is essential for subsequent processes after genomic alignment.

The -3 and -W 5 options allow trimming from the 3’ tail in a sliding window of 5 bp. If the mean quality is below the quality set by -q, the bases are trimmed. In addition, -u specifies what percent of bases are allowed to be unqualified before a read is discarded.

The -x and -g options turn on poly X and poly G tail trimming, respectively. The -l option means the length of the trimmed read must be at least 50 bp. The -c option turns on base correction for read pairs where read1 and read2 overlap.

The fastp application will produce two files. The SAMPLE_umi_extracted_trimmed_R1.fastq and SAMPLE_umi_extracted_ trimmed_R2.fastq contain the reads that are still paired after adapter trimming and quality filtering.

Unpaired reads can optionally be assigned to output files using the --unpaired1 and --unpaired2 options. If the user wants to increase the percentage of passing reads, the quality and length filter thresholds can be lowered.

Select a subsample of reads from a FASTQ file

Random subsampling is useful for normalizing the number of reads per set when making comparisons. With paired end reads, it is important that the two files use the same values for the seed (-s) and number of reads. The seqtk application can write the sampled reads to uncompressed FASTQ files.

The commands above will randomly subsample 40 million matched read pairs from the paired end FASTQ files for a total of 80 million reads per sample. Supplying the same random seed value (-s) ensures that the FASTQ records will remain in synchronized sort order and can be used for mapping, etc.

Note that seqtk requires an amount of Random Access Memory (RAM) proportional to the number of reads being subsampled. As the user increases the size of the subsampled read set, more RAM is needed.

Map reads to the reference genome

Adapter trimmed and quality-filtered reads are mapped to the indexed reference genome using BWA.

In the “Map Reads” step, the -R option defines the read group (“@RG”), which will appear in the BAM header. Within this string is the sample ID (“ID”), description field (“DS”), sequencing platform (“PL”), library name (“LB”), and sample name (“SM”). When a library name, ID, and sample name do not separately exist, a SAMPLE descriptor may be used, as shown in the example above.

The -Y option uses soft clipping for supplementary alignments. It is advisable to keep only reads that are aligned in proper pairs in BAM. The tool samtools_view and the flag -f2 can be used for this purpose.

Add UMI information to the reads in BAM

As UMI information is not retained during BAM to FASTQ conversion, it is necessary to merge the two BAM files containing the UMI information (SAMPLE_unmapped_umi_extracted.bam) and the alignment coordinate information (SAMPLE_ umi_aligned.bam).

The UMI information is now stored in the RX tag of the new umi_extracted_aligned_merged.bam file. MergeBamAlignment requires queryname sorted inputs.

Identify and group reads originating from the same source molecule

The GroupReadsByUmi tool in fgbio utilizes the UMI (UMI1-UMI2) and the genomic alignment start site to assign unique source molecules to each applicable read. GroupReadsByUmi implements the adjacency strategy introduced by UMI- tools.

The user can control how many errors/mismatches are allowed in the UMI sequence when assigning source molecules (--edits=n). UMI group statistics are output to a SAMPLE_umi_group_data.tsv file using the -f flag.

Note: The parameter --edits=1 will account for a single mismatch in the entire UMI sequence (UMI1+UMI2). Altering this parameter to >1 will have a significant impact on the outcome of the UMI grouping algorithm and the resultant UMI groups.

Calculate consensus sequence

The CallMolecularConsensusReads tool in fgbio processes each group of reads identified as originating from the same unique source molecule. The consensus of a group of reads can only be calculated if there are a minimum of two reads in a group.

Reads that occur as singletons are discarded by default, but this can be changed by setting the –min-reads flag to 1, and the single read will be considered the consensus.

Note: Bases with a sequencing quality less than 20 will not be used in the consensus calculation, but this can also be altered with the –min-input-base-quality flag.

Note: Here, reads are defined as those grouped into a UMI family/group, i.e., reads that have the same UMI tag and the same 5’ start position. Furthermore, the two reads minimum could consist of two reads from the same unique source molecule, or one read derived from the forward template of the unique source molecule and one derived from the reverse template of the unique source molecule.

Filter consensus reads

The FilterConsensusReads tool in fgbio is used to filter consensus reads generated by CallMolecularConsensusReads.

The option ––max-read-error-rate is the maximum raw-read error rate across the entire consensus read. The option ––max–base–error–rate is the maximum error rate for a single consensus base. The option ––max–no–call–fraction is the maximum fraction of no-calls in the read after filtering.

 The option ––min– base–quality is the minimum mean base quality across the consensus read. The option ––min–reads is the minimum number of reads supporting a consensus read. The option ––reverse–per–base–tags controls whether per-base tags should be reversed before being used on reads marked as being mapped to the negative strand.

Convert BAM to FASTQ

After consensus calling, the collapsing of the UMI groups results in the loss of alignment coordinate information. To rectify this, the SAMPLE_umi_consensus_unmapped.bam is converted to FASTQ format using SamToFastq in GATK.

Note: Loss of alignment coordinates is an inherent limitation of consensus calling and is related to alignment quality. When base information is statistically extrapolated from two or more molecules, the alignment quality is also statistically averaged. Many downstream variant callers rely on alignment quality; thus, to avoid error, the consensus reads are realigned to ensure correct alignment qualities.

Map consensus reads to the reference genome

A new SAMPLE_umi_consensus_mapped.bam file is generated after aligning the consensus reads to the indexed reference genome using BWA.

Add UMI information to the consensus reads in BAM

The final step is to merge the SAMPLE_umi_consensus_mapped.bam with the SAMPLE_umi_consensus_unmapped.bam to retain the UMI group information. This will yield an aligned BAM file with consensus reads and the UMI information retained in the RX flag. The BAM file will be coordinated and sorted, and a .bai index file will be produced.

Clip overhangs

Read overhangs can be clipped using the fgbio ClipBam tool to avoid double counting variant-supporting reads from the same template. This is recommended prior to variant calling.

Sort BAM and create index

The BAM files need to be sorted and indexed for use in the subsequent steps. Here, the indices for the BAM files are sorted and created both before and after the UMI deduplication. The BAM files can now be used for all downstream applications and analyses described by the particular NGS analysis workflow.

Detect single nucleotide variant (SNV) in baseline sample

After reads are mapped, and duplicates are removed, variants are often called against the reference genome. Somatic variants must be called by callers that are capable of detecting low abundance variants. This next section will outline how to use VarDict to call somatic variants.

VarDict is a widely recognized somatic caller, renowned for its ultra-sensitivity in detecting variants from targeted sequencing data. VarDict’s performance scales linearly to sequencing depth, and it enables ultra-deep sequencing for tumor evolution and liquid biopsy.

Originally written in Perl as vardict.pl, the tool has been developed with a Java-based replacement, which yields ten times more acceleration than the Perl implementation. For longitudinal analysis, the set of SNVs detected in the baseline sample will be used as reporter variant candidates for analysis.

Annotate variants

The variants can be annotated with dbSNP using the SnpSift tool. The dbSNP annotation is used in the longitudinal analysis module to remove common SNPs.

VCF to table

The annotated VCF file can be parsed into a tab-delimited table for downstream processing using GATK VariantsToTable. VCF standard and INFO fields can be extracted using the option -F. The genotype FORMAT fields can be extracted using the option -GF.

The output text file is used for downstream longitudinal analysis. The recommended fields to extract from the vardict VCF file should include the following: AF: Allele frequency; DP: Depth – total depth; VD: AltDepth – variant depth; MQ: mapping quality; QUAL: average base quality at a variant position; SBF: Strand Bias Fisher p-value; NM: mean mismatches in reads.

These fields are used for filtering during the reporter variant selection step in longitudinal mutation analysis.

Longitudinal mutation analysis

Description

Longitudinal mutation analysis involves using a set of reporter variants detected in a baseline cell-free DNA (cfDNA) sample to track and detect the same variants in subsequent follow-up cfDNA sample(s).

Detection in the follow-up sample is based on the number of variant-supporting reads observed. A mutation-positive follow-up sample would have variant-supporting reads that are significantly higher than the background error rate.

An empirical p-value from a Monte-Carlo sampling-based test is used for mutation positivity calling, based on the approach by Newman et al.² The ctDNAtools package provides tools to implement the longitudinal mutation analysis. The main components included in the analysis are 1) Generate background panel and blocklist, 2) Select reporter variants, and 3) Evaluate longitudinal mutations.

Three different samples are used for the longitudinal analysis. The reporter variants are obtained from the baseline cfDNA sample. The germline sample is used to remove reporters observed in the germline.

The follow-up sample is tested for mutation positivity based on the reporter variants. If multiple follow-up samples are subsequently obtained, the analysis can be performed using the same respective baseline and germline samples.

In addition, a set of normal cfDNA samples is required to generate a longitudinal analysis blocklist, which is used to establish a list of error-prone mutations in the panel target regions.

A list of R packages required by ctDNAtools can be found at https://github.com/alkodsi/ctDNAtools/blob/master/DESCRIPTION.

Importantly, the appropriate BSgenome.Hsapiens library needs to be used. (See previous section “Reference Genome and Annotation”).

In addition, the R package, optparse, is required for implementing the modules described in this article using custom scripts. All of the custom scripts described here in the longitudinal analysis can be obtained from https://github. com/Roche-CSI/kapa_nhl. This table summarizes the R packages required:

Generate background panel and blocklist

Longitudinal analysis involves the detection of minor amounts of variant-supporting reads, often below 0.1% allele frequency (AF). A blocklist is highly recommended to filter out variants that are likely to be false positives.

This section will discuss how to generate a blocklist using a set of normal samples with the ctDNAtools package in custom scripts.

The blocklist only needs to be generated once for the longitudinal workflow setup. After the blocklist is obtained, all subsequent analyses will be able to use the same blocklist.txt file. A new blocklist is recommended whenever there are major changes to the workflow, such as new reagents, capture probe panels, or protocols.

The create_bg_panel.R and the create_blocklist.R scripts will utilize the functions create_background_panel and create_black_list from ctDNAtools, respectively.

A background panel file is required prior to generating the blocklist. This is obtained using the create_bg_panel.R script. A .csv file without column headers that contain the path to the bam files from normal samples is used as input, as well as the target region bed file.

The output is a panel_background.RDS file that contains information such as read depth and alt AF for background mutations. The file can be read into the R environment using the function readRDS().

Two types of blocklists are available in the ctDNAtools package. The first type is a list of genomic loci (chr_pos, regardless of substitutions). The second type is substitution-specific (chr_pos_ref_alt), which offers a more specific variant-level list for identifying errors.

The blocklist type can be set using the option ––blist_type during the create_bg_panel.R process. Setting it to “loci” will implement a loci-based list, whereas setting it to “variant” will implement a substitution-specific list.

Lastly, the BSgenome.Hsapiens version should match the reference for the input sample bam files. For example, setting the option ––reference to “BSgenome.Hsapiens.UCSC.hg38” will load BSgenome.Hsapiens.UCSC.hg38.

After the panel_background.RDS file is obtained, the create_blocklist.R script can be used to generate the blocklist.txt file. The selection of parameters in the blocklist creation is critical to identify which loci are considered to be background noise.

The quantile of mean VAF above which the loci are considered noise is set using ––blocklist_ vaf_quantile. A blocklist should be built using a sufficient number of normal samples to capture the error-prone regions in the panel.

The minimum number of normal samples that exhibit at least one non-reference read to be considered noisy is set using ––blocklist_min_samples_one_read. The minimum number of normal samples that exhibit at least two non-reference reads to be considered noisy is set using ––blocklist_min_samples_two_reads.

The minimum number of normal samples that exhibit at least n non-reference reads is set using ––blocklist_min_samples_n_reads. It is recommended to set ––blocklist_min_samples_one_read at greater than 60% of the number of normal samples.

For example, if there are 20 normal samples used for generating the blocklist, ––blocklist_min_samples_one_read should be set to 12. This would allow error-prone variants to be captured efficiently when they occur in at least 12 normal samples.

Similarly, for ––blocklist_min_samples_two_reads, it is recommended to be set at greater than 50% of the number of normal samples. Parameters for blocklist generation should be experimentally determined to obtain optimal results.

Select reporter variants

The list of reporter variant candidates after variant calling by vardict-java should be filtered to obtain confident reporters for longitudinal analysis. Filtering according to variant metrics is performed, followed by germline presence detection.

Only SNVs are used for longitudinal analysis. Reporter variant candidates that do not meet filtering criteria will be removed from downstream longitudinal analysis. The recommended filters for keeping high-confidence somatic SNV reporters are: FILTER=PASS, AF>0.5%, AF<35%, DP>1000, VD>15, MQ>55, QUAL>45.

In addition, variants included in the dbSNP database are removed. More stringent filtering may be applied to reduce the number of false reporter variants. These options are available in the filter_reporters.R script. The default filtering thresholds for the script are shown on the following page.

After variant filtering, the remaining reporter variant candidates are examined for germline sample presence. The germline sample bam file is set using the --germline option. If the same variant is detected in the germline sample, it is removed from downstream analysis.

The germline presence AF threshold can be set using the --germline_af option. A stringent cutoff, such as 0.0005, will remove reporter variant candidates that have detectable reads in the germline sample at 0.05% AF. Reporter candidates included in the blocklist.txt file will also be removed from downstream longitudinal analysis.

The final selected list of baseline reporter variants used for longitudinal analysis are saved in the output file selected_baseline_reporters.txt. To assist in the manual review process, all potential reporters will be annotated with corresponding germline sample data and follow-up sample read counts. This information will be compiled in the ‘all_reporters.txt’ file after the variant filtering step.

The reads to be included for counting can be filtered by ––read_min_bq (minimum base quality) and ––read_ min_mq (minimum mapping quality). The ––read_max_dp option is the maximum depth above which sampling will happen.

Evaluate longitudinal mutations

Finally, the set of selected reporters are used for evaluating longitudinal mutations and determining mutation positivity in the follow-up cfDNA sample. The longitudinal_analysis.R script will run the test_ctDNA function from ctDNAtools.

The reporters listed in the selected_baseline_reporters.txt file will be used to query the follow-up sample BAM file FOLLOWUP_SAMPLE_umi_deduped.clipov.sorted.bam to evaluate mutation positivity.

The primary steps are as follows:

1. Quantifying the read counts for both reference and variant alleles of the reporter variants.
2. Estimating the background rate of the tested sample.
3. Utilizing a Monte Carlo-based sampling test to calculate an empirical p-value.

The p-value will only be used to determine positivity if the number of informative reads (number of all unique reads covering the mutations) exceeds the specified threshold. Otherwise, the sample will be considered undetermined.

The threshold for the number of informative reads is established through the ––reads_threshold option. The p-value threshold for determining ctDNA positivity or negativity is set using the ––pvalue_threshold option. To optimize sensitivity and specificity, it is important to carefully determine the appropriate p-value threshold.

A set of normal cfDNA samples can be used to test the range of p-values to call ctDNA negativity. When calculating the background rate, the ––vaf_threshold is used to ignore the bases with higher than the VAF threshold (likely true mutations).

The ––read_min_bq, ––read_min_mq, and ––read_max_dp are options used for read counting, which are the same options as the select_reporters.R script in the previous step. The ––n_sim option is the number of Monte Carlo simulations.

Basic mapping metrics

Basic mapping metrics can be calculated using GATK CollectAlignmentSummaryMetrics.

Calculate base mismatch error rate

The custom script mismatch_rate.R (available at https://github.com/Roche-CSI/kapa_nhl) is used to obtain the base mismatch rate from a BAM file. The script will run the function get_background_rate from ctDNAtools.

If the blocklist.txt file is used, the mismatch rate will be computed on the target regions after excluding the loci in the blocklist. The output result file SAMPLE_mismatch_rate.csv will include the general mismatch rate and substitution-specific rates.

Count optical duplicates

GATK MarkDuplicates is used to count the number of optical duplicates.

Estimate insert size distribution

The DNA used in sequence capture undergoes a process of random fragmentation, followed by a subsequent size selection stage.

It is normal to observe a range of fragment sizes, but if skewed too large or too small, the on-target rate and/or percent of bases covered with at least one read can be adversely affected. The size of these fragments can be estimated from paired end sequencing reads (will not work for single end reads).

Count on-target reads

Use GATK CountReads to calculate the number of reads that overlap a target BED file by at least 1 bp. Calculation of on-target reads is one measure of the success of a Roche TE experiment, though optimal on-target is design-specific. The on-target metric is affected by library insert size, hybridization and wash stringency, and laboratory protocol.

It is necessary to apply filters to include specific reads for analysis. In the example commands above, “MappedReadFilter” and “NotSecondaryAlignmentReadFilter” are used to filter out reads that are unmapped or represent secondary alignments.

See https://gatk.broadinstitute.org/hc/en-us/articles/360057438571--Tool-Documentation-Index#ReadFilters for a full list of available read filters.

Divide “the number of on-target reads” found in ontarget_reads_umi_deduped.txt by “the total number of mapped, non-duplicate reads” to get “the percentage of on-target reads after duplicates removal”.

Similarly, divide “the number of on-target reads” found in ontarget_reads_sorted.txt by “the total number of mapped reads” to get “the percentage of on-target reads before duplicates removal”. See Basic Mapping Metrics for reporting of the total number of mapped and non-duplicate reads.

Target-adjacent coverage is typical for target enrichment due to the capture of partially on-target DNA library fragments that also extend outside the capture region.

To optionally assess the number of reads that are target adjacent, add --interval_ padding 100 to the commands above to add 100 bp to both sides of all targets.

Although 100 bp is commonly used for this kind of padding, shorter or longer lengths may also be appropriate depending on expected library fragment sizes. Please note: all remaining steps in this article use non-padded targets.

Create genomic interval lists

Interval lists are genomic interval description files required by GATK CollectHsMetrics that contain a SAM-like header describing the reference genome and a set of coordinates with strand and name for each interval.

The Roche-provided “primary target” files can be provided as GATK “target interval” inputs, and the Roche-provided “capture target” files can be provided as GATK “bait interval” inputs.

However, here, the focus is on assessing the performance of sequence capture on the ‘capture target.’ Specifically, the aim is to apply the --PER_BASE_COVERAGE and --PER_TARGET_COVERAGE options in GATK’s CollectHsMetrics to the ‘capture target.’ This will provide detailed coverage data for each base and each target region.

An interval list for the “capture target” file will be provided as both the “target interval” input and the “bait interval” input in GATK CollectHsMetrics.

Use the GATK BedToIntervalList command to create Interval List files from target BED files.

Create padded target bed files

A padded target bed file can be useful when the analysis includes on-target DNA fragments that extend slightly beyond the capture regions. For example, variant calling can be expanded to those adjacent target regions by using a padded target bed file.

bedtools slop will be used to add padding to the target regions. The number of base pairs to add is set using the option -b. For example, setting -b to 100 will add 100 bp to both directions of a target region. bedtools merge combines overlapping or “book-ended” target regions into one region.

Hybrid selection (HS) analysis metrics

The CollectHsMetrics command calculates a number of metrics assessing the quality of target enrichment reads.

Here, the same interval file is supplied to both --TARGET_INTERVALS and --BAIT_INTERVALS parameters of GATK CollectHsMetrics, as the overall aim is to leverage the --PER_BASE_COVERAGE and --PER_TARGET_COVERAGE options to examine the per base and per target coverages in those capture regions.

See https://broadinstitute.github.io/picard/picard- metric-definitions.html#HsMetrics for a description of output metrics.

It is important to mention that some metrics are not directly comparable as some are obtained before read or base filters are applied (e.g., capture bases metrics) while others are calculated after (e.g., target coverage metrics).

Note: The CollectHsMetrics tool reports the percent of bases covered at certain sequencing depths (e.g., 1X, 10X, 20X, 30X, 40X, and 50X). To obtain coverage for additional sequencing depths ≥N use the command “gawk ‘$1 != “chrom” && $4 >= N’ SAMPLE_per_base_coverage_umi_deduped.txt | wc -l” to obtain the number of bases with at least N coverage.

Divide this number by the “TARGET_TERRITORY” value from SAMPLE_hs_metrics_umi_deduped.txt to calculate “% bases ≥ N”. The “SAMPLE_per_base_coverage*.txt” files can be quite large for large designs and can be compressed once sequencing depths have been calculated using gzip as described earlier.

Calculate coverage in exonic target regions

The percentage of exonic positions with coverage higher than a cutoff is an important metric to evaluate the quality of target enrichment in exonic target regions. It can be calculated by leveraging the “SAMPLE_per_base_coverage_umi_deduped.txt” file generated by GATK CollectHsMetrics.

Note: For coverage at sequencing depths of N or higher, the command gawk ‘$4 >= N’ SAMPLE_per_base_coverage_umi_deduped_exon.bed | wc -l can be used. This command calculates the number of bases in exonic regions with a minimum of N coverage. To determine the “% Panel exon region ≥N,” divide this number by the “TARGET_TERRITORY” value found in the SAMPLE_hs_metrics_umi_deduped.txt file.

Description of metrics

The tools used in this article generate output files that contain many metrics. There are some metrics that are frequently monitored to assess capture experiment performance. The table below describes various metrics, the tool(s) used for generating these metrics, and any additional mathematical or string parsing operations that might be required to derive the final values.

Table 3. Description of important metrics for DNA samples. Source: Roche Sequencing and Life Science

References

1. Alkodsi A, Meriranta L, Pasanen A, Leppä S. ctDNAtools: An R package to work with sequencing data of circulating tumor DNA. bioRxiv (2020).
2. Newman, A., Lovejoy, A., Klass, D. et al. Integrated digital error suppression for improved detection of circulating tumor DNA. Nat Biotechnol 34, 547 – 555 (2016).

About Roche Sequencing and Life Science

Roche Sequencing & Life Science is part of Roche Diagnostics, which, along with Roche Pharmaceuticals, plays an important role in modern healthcare. Roche Diagnostics’ broad range of innovative diagnostic tests and systems play a pivotal role in the groundbreaking area of integrated healthcare solutions and cover the early detection, targeted screening, evaluation and monitoring of disease. Roche Diagnostics is active in all market segments, from scientific research and clinical laboratory systems to patient self-monitoring.

Sponsored Content Policy: News-Medical.net publishes articles and related content that may be derived from sources where we have existing commercial relationships, provided such content adds value to the core editorial ethos of News-Medical.Net which is to educate and inform site visitors interested in medical research, science, medical devices and treatments.