Describe the issue
I have WES data from cancer cell lines established from tumor samples in our lab. I dont have corresponding, matched normal healthy samples, but we have sequenced a number of non-malignant, cancer associated fibroblasts, that I use as process-matched normal samples. I realized that in the genes.tsv files there is not a single gene with a value of 0 in the C column (the minimum value is 0.42; all other values are > 0.5).
Also, in the "type" column of genes.tsv files, there is just one deletion, in all samples.
table(tib_v_tumorSampleId_c_genes$type, useNA = "ifany")
# AMPLIFICATION DELETION <NA>
# 6389 1 979333
Related to that: I know from some cell lines based on aCGH data that they have a homozygous deletion in a gene, but when I look for that gene in genes.tsv file of the respective samples, this gene is missing in the genes.tsv file.
To Reproduce
This is the PureCN.R command:
Rscript "$purecn/PureCN.R" \
--out "$outdir" \
--tumor "$cov" \
--sampleid "$s" \
--vcf "$VCF" \
--normaldb "$normaldb" \
--mapping-bias-file "$mapping_bias" \
--intervals "$intervals" \
--snp-blacklist "$SNP_BLACKLIST" \
--genome hg38 \
--min-base-quality 20 \
--min-purity 0.9 \
--max-purity 0.99 \
--max-copy-number 8 \
--cosmic-cnt-info-field "CNT" \
--cosmic-vcf-file "$COSMIC_VCF" \
--min-cosmic-cnt 6 \
--fun-segmentation PSCBS \
--model-homozygous \
--post-optimize \
--max-segments 500 \
--out-vcf \
--cores 4 \
--force \
--seed 123
Expected behavior
I would have expected more homozygous deletions, based on previous aCGH from replicates.
Log file
cat ./results/VUC133_946_T/VUC040_R1_TI_C/VUC040_R1_TI_C.log
INFO [2026-03-01 11:20:55] ------------------------------------------------------------
INFO [2026-03-01 11:20:55] PureCN 2.14.1
INFO [2026-03-01 11:20:55] ------------------------------------------------------------
INFO [2026-03-01 11:20:55] Arguments: -tumor.coverage.file /scratch/rgraus29/projects/031_wesB1_b2_afterSarek/analysis/rgraus29/issue7_sarekDownstreamMutectNew/022_cl_pureCn_VucVufsCemm/coverage_files/cellLines/VUC040_R1_TI_C.recal_coverage_loess.txt.gz -log.ratio -seg.file -vcf.file ../021_sarekMutectCemmVucVufsTumOnly/results/sarek_VucVufCemmTumOnlyPureCnSet/variant_calling/mutect2/VUC040_R1_TI_C/VUC040_R1_TI_C.mutect2.filtered.vcf.gz -genome hg38 -sex ? -args.setPriorVcf 6 -args.setMappingBiasVcf reference_files/mapping_bias_twist_hg38.rds -args.filterIntervals 100,0.05 -args.segmentation 0.005,NULL, -sampleid VUC040_R1_TI_C -min.ploidy 1.4 -max.ploidy 6 -test.purity 0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99 -max.non.clonal 0.2 -max.homozygous.loss 0.05,1e+07 -log.ratio.calibration 0.1 -model.homozygous TRUE -error 0.001 -interval.file ../014_cl_pureCn_muParamPaddingVarNotInDb/reference_files/twist_baits_hg38_intervals_p100_noOffTarget.txt -min.logr.sdev 0.15 -max.segments 500 -plot.cnv TRUE -vcf.field.prefix PureCN. -cosmic.vcf.file -DB.info.flag DB -POPAF.info.field POP_AF -Cosmic.CNT.info.field Cosmic.CNT -model beta -post.optimize TRUE -log.file results/VUC133_946_T/VUC040_R1_TI_C/VUC040_R1_TI_C.log -normal.coverage.file <data> -normalDB <data> -args.filterVcf <data> -fun.segmentation <data> -test.num.copy <data> -speedup.heuristics <data> -BPPARAM <data>
INFO [2026-03-01 11:20:55] Using BiocParallel for parallel optimization.
INFO [2026-03-01 11:20:55] Loading coverage files...
INFO [2026-03-01 11:21:00] Mean target coverages: 103X (tumor) 100X (normal).
INFO [2026-03-01 11:21:01] Mean coverages: chrX: 89.59, chrY: 3.31, chr1-22: 101.94.
INFO [2026-03-01 11:21:01] Mean coverages: chrX: 120.71, chrY: 3.17, chr1-22: 101.45.
INFO [2026-03-01 11:21:13] Removing 13451 intervals with missing log.ratio.
INFO [2026-03-01 11:21:14] Removing 3 low/high GC targets.
INFO [2026-03-01 11:21:14] Removing 1745 intervals excluded in normalDB.
INFO [2026-03-01 11:21:14] Removing 5 intervals with low total coverage in normal (< 150.00 reads).
INFO [2026-03-01 11:21:14] normalDB provided. Setting minimum coverage for segmentation to 0.0015X.
INFO [2026-03-01 11:21:14] Removing 80 low count (< 100 total reads) intervals.
INFO [2026-03-01 11:21:14] Removing 9 non-diploid allosome intervals.
INFO [2026-03-01 11:21:14] Using 232317 intervals (232317 on-target, 0 off-target).
INFO [2026-03-01 11:21:15] No off-target intervals. If this is hybrid-capture data, consider adding them.
INFO [2026-03-01 11:21:16] AT/GC dropout: 1.01 (tumor), 1.06 (normal), 0.97 (coverage log-ratio).
WARN [2026-03-01 11:21:16] High GC-bias in normalized tumor vs normal log2 ratio.
INFO [2026-03-01 11:21:16] Loading VCF...
INFO [2026-03-01 11:21:19] Found 75993 variants in VCF file.
INFO [2026-03-01 11:21:19] Removing 2564 triallelic sites.
INFO [2026-03-01 11:21:19] Maximum of POPAF INFO is > 1, assuming -log10 scaled values
WARN [2026-03-01 11:21:19] vcf.file has no DB info field for membership in germline databases. Found and used valid population allele frequency > 0.001000 instead.
INFO [2026-03-01 11:21:20] 36609 (49.9%) variants annotated as likely germline (DB INFO flag).
WARN [2026-03-01 11:21:21] Found 1 variants with missing allelic fraction starting with chr20:2652734_G/A. Removing them.
INFO [2026-03-01 11:21:21] VUC040_VUC040_R1_TI_C is tumor in VCF file.
INFO [2026-03-01 11:21:23] 545 homozygous and 76 heterozygous variants on chrX.
INFO [2026-03-01 11:21:23] Sex from VCF: M (Fisher's p-value: < 0.0001, odds-ratio: 4.77).
WARN [2026-03-01 11:21:23] Sex mismatch of coverage and VCF. Could be because of noisy data, contamination, loss of chrY or a mis-alignment of coverage and VCF.
INFO [2026-03-01 11:21:23] Detected MuTect2 VCF.
INFO [2026-03-01 11:21:23] Removing 28683 Mutect2 calls due to blacklisted failure reasons.
INFO [2026-03-01 11:21:24] Removing 48 low quality variants with non-offset BQ < 20.
INFO [2026-03-01 11:21:24] Base quality scores range from 19 to 42 (offset by 1)
INFO [2026-03-01 11:21:24] Minimum number of supporting reads ranges from 2 to 27, depending on coverage and BQS.
INFO [2026-03-01 11:21:28] Initial testing for significant sample cross-contamination: unlikely
INFO [2026-03-01 11:21:28] Removing 7571 variants with AF < 0.030 or AF >= Inf or insufficient supporting reads or depth < 15.
INFO [2026-03-01 11:21:36] Removing 2160 blacklisted variants.
INFO [2026-03-01 11:21:36] Total size of targeted genomic region: 70.11Mb (86.74Mb with 50bp padding).
INFO [2026-03-01 11:21:37] 14.1% of targets contain variants.
INFO [2026-03-01 11:21:37] Removing 675 variants outside intervals.
INFO [2026-03-01 11:21:37] Setting somatic prior probabilities for likely germline hits to 0.000500 or to 0.500000 otherwise.
INFO [2026-03-01 11:21:37] Loading mapping bias file mapping_bias_twist_hg38.rds...
INFO [2026-03-01 11:21:38] Found 709342 variants in mapping bias file.
INFO [2026-03-01 11:21:42] Imputing mapping bias for 732 variants...
INFO [2026-03-01 11:21:55] Excluding 10267 novel or poor quality variants from segmentation.
INFO [2026-03-01 11:21:55] Sample sex: F
INFO [2026-03-01 11:21:55] Segmenting data...
INFO [2026-03-01 11:21:55] Interval weights found, will use weighted PSCBS.
INFO [2026-03-01 11:21:55] MAPD of 23990 allelic fractions: 0.00 (0.03 adjusted).
INFO [2026-03-01 11:21:56] Setting undo.SD parameter to 0.750000.
INFO [2026-03-01 11:22:51] Found 713 segments, exceeding max.segments threshold of 500.
INFO [2026-03-01 11:22:51] Setting undo.SD parameter to 1.125000.
INFO [2026-03-01 11:23:37] Setting prune.hclust.h parameter to 0.200000.
INFO [2026-03-01 11:23:39] Found 419 segments with median size of 0.56Mb.
INFO [2026-03-01 11:23:39] Using 34291 variants.
INFO [2026-03-01 11:23:40] Mean standard deviation of log-ratios: 0.23 (MAPD: 0.18)
INFO [2026-03-01 11:23:40] 2D-grid search of purity and ploidy...
INFO [2026-03-01 11:23:49] Local optima: 0.9/4.4, 0.97/5.6, 0.9/2.2, 0.9/3.4, 0.9/1.4
INFO [2026-03-01 11:23:50] Testing local optimum 1/5 at purity 0.90 and total ploidy 4.40...
INFO [2026-03-01 11:23:50] Testing local optimum 2/5 at purity 0.97 and total ploidy 5.60...
INFO [2026-03-01 11:23:50] Testing local optimum 3/5 at purity 0.90 and total ploidy 2.20...
INFO [2026-03-01 11:23:50] Testing local optimum 4/5 at purity 0.90 and total ploidy 3.40...
INFO [2026-03-01 11:23:50] Testing local optimum 5/5 at purity 0.90 and total ploidy 1.40...
INFO [2026-03-01 11:23:55] Recalibrating log-ratios...
INFO [2026-03-01 11:23:55] Testing local optimum 2/5 at purity 0.97 and total ploidy 5.60...
INFO [2026-03-01 11:23:55] Recalibrating log-ratios...
INFO [2026-03-01 11:23:55] Testing local optimum 5/5 at purity 0.90 and total ploidy 1.40...
INFO [2026-03-01 11:24:01] Recalibrating log-ratios...
INFO [2026-03-01 11:24:01] Testing local optimum 5/5 at purity 0.90 and total ploidy 1.40...
INFO [2026-03-01 11:24:08] Recalibrating log-ratios...
INFO [2026-03-01 11:24:08] Testing local optimum 5/5 at purity 0.90 and total ploidy 1.40...
INFO [2026-03-01 11:24:14] Recalibrating log-ratios...
INFO [2026-03-01 11:24:14] Testing local optimum 4/5 at purity 0.90 and total ploidy 3.40...
INFO [2026-03-01 11:24:17] Recalibrating log-ratios...
INFO [2026-03-01 11:24:17] Testing local optimum 2/5 at purity 0.97 and total ploidy 5.60...
INFO [2026-03-01 11:24:39] Recalibrating log-ratios...
INFO [2026-03-01 11:24:39] Testing local optimum 2/5 at purity 0.97 and total ploidy 5.60...
INFO [2026-03-01 11:24:40] Recalibrating log-ratios...
INFO [2026-03-01 11:24:40] Testing local optimum 4/5 at purity 0.90 and total ploidy 3.40...
INFO [2026-03-01 11:24:46] Recalibrating log-ratios...
INFO [2026-03-01 11:24:46] Testing local optimum 4/5 at purity 0.90 and total ploidy 3.40...
INFO [2026-03-01 11:25:03] Skipping 2 solutions exceeding max.non.clonal (0.20): 0.99/1.1, 0.9/1.33 (purity/tumor ploidy)
INFO [2026-03-01 11:25:03] Fitting variants with beta model for local optimum 1/5...
INFO [2026-03-01 11:25:04] Fitting variants with beta model for local optimum 2/5...
INFO [2026-03-01 11:25:04] Fitting variants with beta model for local optimum 3/5...
INFO [2026-03-01 11:25:05] Fitting variants for purity 0.95, tumor ploidy 3.69 and contamination 0.01.
INFO [2026-03-01 11:25:05] Fitting variants for purity 0.97, tumor ploidy 4.93 and contamination 0.01.
INFO [2026-03-01 11:25:06] Fitting variants for purity 0.90, tumor ploidy 2.43 and contamination 0.01.
INFO [2026-03-01 11:27:35] Fitting variants for purity 0.91, tumor ploidy 3.69 and contamination 0.01.
INFO [2026-03-01 11:27:38] Rare karyotype solution. Skipping post-optimization.
INFO [2026-03-01 11:27:38] Optimized purity: 0.97
INFO [2026-03-01 11:27:39] Fitting variants for purity 0.91, tumor ploidy 2.43 and contamination 0.01.
INFO [2026-03-01 11:30:06] Fitting variants for purity 0.92, tumor ploidy 3.69 and contamination 0.01.
INFO [2026-03-01 11:30:09] Fitting variants for purity 0.92, tumor ploidy 2.43 and contamination 0.01.
INFO [2026-03-01 11:32:38] Fitting variants for purity 0.93, tumor ploidy 3.69 and contamination 0.01.
INFO [2026-03-01 11:32:41] Fitting variants for purity 0.93, tumor ploidy 2.43 and contamination 0.01.
INFO [2026-03-01 11:35:06] Fitting variants for purity 0.94, tumor ploidy 3.69 and contamination 0.01.
INFO [2026-03-01 11:35:10] Fitting variants for purity 0.94, tumor ploidy 2.43 and contamination 0.01.
INFO [2026-03-01 11:37:38] Fitting variants for purity 0.96, tumor ploidy 3.69 and contamination 0.01.
INFO [2026-03-01 11:37:42] Optimized purity: 0.94
INFO [2026-03-01 11:40:09] Fitting variants for purity 0.97, tumor ploidy 3.69 and contamination 0.01.
INFO [2026-03-01 11:42:37] Fitting variants for purity 0.98, tumor ploidy 3.69 and contamination 0.01.
INFO [2026-03-01 11:45:07] Fitting variants for purity 0.99, tumor ploidy 3.69 and contamination 0.01.
INFO [2026-03-01 11:47:38] Optimized purity: 0.97
INFO [2026-03-01 11:47:38] Done.
INFO [2026-03-01 11:47:38] ------------------------------------------------------------
INFO [2026-03-01 11:47:43] Generating output files...
INFO [2026-03-01 11:47:59] Tumor/normal noise ratio: 13.362
WARN [2026-03-01 11:47:59] Extensive noise in tumor compared to normals.
INFO [2026-03-01 11:47:59] Using BiocParallel for parallel optimization.
B-allele frequency plot
This is the screenshot from the plots on page 2 of the sample.pdf file of the most likely solution:
Session Info
R version 4.5.1 (2025-06-13)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] qs2_0.1.4 SummarizedExperiment_1.38.1
[3] Biobase_2.68.0 GenomicRanges_1.60.0
[5] GenomeInfoDb_1.44.3 IRanges_2.42.0
[7] S4Vectors_0.48.0 BiocGenerics_0.54.1
[9] generics_0.1.3 MatrixGenerics_1.20.0
[11] matrixStats_1.5.0 here_1.0.1
[13] rio_1.2.3 data.table_1.17.0
[15] magrittr_2.0.3 lubridate_1.9.4
[17] forcats_1.0.0 stringr_1.5.1
[19] dplyr_1.1.4 purrr_1.0.4
[21] readr_2.1.5 tidyr_1.3.1
[23] tibble_3.2.1 ggplot2_3.5.1
[25] tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] gtable_0.3.6 xfun_0.51 lattice_0.22-7
[4] tzdb_0.4.0 vctrs_0.6.5 tools_4.5.1
[7] parallel_4.5.1 pkgconfig_2.0.3 R.oo_1.27.0
[10] Matrix_1.7-4 RcppParallel_5.1.10 lifecycle_1.0.4
[13] GenomeInfoDbData_1.2.14 compiler_4.5.1 munsell_0.5.1
[16] DESeq2_1.48.2 codetools_0.2-20 pillar_1.10.1
[19] crayon_1.5.3 R.utils_2.13.0 BiocParallel_1.42.2
[22] DelayedArray_0.34.1 abind_1.4-8 locfit_1.5-9.11
[25] zip_2.3.2 tidyselect_1.2.1 stringi_1.8.4
[28] rprojroot_2.0.4 grid_4.5.1 colorspace_2.1-1
[31] cli_3.6.4 SparseArray_1.8.1 S4Arrays_1.8.1
[34] utf8_1.2.4 withr_3.0.2 scales_1.3.0
[37] UCSC.utils_1.4.0 bit64_4.6.0-1 timechange_0.3.0
[40] XVector_0.48.0 httr_1.4.7 bit_4.5.0.1
[43] R.methodsS3_1.8.2 hms_1.1.3 stringfish_0.16.0
[46] openxlsx_4.2.8 evaluate_1.0.3 knitr_1.49
[49] rlang_1.1.5 Rcpp_1.0.14 glue_1.8.0
[52] pkgload_1.4.1 vroom_1.6.5 jsonlite_1.9.0
[55] R6_2.6.1
Describe the issue
I have WES data from cancer cell lines established from tumor samples in our lab. I dont have corresponding, matched normal healthy samples, but we have sequenced a number of non-malignant, cancer associated fibroblasts, that I use as process-matched normal samples. I realized that in the genes.tsv files there is not a single gene with a value of 0 in the C column (the minimum value is 0.42; all other values are > 0.5).
Also, in the "type" column of genes.tsv files, there is just one deletion, in all samples.
Related to that: I know from some cell lines based on aCGH data that they have a homozygous deletion in a gene, but when I look for that gene in genes.tsv file of the respective samples, this gene is missing in the genes.tsv file.
To Reproduce
This is the PureCN.R command:
Expected behavior
I would have expected more homozygous deletions, based on previous aCGH from replicates.
Log file
B-allele frequency plot
This is the screenshot from the plots on page 2 of the sample.pdf file of the most likely solution:
Session Info