tsoppy input standardization - 2. CNV data

[danielvo]

Inputs considered so far:

• Dragen TSO500 analysis pipeline v. 2.6.2
  • the information below is based on output generated for samples sequenced with TSO500 assay version 2, but the current
  • unclear: is the Dragen server version a factor?
    • OUS assumes it does not, and both OUS and UNN use the same Dragen server version.
      • It looks like both assay v1 and assay v2 samples can now be analyzed with the same Dragen software.
    • Dragen server version [07.021.624.]3.11.2 was used when generating the output referenced in the text below
    • UNN utilizes the same Dragen server version and Dragen TSO500 analysis pipeline version
    • unclear: When analyzing samples with ICA (Illumina Cconnected Analysis), can one pick exactly the same Dragen server version and TSO500 pipeline version?
  • unclear: can the analysis output structure vary based on the sample sheet content/format choices?
    • different Pair_ID and Sample_ID values were used when generating the output referenced in the text below
    • we should aim to standardize the sample sheets across InPreD nodes either way
• Local App version 2.2.0.12
  • the sample sheet content/format (e.g., the Pair_ID value choice) can impact the analysis output structure
    • different Pair_ID and Sample_ID values were used when generating the output referenced in the text below

2. CNV data

2.1. Input files

2.1.1. Target-wise read coverage data

Dragen TSO500 2.6.2:

• file path for raw read counts: Logs_Intermediates/DnaDragenCaller/<sample_ID>/<sample_ID>.target.counts.gz

#TARGET COUNTS FILE
#Format: raw coverage
#Input alignment file: NA
#Input target file: resource_bundle/cnv/TSO500_DRAGEN_BASELINE_SOLID_matched_targets_v2.combined.counts.txt.gz
contig  start   stop    name    TVD0016-D01-P01-A00     improper_pairs
chr1    2488102 2488174 TNFRSF14-chr1.0 1778    6
chr1    2489163 2489275 TNFRSF14-chr1.1 2904    5
<...>

• file path for normalized read counts: Logs_Intermediates/DnaDragenCaller/<sample_ID>/<sample_ID>.tn.tsv.gz

#title = Tangent normalized coverage profile
#sex = MALE
contig  start   stop    name    TVD0016-D01-P01-A00     improper_pairs
chr1    2488102 2488174 TNFRSF14-chr1.0 0.27543534878419546     6
chr1    2489163 2489275 TNFRSF14-chr1.1 0.46692604351802625     5
<...>

Local App 2.2.0.12:

• file path for raw and normalized read counts: Logs_Intermediates/CnvCaller/<sample_ID>/<sample_ID>_normalizedBinCount.tsv

Chr     Start   End     GC      Gene    binName Raw     NormalizedCount
chr1    2488102 2488174 0.643835616438356       TNFRSF14        chr1.0  1149    0.968606613912893
chr1    2489163 2489275 0.646017699115044       TNFRSF14        chr1.1  1674    1.15388179631073
<...>

2.1.2. Gene-level coverage-based data (fold changes)

Dragen TSO500 2.6.2:

• file path: Logs_Intermediates/DnaDragenCaller/<sample_ID>/<sample_ID>.cnv.vcf.gz
• information about 500 genes

##fileformat=VCFv4.2
<...>
##ALT=<ID=CNV,Description="Copy number variant region">
##ALT=<ID=DEL,Description="Deletion relative to the reference">
##ALT=<ID=DUP,Description="Region of elevated copy number relative to the reference">
##INFO=<ID=REFLEN,Number=1,Type=Integer,Description="Number of REF positions included in this record">
##INFO=<ID=SVLEN,Number=.,Type=Integer,Description="Difference in length between REF and ALT alleles">
##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=CIPOS,Number=2,Type=Integer,Description="Confidence interval around POS">
##INFO=<ID=CIEND,Number=2,Type=Integer,Description="Confidence interval around END">
##INFO=<ID=SEGID,Number=1,Type=String,Description="Segment ID">
##FILTER=<ID=cnvQual,Description="CNV with quality below 0">
##FILTER=<ID=cnvBinSupportRatio,Description="CNV with low supporting number of bins with respect to event length">
##FILTER=<ID=cnvCopyRatio,Description="CNV with copy ratio within +/- 0 of 1.0">
##FILTER=<ID=LoDFail,Description="CNV is closer to neutral PON counts than the specified limit of detection threshold (PON * (1 +/- 0.2))">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=SM,Number=1,Type=Float,Description="Linear copy ratio of the segment mean">
##FORMAT=<ID=BC,Number=1,Type=Integer,Description="Number of bins in the region">
##FORMAT=<ID=PE,Number=2,Type=Integer,Description="Number of improperly paired end reads at start and stop breakpoints">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	TVD0016-D01-P01-A00
chr1	2489163	DRAGEN:REF:chr1:2489163-2494714	N	.	3	PASS	END=2494714;REFLEN=5552;SEGID=TNFRSF14	GT:SM:BC:PE	./.:1.11059:7:5,14
chr1	8073268	DRAGEN:REF:chr1:8073268-8085120	N	.	3	PASS	END=8085120;REFLEN=11853;SEGID=ERRFI1	GT:SM:BC:PE	./.:0.966845:5:68,11
<...>

Local App 2.2.0.12:

• file path Logs_Intermediates/CnvCaller/<sample_ID>/<sample_ID>_CopyNumberVariants.vcf
• information about 59 genes

##fileformat=VCFv4.1
<...>
##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=ANT,Number=1,Type=String,Description="Gene associated with copy number variant">
##FORMAT=<ID=FC,Number=1,Type=Float,Description="Copy number variant fold change, calculated as ratio of sample to control">
##FILTER=<ID=LowValidation,Description="Algorithm for copy number loss has been validated using only simulation datasets">
##ALT=<ID=DEL,Description="Identified deletion in target gene">
##ALT=<ID=DUP,Description="Identified duplication in target gene">
<...>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	TVD0016-D01-P01-A00
chr1	40356095	.	C	.	5	PASS	END=40372764;ANT=MYCL	FC	1.000
chr1	115245084	.	T	.	28	PASS	END=115261621;ANT=NRAS	FC	0.910
<...>

• file path: Logs_Intermediates/CnvCaller/<sample_ID>/<sample_ID>_foldChange.tsv
• information about 521 genes

# Gender = Male
# Gender score = 0.474
GeneName        FC      qScore
TNFRSF14        1.02589443831505        0.468162465039736
ERRFI1  0.958145227608634       0.3269275059854
<...>

2.1.3. Gene-level copy number data

• Only generated by Dragen, and only for samples with HRD probes
  • OUS assumes that the results will be compatible for assay v2 samples and assay v1 samples that utilized the HRD add-on (we can run the Dragen TSO500 2.6.2 analysis on old v1 samples that utilized the HRD add-on and check the assumption)

Dragen TSO500 2.6.2:

• file path: Logs_Intermediates/Gis/<sample_ID>/<sample_ID>.abcn_genes.tsv
• CN stands for (total) copy number
• MCN standas for minor copy number

#Chromosome     start   end     gene    CN      MCN
chr1    2489162 2494714 TNFRSF14        2       1
chr1    8073267 8085120 ERRFI1  2       1
<...>

• file path: Logs_Intermediates/Gis/<sample_ID>/<sample_ID>_bAllele.tsv
• in case of the available AcroMetrix test sample, information for 23617 sites is included in the output (each site is defined by an rs/ss ID)
• this data can be directly utilized for extended BAF plotting

<...>
chr1    1868712 rs2748970       0.517961
chr1    1967954 rs2376803       1
chr1    2046938 rs3128339       0
<...>

• file path: Logs_Intermediates/Gis/<sample_ID>/<sample_ID>_logRatio.tsv
• in case of the available AcroMetrix test sample, the following site information is included
  • sites with 23616 rs/ss IDs that were also captured in the *bAllele.tsv file mentioned above
  • additional sites defined by rs IDs (24600 rs/ss variant sites in total)
  • 18400 additional elements (IDs named after WES targets and specific genes
    • these are within the standard TSO500 targets

chr1    2488119 2488120 target-wes-chr1-2488101:2488138/1       -0.0224617
<...>
chr1    1868712 1868713 rs2748970       -0.00351495
chr1    1967954 1967955 rs2376803       -0.260433
chr1    2046938 2046939 rs3128339       0.19868
<...>
chr17   41261513        41261514        BRCA1_336       0.183409
<...>

2.2. Notable differences between Dragen and Local App output files

• The target-level read normalization differs between Dragen and Local App (the "normal" baseline value seemingly being 0 and 1, respectively)
• Some of the key information is located in different output files
  • The patient sex ("Gender") is reported in different output files (target-wise vs. gene-level)
  • The complete gene-wise information is available in Dragen's final VCF and Local App's intermediate TSV file
  • "Fold Change" nomenclature: Dragen's "Linear copy ratio of the segment mean" corresponds to Local App's "Copy number variant fold change, calculated as ratio of sample to control"

2.3. CNV data usage in tsoppy modules

2.3.1. Modules that would utilize the directly available CNV data

• NTRC (Normalized Tumor Read Coverage) plots
• Genome-wide/chromosome-wise CNV plots

2.3.2. Modules that require additional pre-processing

• PureCN and CANCAN allele-specific CN analysis
• Small variant table/predisposition variant table modules (a gene's copy number, as well as expected germline/somatic VAFs could be reported based on CANCAN's allele-specific CNV output)

2.3.3. Suggestions for implementation

• The Local App analysis utilizes coverage data normalized against internal controls (2.1.1. and 2.1.2.) in order to infer calculate fold changes that can converted into total copy number if tumor content and tumor ploidy are known. This approach has been utilized by TSOPPI and its functionality and corresponds to modules named in section 2.3.1.
  • The same coverage profiles along the genome (and sorted lists of genes according to observed FC values, etc.) could be created also from Dragen output, but the Dragen analysis pipeline apparently employs a different coverage normalization method on the target level.
    • It needs to be tested what the plots would look like for Dragen's "Tangent normalized coverage profile" values when compared to plots generated for Local App's log(NormalizedCount) values.
      • OUS can do this exercise.
• Local App-generated FC values and Dragen's "Linear copy ratio of the segment mean" should correspond to each other, they just need to be extracted from different files/file formats.
• For samples that utilize the HRD probes (this includes assay v2 samples), total copy number and minor copy number results are available on gene level.
  • Values for a single solution are presented (including purity and ploidy).
    • Coverage ration and VAF information is provided for a large number of sites, and these sites are not excluded from the one BAM file generated for DNA samples by the Dragen pipeline.
      • Custom analysis could take advantage of the data in extended CANCAN/PureCN analysis.

[tinavisnovska]

Hi @danielvo, regarding 2.3. CNV data usage in tsoppy modules: Do we still want to keep the possibility open for distributing cnv data via variants for interpretation / predisposition table similarly to what https://github.com/InPreD/TSOPPI_code/blob/main/user_scripts/libs/05_PCGR_to_variant_interpretation_table.py does? If so, maybe we can add that usecase into the 2.3 section in this document too.

[danielvo]

Hello @tinavisnovska - let me add it to the list above. Having CNV information in the small variant tables was my own initiative (it wasn't based on a request from molecular biologists), so I am not sure whether anyone has found it useful in practice, even though it sounded like a good idea. We could ask the interpretation teams for feedback on this. I think the feature would be more useful if it was based on a proper ploidy/purity estimate (one could use CANCAN's output instead of the Local App/Dragen output perhaps, even though nothing will be reliable 100 % of the time). I will thus put it under 2.3.2. for now.

[marrip]

Files required for NRC data based on LocalApp output:

• Logs_Intermediates/CnvCaller/*/*_normalizedBinCount.tsv
• Logs_Intermediates/CnvCaller/*/*_foldChange.tsv
• default_NRC_gene_list.tsv
• valid_NRC_gene_targets.tsv
• TST500C_manifest.bed
• PAR_genes.tsv

We would need to define all files from both sources needed for format_NRC_data.py to get the table necessary for plotting. For me it would be useful to get List of Dragen input and LocalApp input plus list of reference files that should go in and then a clearly defined table format. Also how to convert values from one source or the other to match expected parameter. Should we place the reference bed and tsv files into our reference repo? If they are small enough, it would be nice to have them available. Or are they illumina owned?

[danielvo]

• TST500C_manifest.bed is from Illumina, defining the TSO500 DNA assay target regions.
• PAR_genes.tsv is a custom list of genes that are on both the X and Y chromosomes, so having two copies despite being located on chromosomes that males have only one copy of.
• default_NRC_gene_list.tsv and valid_NRC_gene_targets.tsv are custom lists for deciding what will be plotted.

This leaves *_normalizedBinCount.tsv and *_foldChange.tsv files as the necessary Local App outputs. Dragen unfortunately uses a different coverage normalization method, and it's unclear to me as of now whether it's possible to simply convert one to the other. Judging by the DNA sample BAM file headers generated by the Local App and Dragen, different reference genomes might be used for DNA alignment by the two pipelines, which makes the coverage profiles compatibility questionable (analyses run on the same FASTQ files result in different reported coverages, but the full reasoning for this is also unclear to me now).

It looks like fold change value equivalents can be easily recovered for 500 genes from both pipelines, albeit from different sources (Logs_Intermediates/DnaDragenCaller/<sample_ID>/<sample_ID>.cnv.vcf.gz and Logs_Intermediates/CnvCaller/<sample_ID>/<sample_ID>_foldChange.tsv).

[marrip]

• where did we get PAR_genes.tsv from? Is it something we extracted or is it part of the illumina reference bundle?
• How were the custom lists constructed?
• Is the bed file freely available on illumina's website?

[danielvo]

The PAR gene list is manually extracted from a publication referenced in the PAR_genes.tsv file header. I can update #21 regarding the custom list construction. I am not aware of the target BED file being made available via Illumina's website.

[marrip]

Is 06_plot_CNVs.R this the CNV plotting script for "Genome-wide/chromosome-wise CNV plots"? If yes, here are the LocalApp inputs:

• Logs_Intermediates/CnvCaller/*/*_CopyNumberVariants.vcf
• Logs_Intermediates/CnvCaller/*/*_normalizedBinCount.tsv
• Logs_Intermediates/CnvCaller/*/*_foldChange.tsv
• merged_normal_sample_CNV_data.tsv
• PAR_genes.tsv
• CnvGeneThresholds.csv
• centromere_data/hg19_chromosome_sizes.tsv

Then 03_TSO_500_copy_number_aggregation.py is the script we need to consider the input and table format for similar to NRC.

[danielvo]

Yes, the 06_plot_CNVs script generates the Genome-wide/chromosome-wise CNV plots.

• merged_normal_sample_CNV_data.tsv are coverage values observed in InPreD's normal samples, it looks like we will need to generate a Dragen-specific version of the file
• PAR_genes.tsv is a custom list of genes that are on both the X and Y chromosomes (having two copies despite being located on chromosomes that males have only one copy of)
• CnvGeneThresholds.csv is a Local App resource with Illumina's thresholds for deletions and amplifications - it's only used to highlight the non-abnormal value ranges (according to Illumina, before any adjustments..) in some plots; Dragen will have its own counterpart, quite possibly different from the Local App one; if we want to refer to the Illumina pipeline thresholds in our text/graphical output, we will want access to the values
• centromere_data/hg19_chromosome_sizes.tsv is a custom resources that helps plotting the whole genome as a single continuous plot

That leaves us with the Illumina pipeline files for individual samples, with normalized target-wise coverage values and gene-level fold changes (as discussed above).

The 03_TSO_500_copy_number_aggregation.py script takes all the CNV/coverage data for both the tumor (and a matched normal, if available), it calculates tumor content adjusted values and prints a table with details over everything down to target level. The 06_plot_CNVs script can then extract the necessary values for the individual plots it generates.

[marrip]

As discussed previously, for the PON we would need a description on how it is created so we can develop a workflow for it. That would be the best so we can basically redo the PON whenever we like with any samples we choose.

The chromosome sizes, how was that created? Should probably go into reference repo as well.

[danielvo]

Agreed.

I have started describing the PON-creation process here: #25 .
The chromosome sizes file puts together some public information (from genome reference files and the UCSC Genome Browser) utilizing custom scripts (these are actually also available in the TSOPPI resources directory). We can have both the source files and the scripts (after any adjustments that you would deem appropriate) in the reference repo, yes.

[marrip]

I would keep scripts out of the reference repo as it should be only for the actual files and metadata but I would reference the scripts and commands in the Method section of the README.md.

[marrip]

For PureCN, we only need:

• Logs_Intermediates/StitchedRealigned/*/*.bam
• Logs_Intermediates/Tmb/*/*_TMB_Trace.tsv
• PON and reference files

which Dragen files are corresponding, do they differ and, if yes, how?

Also, what format does the input vcf to PureCN need to?

[tinavisnovska]

• *TMB_trace.tsv files are compared in 5. Small variants, in 5.1.2
• *.bam files in 4. dna and rna alignment data

[danielvo]

When running PureCN and CANCAN, one wishes to maximize the amount of usable coverage and germline variant AF data points. It's therefore desirable to do some custom pre-processing on extended assay targets utilizing tools like GATK (one can then roughly double the number of utilized genomic positions with the Local App generated sample output). Assay v2 (and the HRD-add on for assay v1) in addition target over 20 000 SNP locations, which can add an order of magnitude more variant AF and coverage data points. I have been focusing on recurrent germline variants within our limited panel of normals when setting up these custom analyses, so I don't believe the data would be useful in somatic analysis context other than copy number determination (which makes me think that PureCN/CANCAN pre-processing and analysis are largely a separate analysis task - not a simple utilization of the Local App/Dragen output).

[tinavisnovska]

Does your text in the last bracket have any possible implications for whether to include purecn/cancan analysis into tsoppy?

[danielvo]

I meant it more as a distinction along the lines of sections 2.3.1. and 2.3.2. above - PureCN and CANCAN won't simply utilize TSV files already made available by the Local App and Dragen, they need to generate own small variant and coverage data from the BAM files. I do think that we want them to be part of the tsoppy NextFlow pipeline, in order to have their results available when tsoppy is done.