Sample processing for CANCAN and PureCN analyses

[danielvo]

under construction

1. Normal sample pre-processing and panel of normals (PON) creation

The importance of PON data

In order to be able to interpret copy number data from a single tumor sample (i.e., without a matched control), appropriate reference data is necessary for defining the expected normal baseline. Besides containing artifacts/noise, targeted assays are also naturally biased in their coverage profile (there will be a coverage peak in the middle of a target probe location, and the coverage will gradually fall off with increasing distance from the probe center; a series of nearby probes will form patterns with hills and valleys in between). Constructing a PON for a given assay allows for:

• normalization of single-sample coverage profiles (by setting the expectations for non-aberrant genomes);
• defining a set of reliable germline variants: germline variants with stable variant allelic fractions (AFs) that can be utilized for allelic composition calculations.

CANCAN and PureCN both rely on PON data in their analyses. Despite each of them utilizing their custom PON format, the underlying data (coverage and small variant information from a selected set of normal sample) is largely the same. Below is an outline for a PON-creation process that aims to maximize the utilization of information available in TSO500 DNA samples. The approach has been previously tested at OUS.

1.1. Extension of assay target regions

There is coverage data available outside of the declared TSO500 targets, which are mostly limited to exonic regions. Utilizing a set of high-coverage control samples, the core TSO500 targets can be expanded ('grown' upstream and downstream) to cover additional bases at usable coverage (from previous experience, the number of utilized bases can nearly double if the required read depth is >= 50 bp). The extended targets will contain information about additional germline variants, providing additional variant AF data points for the copy number (CN) segmentation process. The HRD add-on (included in the TSO500 assay v2 and available separately for the TSO500 assay v1) has the potential to offer information for variant counts on a higher order of magnitude compared to the core TSO500 regions (OUS has not tested HRD region utilization in its PON efforts so far).

Previously tested target expansion process, utilizing a set S of normal samples and a selected maximum padding interval of p base-pairs:

• extend each target region with the selected padding interval in each direction (suggested software: bedtools slop)
• define the expansion-only areas by subtracting the original targets from the extended targets (bedtools subtract)
• merge overlapping expansion-only areas (bedtools merge)
• determine the base-pair level coverage within the expansion-only areas for each sample within S (bedtools coverage)
• extend the original targets based on the normal sample coverage data (custom script used previously, the original targets have been extended one bp at a time in each direction, as long as >= 70 % of the utilized normal samples maintained coverage >= 50 reads)
  • a challenge here is dealing with male vs. female samples, as the two groups have different expected coverages at certain regions of the X chromosome
• merge overlapping extended targets (bedtools merge)
• split the merged extended targets into smaller chunks, providing convenient bins for later coverage analysis (custom script used previously, with the aim to create 400 bp chunks while avoiding creation of left-over chunks smaller than 250 bp)

The whole process operates on BED files.

1.2. Normal sample GATK pre-processing/GATK PON data creation

The goal: Gather information about coverage profiles and recurrent germline variants in a given set of normals samples (S).

Previously tested process:

• convert the extended targets BED file into a GATK interval format (GATK PreprocessIntervals)
• add read group information to the normal sample BAM files (necessary for BAM files generated by the Local App) a prerequisite for running GATK tools on the data (picard AddOrReplaceReadGroups)
• collect coverage information within the extended targets for each normal sample from S (GATK CollectReadCounts)
• create a GATK coverage data PON for all samples from S (GATK CreateReadCountPanelOfNormals)
• remove noise from the normal sample coverage files utilizing the GATK coverage data PON (GATK DenoiseReadCounts)
• call germline variants for each sample within S (GATK HaplotypeCaller)
• import single-sample GVCFs into GenomicsDB before joint genotyping (GATK GenomicsDBImport)
• perform joint genotyping on samples from S (GATK GenotypeGVCFs)
• perform variant quality score recalibration for variants called within S (GATK VariantRecalibrator and GATK ApplyVQSR)
• split multi-allelic germline variant entries into single variant allele entries (bcftools norm)
• define a set of initial germline variant candidates for S: extract PASSing germline variants from the variant set generated upstream, create BED and VCF output files, with the information necessary in downstream steps (custom script)
• collect reference and alternative allele support counts for initial germline variant candidates for all samples within S (GATK CollectAllelicCounts)
• generate segmented copy ratios and minor AFs for each sample from S, utilizing denoised read counts and allelic counts (GATK ModelSegments)
• for each sample within S, filter the heterozygous germline variant AF data generated in the previous step, utilizing a simple read-count threshold in order to limit noise (custom script)

1.3. CANCAN- and PureCN-specific reference and PON data

CANCAN-specific data

Regarding coverage normalization, CANCAN can rely on the tumor sample pre-processing performed by GATK (as described in section 2. below). When it comes to germline variant AF data, CANCAN provides functionality for creating a centroid - a filtered set of germline variants with high reliability (variants with stable VAF in the normal sample set) - CANCAN can at the same time account for systematic stable variant AF shifts (e.g., when observing that a certain variant always tends to have AFs around 0.6 when being heterozygous in the normals, rather than the expected 0.5).

PureCN-specific data

PureCN provides a set of commands/tools for creating the necessary reference datasets based on the selected reference genome, the relevant assay targets and the available normal samples:

• target interval definition and pre-processing, utilizing optional data such as mappability and replication timing; the extended targets would be utilized here (PureCN IntervalFile.R)
• preparation of the NormalDB PureCN PON resource, which aggregates and normalizes the normal sample coverage and germline variant data; based on how one wishes to run PureCN during the tumor sample analysis, one could utilize either GATK-generated coverage data or PureCN's own coverage profile data; the variant input could be based either on the CANCAN-generated centroid variants or the initial germline variant candidates (PureCN NormalDB.R)

Multiple PureCN setups/analysis modes have been tested at OUS on a limited number of samples. It is difficult to say which setup is optimal for TSO500 data.

2. Tumor sample pre-processing and analysis

2.1. Tumor sample GATK pre-processing

The tumor sample needs to be pre-processed in order to generate coverage and germline variant data compatible/usable with the PON data:

• add read group information to the tumor sample BAM file (picard AddOrReplaceReadGroups)
• collect coverage information within the extended targets for the tumor sample (GATK CollectReadCounts)
• remove noise from the tumor sample coverage data utilizing the GATK coverage PON (GATK DenoiseReadCounts)
• collect reference and alternative allele support counts for initial germline variant candidates for the tumor sample (GATK CollectAllelicCounts)
• generate segmented copy ratios and minor AFs for the tumor sample, utilizing its denoised read counts and allelic counts (GATK ModelSegments)
• filter the heterozygous germline variant AF data generated in the previous step, utilizing a simple read-count threshold in order to limit noise (custom script)
• extract information about somatic variant calls available for the sample (custom script)

2.2. CANCAN and PureCN copy number analysis

Both tools need to be run using R. We can utilize standardized R scripts that perform the main CANCAN and PureCN analyses, as well as any R object -> text output conversions that we'll deem useful for later reference (e.g., generating PureCN solution overview tables and CANCAN gene-wise copy number tables).

[marrip]

wow, the PON creation is more complex than I thought but I think the listing of the different tools and subcommands helps tremendously with understanding what is going on 👍 Would it be possible get extend the tool subcommand sections with actual parameters and placeholders where necessary? Then we could use that as a blueprint for a PON pipeline that can be part of the TSO500 universe. Just to clarify, this process is independent of dragen v2 and LocalApp, right? We would run the same for both TSO500 v1 and v2 samples?

[danielvo]

As everything described above has been previously run on consented samples within TSD, scripts with complete commands (as well as all outputs) are available there - I can share everything with whoever implements the preprocessing in tsoppy (so they won't need to look up things here). :-) Should I include additional details here anyway? I thought this thread could provide a general commented overview of the process, containing the necessary information for allowing a discussion in case we'd like to tweak the pipeline.

I am not sure about the Dragen/LocalApp independence. The two pipelines report very different coverages.* Unless everything is just proportionately scaled down by one of the pipelines compared to the other, we will need separate PONs. We have the LocalApp-based PON available inside TSD for 166 normals; we could run Dragen on them all and create a Dragen-based PON from the same sequence data there. V2 samples contain the additional HRD targets, so we'd definitely want to have a separate PON for those - that one will only be for Dragen data, as the LocalApp does not create output for HRD regions.

* In case of the one test sample (TVD0016*), both pipelines report the same starting total read count, quite similar median target coverages (947 vs. 918), but a big difference in a random KRAS position that I looked up (1263 vs. 887). I therefore suspect that separate LocalApp and Dragen PONs will be necessary.

[marrip]

Ah ok, no, if you already have it documented somewhere it should be fine. Whoever implements should get access, I agree. But the pipeline will not be in tsoppy but rather its own thing and tsoppy will maybe carry potential python scripts that the pipeline needs. Maybe I should talk about the differences in the work meeting. I feel like it might not be clear yet what is covering what. My apologies if I was unclear.

For me it was more about are the processing steps for LocalApp and Dragen different due to the difference in coverage or can we run the same steps for both? What would be different? I didn't want to mix the data in any way. Sorry, my question might have been ambiguous.

[danielvo]

Yes, let's discuss the vision. :-D

We should be able to run the same steps for all scenarios, as far as I am aware. For the V2 assay we will need to provide a different target list (taking into the account the HDR probe sites). The Dragen DNA BAM files already have read groups, but the AddOrReplaceReadGroups tools could likely be used uniformly on data from both pipelines (adding read groups in case of the LocalApp data and replacing read groups in case of Dragen data). The expected sample data input would always be BAM files and somatic variant calls (we should have the latter available once we parse the small variant input).

[danielvo]

@marrip How do we go about implementing this, Martin? I assume the different processes (e.g., tumor-sample preprocessing) will be workflows in NextFlow, with some user inputs that are being channeled through standard external tools most of the time. We didn't get any feedback from Arne yet regarding his wishes for the CANCAN Docker image placement etc. I guess a NextFlow workflow for PON creation wouldn't necessary fit into the current CANCAN repository anyway.. so we can assume it being a tso500_v2_post_processing ("TSOPPI_2"?) workflow. The necessary custom scripts will be a part of tsoppy, I assume.

[marrip]

exactly, the custom scripts should go into tsoppy. I think PON creation should get its own small workflow as it is more a one-off job (at least in my understanding). The created PON can then be put where the user wants and be accessed during the tso500_v2_post_processing. Also, we can combine creation of CANCAN and PureCN PON creation and potentially other future PONs and have switches to decide which ones we need.