DNA Analysis Methods
Contamination Detection
The contamination analysis step detects foreign human DNA contamination using the SNP error file and Pileup file that are generated during the small variant calling, and the TMB trace file. The software determines whether a sample has foreign DNA using the contamination score. The contamination score is the sum of all the log likelihood scores across the pre-defined SNP positions whose minor allele frequency is <25% in the sample and not likely due to CNV events.
In contaminated samples, the variant allele frequencies in SNPs shift away from the expected values of 0%, 50%, or 100%. The algorithm collects all the positions that overlap with common SNPs that have variant allele frequencies of < 25% or > 75%. Then, the algorithm computes the likelihood that the positions are an error or a real mutation by using the following qualifications:
Estimates error rate per sample
Counts mutation support
Counts total depth
Copy Number Variant Calling
The copy number variant caller performs amplification, reference, and deletion calling for CNV targets within the assay. It counts the coverage of each target interval on the panel, performs normalization, calculates fold change values, and determines the CNV status for each CNV target. During normalization steps, coverage biases are caused by the following potential variables: sequencing depth, target size, PCR duplicates, probe efficiency, GC bias, and DNA type. A collection of normal cell-free DNA and contrived samples is used to correct for some of these variables. The outputs is VCF file. Amplification and deletion are annotated as DUP and DEL in the VCF file.
DNA Alignment and Read Collapsing
The alignment step uses DRAGEN Aligner with UMI collapsing to align DNA sequences in FASTQ files to the hg19_decoy genome. This step combines sets of reads (ie, families) that are grouped together based on genomic locations and UMI tags into representative sequences. This process accurately removes duplicate reads and sequencing errors without losing the signal of very low frequency (< 1%) sequence variations.
This alignment step generates BAM files (*.bam) and BAM index files (*.bam.bai) that are saved to the alignment folder. A BAM file is the compressed binary version of a SAM file that is used to represent aligned sequences.
Read collapsing adds the following BAM tags:
RX/XU: UMI combination. RX is duplicated from XU to satisfy the BAM/SAM format
XV: Number of reads in the family on one strand.
XW: Number of reads in the duplex-family or 0 if not a duplex family.
Fusion Calling
DRAGEN fusion caller performs the fusion calling. The outputs are VCF files. It first scans the genome to discover evidence of possible structural variants (SV) and large indels based on split- and spanning-reads. The evidence is enumerated into a graph with edges that connect all regions of the genome that have a possible breakend association. It then analyzes individual graph edges to discover and score SVs associated with the edges. The substeps of this process include the following items:
Inference of SV candidates associated with the edge.
Attempted assembly of the SVs breakends.
Scoring, genotyping, and filtering of the SV.
Output to VCF.
Fusion Filtering
In assays, there are hundreds to thousands of fusion candidates in a single sample. Most fusion candidates (~99%) are false positives. The fusion filtering tool, DNA Fusion Filter (DNAFF), distinguishes true fusion calls from the false positives. DNAFF performs the following functions:
Removes spurious fusions including fusions with only one supporting read pair and fusions that overlap with repeat regions, which are more likely to have sequencing errors.
Filters nonconfident supporting reads for all fusion candidates based on the following criteria:
Filter reads with low-sequence identity with the fusion contig.
Filter erroneous reads, which are reads that do not support the fusion. For example, reads that have suspicious supplementary alignment.
Deduplicate reads based on UMI information.
Applies the following rules for the final fusion output
If the fusion gene pair has been reported in COSMIC, it must have ≥ 2 supporting reads.
If the fusion gene pair has not been reported in COSMIC, it must have ≥ 3 supporting reads.
At least one fusion breakpoint must fall within the 23 target genes.
Indel Realignment and Read Stitching
The Gemini software component performs local indel realignment, paired‑read stitching, and read filtering. A stitched read is a single read that has been combined from a pair of reads. Reads near detected indels are realigned to remove alignment artifacts. The input is a single BAM file and the reference genome FASTA used to align it; the output is a corresponding single BAM file with stitched, pair‑realigned reads. Read pairs with poor map quality or supplementary and secondary alignments from the input BAM are ignored.
The following BAM tags are added to the stitched reads:
XD—Directional support string indicating forward, reverse, and stitched positions.
XR—Read pair orientation, which can be either forward-reverse (FR) or reverse-forward (RF).
Max Somatic VAF
MaxSomaticVaf is a surrogate variant for tumor fraction in ctDNA. The TMB step flags the variants by potential somatic status using database, VAF and clonal hematopoiesis information. The remaining variants are ranked by variant allele frequency in descending order. The variant allele frequency of first COSMIC hotspot (count >50) or confident somatic variant (having significantly shorter fragment size) is reported as the MaxSomaticVaf for each sample. If no such variant exists, the 4th variant is reported.
Microsatellite Instability Status (MSI)
The DRAGEN microsatellite instability status (MSI) module assesses microsatellite sites for evidence of microsatellite instability, relative to a set of baseline normal samples that are based on Jensen– Shannon (JS) distance (an information entropy metric). In total, there are 2343 selected MSI sites with 6 or 7 single nucleotide repeats in the panel. For MSI sites with 500 or more spanning duplex collapsed reads, JS distance is calculated using a test sample vs baseline normal samples, and then any two baseline normal samples. If the JS distance is significantly higher in the test sample vs baseline normal with p-value ≤ 0.01 and the JS distance difference is ≥ 0.02, the site is considered unstable. The final MSI score aggregates all JS distance across all unstable sites. The input is the BAM file from the DNA alignment and read collapsing step, and the output is an MSI metric file.
Phased Variant Calling
Scylla rapidly detects multiple nucleotide variants (MNVs) in a given sample. The software uses Scylla to detect specific, clinically relevant mutations in EGFR exon 19 and mutations in RET that would otherwise be out of scope for the variant caller. Psara filters the small variant gVCF to a small region in exon 19 of EGFR. Candidate SNPs, MNVs, and indels from this subset of the gVCF are given to Scylla along with the BAM output from Gemini. Scylla uses the original BAM to determine which of these small variants should be phased together into longer MNVs.
At a high level, Scylla identifies variants that are candidates for phasing in the input gVCF and arranges the variants into local neighborhoods. Scylla then mines the sample BAM file for any evidence that these small variants occur in the same clonal sub-populations with each other. This is done by clustering overlapping reads in the neighborhood into a minimal set of clusters, which contain the same variants.
Unlike Pisces, Scylla does not require that variants be on the same read to be phased. Once the phasing is complete, a new VCF is generated.
Small Variant Calling
DRAGEN somatic small variant caller performs somatic variant calling to identify variants in DNA samples, via local de novo assembly of haplotypes in active regions. Haplotypes are first generated with de Bruijn graph. The likelihood of a read supporting a haplotype is calculated using a Paired Hidden Markov Model. Somatic Score (SQ) is calculated as the joint posterior probability that a variant is present in the tumor sample. For each variant candidate, background noise at the same site is estimated from normal baseline samples of varying qualities. A p-value is calculated using the observed mutant depth, total depth, and background noise using binomial distribution and then converted to a variant quality score (AQ).
Variants are called if SQ >= 2 and AQ >= 20 for Catalogue of Somatic Mutations in Cancer (COSMIC) with count > 50 or SQ >= 2 and AQ >= 60 for remaining sites.
DRAGEN somatic small variant caller also performs variant merging/phasing for multiple nucleotide variant (MNVs) and complex variant detection. Variants residing on the same graph haplotype within 15 bp can be combined into one variant. The somatic variant caller generates .hard-filtered.vcf and .hard-filtered.gvcf files.
The net effect of the read collapsing and variant calling reduces false positives in a typical cell-free DNA sample from ~1500 per Mb to < 5 per Mb.
Small Variant Calling Filtering
The software component, Pepe, performs post‑processing on the small variant calling genome VCFs to polish backgrounds and adjust quality scores. The software filters out variants when error rates do not meet quality thresholds. This analysis step produces genome VCF files and associated error rate files. The minimum read depth for reference calls is 1000. The limit of detection for VAF is .05% at the minimum read depth.
Pepe computes two quality scores for each candidate that dynamically adjust for the following conditions:
Background noise
Trinucleotide change
Read support type
Filter Information for Small Variant Filtering:
.
PASS
WT
., A, C, G, etc¹
LowDP²
Reference or filtered variant candidate with depth <1000.
A, C, G, etc¹
PASS
PASS variants.
A, C, G, etc¹
LowSupport, LowVarSupport
Filtered variant candidate with low-quality score.
A, C, G, etc¹
Blacklist
Position with high background noise. Not available for variant detection.
A, C, G, etc¹
VarBias
Filtered variant candidate showed bias clustered at fragment ends.
¹ Etc refers to other variants types not mentioned in the table.
² For reference positions, a coverage below 1000X directly translates into LowDP, but if a position has a PASS filter for variant calls, LowDP is not applied. This depends on the LQ/AQ thresholds for COSMIC or non-COSMIC variants, and allele frequencies.
The net effect of the read collapsing and TMB analysis steps reduces false positives in a typical cell‑free DNA sample from ~1500 per Mb to < 5 per Mb.
For more information on DNA alignment and read collapsing, refer to Illumina DRAGEN Bio-IT Platform Product Documentation.
For each variant candidate, background noise at the same site is estimated from normal baseline samples of varying qualities. A p-value is calculated using the observed mutant depth, total depth, and background noise using binomial distribution. The p-value is then converted to a variant quality score (AQ). The sample-specific error rate of each trinucleotide change is estimated from different support categories in each sample by using all the positions with an allele frequency less than 1%. For each variant candidate, a likelihood ratio score (LQ) is computed by the corresponding error rate of the observed total and mutant read. A bias score (BFQ) is computed on each variant candidate to evaluate the imbalance of mutant vs total read support between different support groups.
For variants with a Catalogue of Somatic Mutations in Cancer (COSMIC) count > 50, the LQ and AQ thresholds are 20 and the remaining sites are 60. For indel, at least one stitched mutant support is required. For non-COSMIC variant, threshold for BFQ is < 20. In addition, positional information of mutant and WT allele in fragment will be extracted for each variant candidate. A Kolmogorov-Smirnov test will be applied to compute p-value between mutant and WT position. Variants with p-value < 0.05 and median difference > = 0.5 will be filtered and labeled VarBias. The net effect of the read collapsing and variant filtering significantly reduces false positives. For example, false positives in a typical cell-free DNA sample were reduced to < 5 per Mb from ~1500 per Mb.
In addition to the evaluation of the quality scores, certain regions covered in the product manifest are excluded from analysis due to high background noise. All excluded variants are identified in the VCF using a flag. The block list of excluded sites can be obtained on request from your local Illumina representative.
Some regions are known to be difficult to sequence. One example region is the TERT promoter region. Although sequencing can occur at the TERT promoter region, this location might result in low coverage due to the GC rich content of the sequenced region.
Tumor Mutational Burden
DRAGEN tumor mutational burden (TMB) analysis step generates TMB metrics from the annotated small variant JSON file, the callability file generated from DRAGEN metrics and intermediate files generated from small variant calling. To remove germline variants from the TMB calculation, the software uses a combination of public database filtering and a post‑database filtering strategy that uses allele frequency information and variants in close proximity. Database filtering uses the GnomAD exome, genome, and 1000 genomes database.
The TMB is calculated using the following equation:
The formula for TMB calculation is:
The eligible variants and effective panel size of the TMB calculation are detailed in the following table:
Eligible variants (numerator)
Variants in the coding region (RefSeq Cds).
Variant frequency ≥ 0.2%.
Coverage ≥ 1000X.
SNVs and Indels (MNVs excluded).
Nonsynonymous and synonymous variants.
Variants with cosmic count > 50 are excluded.
Mutations in ASXL1, DNMT3A, PPM1D, and TET2 are excluded.
Fragment-size based potential clonal hematopoiesis mutations are excluded.
Effective panel size (denominator)
Total coding region with coverage ≥ 1000X.
Outputs are captured in a *_tmb_trace.tsv file that contains the variants used in the TMB calculation, and a *.tmb.metrics.csv file that contains the TMB score calculation and configuration details.
Variant Merging
The software merges the phased variants with the other small variants generated from the small variant filtering step. In this process, exact duplicates that match chromosome, position, Ref, and Alt are removed.
Only the following Epidermal Growth Favor Receptor (EGFR) variants are added if found from Phased Variant Calling. All other EGFR variants are filtered out in variant merging.
EGFR Variants
chr7
55242482
CATCTCCGAAAGCCAACAAGGAAAT
C
chr7
55242466
GAATTAAGAGAAGCAACAT
G
chr7
55242465
GGAATTAAGAGAAG
AATTC
chr7
55242465
GGAATTAAGAGAAGCAAC
AAT
chr7
55242469
TTAAGAGAAGCAACATCTC
T
chr7
55242467
AATTAAGAGAAGCAACATC
A
chr7
55242469
TTAAGAGAAG
C
chr7
55242467
AATTAAGAGAAGCAACATC
T
chr7
55242465
GGAATTAAGA
G
chr7
55242467
AATTAAGAGAAGCAACATCTC
TCT
chr7
55242467
AATTAAGAGAAGCAAC
T
chr7
55242464
AGGAATTAAGAGAAGC
A
chr7
55242466
GAATTAAGAGAAGCAA
G
chr7
55242464
AGGAATTAAGAGA
A
chr7
55242469
TTAAGAGAAGCAA
T
chr7
55242465
GGAATTAAGAGAAGCAACATC
AAT
chr7
55242469
TTAAGAGAAGCAACATCT
CAA
chr7
55242463
AAGGAATTAAGAGAAG
A
chr7
55242468
ATTAAGAGAAGCAACATCT
A
chr7
55242462
CAAGGAATTAAGAGAA
C
chr7
55242465
GGAATTAAGAGAAGCAA
AATTC
chr7
55242469
TTAAGAGAAGCAA
C
chr7
55242467
AATTAAGAGAAGCAAC
A
chr7
55242469
TTAAGAGAAGCAACATCTCC
CA
chr7
55242468
ATTAAGAGAAG
GC
chr7
55242465
GGAATTAAGAGAAGCA
G
chr7
55242468
ATTAAGAGAAGCAAC
GCA
chr7
55242465
GGAATTAAGAGAAGCAACA
G
chr7
55249011
AC
CCAGCGTGGAT
Regions of the RET gene in exons 11, 13, and 15 are evaluated and added if detected from Phased Variant Calling.
Last updated