Hi Markus,

Thank you for your reply. I do not have the option of upgrading my current R version which is 3.5.1 and I can only install pureCN 1.12.2 on this version of R. Is it possible to work around this? I am also not clear on whether my workflow is correct. I would like to run PureCN with Mutect2, GATK4 and integrate cnvkit results with PureCN.

Currently what I do is I create a pooled reference by: 1) running each germline (unmatched but process matched) sample in Mutect2 tumor mode, and then run SomaticPanelOFNormals on this.

$GATK Mutect2 \
 -R $REF_GENOME_b37 \
 -I  ${GERMLINE_PFX}_G.bam \
 -tumor $GERMLINE_PFX\
 --germline-resource af-only-gnomad.raw.sites.vcf \
 -O  ${GERMLINE_PFX}_for_pon.vcf.gz

ls  *_pon.vcf.gz >> normals_for_pon_vcf.args

 $GATK CreateSomaticPanelOfNormals \
 -vcfs normals_for_pon_vcf.args \
 -O pon.vcf.gz

2) using bcftools merge to merge all the Mutect2 vcf files for the germlines to create one merged germline.vcf. I then compress using bgzip and index using tabix.

bcftools-1.9/bcftools merge $(ls *_pon.vcf.gz) -Oz -o pon.vcf.gz
bcftools-1.9/bcftools index -t pon.vcf.gz

3) I run this step:

Rscript $PURECN/NormalDB.R --outdir $OUT --normal_panel pon.vcf.gz --assay agilent_v6 --genome hg19 --force

4) I create a mapping bias file using an R script

 Rscript  create_mapping_bias_file.r

contents of create_mapping_bias_file.r:

library("PureCN")
normal.panel.vcf.file <- "pon.vcf.gz"
bias <- calculateMappingBiasVcf(normal.panel.vcf.file, genome = "h19")
saveRDS(bias, "mapping_bias.rds")

5) I run this command for every tumor sample:

cnvkit.py export seg ${TUMOR_PFX}_T.cns --enumerate-chroms -o ${TUMOR_PFX}.seg

5) I run mutect2 on every single tumor sample

$GATK Mutect2 \
 -R $REF_GENOME_b37 \
 -I ${TUMOR_PFX}_T.bam \
 -tumor $TUMOR_PFX\
 --panel-of-normals pon.vcf.gz \
 --germline-resource af-only-gnomad.raw.sites.vcf \
-O ${TUMOR_PFX}.vcf.gz

6) I run pureCN

Rscript $PURECN/PureCN.R \
                        --out  ${TUMOR_PFX}  \
                        --sampleid ${TUMOR_PFX}\
                        --tumor ${TUMOR_PFX}_T.cnr \
                        --segfile ${TUMOR_PFX}.seg \
                        --vcf ${TUMOR_PFX}.vcf.gz \
                        --genome hg19 \
                        --force --postoptimize --seed 123
                        --normaldb mapping_bias_agilent_v6_hg19.rds\

In the last step, somehow I get the error that the

seg file does not match the normal db

and if I remove that, I get the error above ((list) object cannot be

coerced to type 'double' Calls: runAbsoluteCN ... filterVcfMuTect2 -> filterVcfBasic -> .filterVcfByBQ -> which Execution halted).

I have searched extensively for documentation on how to integrate Mutect2 and GATK4 with PureCN while using cnvkit output and I could not get much information and so I thought documenting it here and asking for advice might help others in future. Thank you.

Hi Markus,

Thank you, I updated to GATK 4.1.9.0. Here are the steps that I have followed and the error I get:

1) I create a pooled reference by: 1) running each germline (unmatched but process matched) sample in Mutect2 tumor mode, and then run SomaticPanelOFNormals on this. Here I used the link you provided above.

For every germline sample:

$GATK Mutect2 \
 -R $REF_GENOME_b37 \
 -I ${GERMLINE_PFX}_G.bam \
 --max-mnp-distance 0 \
 -O ${GERMLINE_PFX}.vcf.gz

then:

$GATK GenomicsDBImport \
 -R $REF_GENOME_b37 \
-L $BEDFILE \
 --genomicsdb-workspace-path pon_db \
 -V ${GERMLINE_PFX}.vcf.gz \ #sample 1 germline
.
.
.
 -V${GERMLINE_PFX}.vcf.gz \ #sample 9 germline

and

$GATK CreateSomaticPanelOfNormals \
 -R $REF_GENOME_b37 \
 --germline-resource af-only-gnomad.raw.sites.vcf \
 -V gendb://pon_db \
 -O pon.vcf.gz

2) I run:

Rscript $PURECN/NormalDB.R --outdir $OUT --normal_panel pon.vcf.gz --assay agilent_v6 --genome hg19 --force

Here I get the error:

INFO [2020-11-13 12:48:51] Loading PureCN 1.21.0... INFO [2020-11-13 12:49:03] Creating mapping bias database. INFO [2020-11-13 12:49:04] Processing variants 1 to 5000... FATAL [2020-11-13 12:49:04] The normal.panel.vcf.file contains only a single sample.

FATAL [2020-11-13 12:49:04]

FATAL [2020-11-13 12:49:04] This is most likely a user error due to invalid input data or

FATAL [2020-11-13 12:49:04] parameters (PureCN 1.21.0).

Error: The normal.panel.vcf.file contains only a single sample.

This is most likely a user error due to invalid input data or parameters (PureCN 1.21.0). In addition: Warning message: In .vcf_usertag(map, tag, nm, verbose) : ScanVcfParam ‘geno’ fields not found in header: ‘AD’ Execution halted

Would you happen to know why this error occurs? Thank you.