Hello everyone, I am currently using the Docker image for PureCN and this is my workflow. Could anyone advise me on whether this workflow is right? The reason I ask this is because the results generated using this workflow are significantly different from the results I get from using Sequenza for 10 paired tumor-germline samples:

#Step 1 - Generating interval files

Rscript $PURECN/IntervalFile.R --in-file $BAIT_COORDINATES \
    --fasta $REF_GENOME_b37 --out-file $OUT_REF/baits_hg19_intervals.txt \
    --off-target --genome hg19 \
    --export $OUT_REF/baits_optimized_hg19.bed  \
    --mappability $MAPPABILITY_FILE_B37

#Step 2 - Run mutect on 10 germline samples, create genomics db and use CreateSomaticPanelOfNormals

> $PURECN_G/G_sample_map

        count=1
        while IFS='' read -r line
        do
        if [ ! -z $GERMLINE_PFX ] && [ ! -f $PURECN_G/${GERMLINE_PFX}.vcf.gz ] && [ $count -le 10 ]; then
                count=$((count+1))
                echo $GERMLINE_PFX
                #For each germline bam files
                echo "Running Mutect2 on germline bams.................."

                #Run mutect on 10 germline samples
                gatk Mutect2 \
                -R $REF_GENOME_b37 \
                -I $SAVEROOT/$TUMOR_PFX/recal/bqsr_${GERMLINE_PFX}*.bam \
                --max-mnp-distance 0 \
                --genotype-germline-sites true \
                --genotype-pon-sites true \
                --interval-padding 75 \
                --L $BEDFILE_0B \
                -O $PURECN_G/${GERMLINE_PFX}.vcf.gz

                #Create a map file of all germline vcfs
                echo -e "$GERMLINE_PFX\t$PURECN_G/${GERMLINE_PFX}.vcf.gz" >> $PURECN_G/G_sample_map
        fi

       done<$LIST_OF_SAMPLES

        #Create a Genomics DB from the normal Mutect2 calls 
        gatk GenomicsDBImport \
        -R $REF_GENOME_b37 \
        -L $BEDFILE_0B \
        --genomicsdb-workspace-path $PURECN_G/pon_db \
        --sample-name-map $PURECN_G/G_sample_map

        #Combine the normal calls using CreateSomaticPanelOfNormals
        gatk CreateSomaticPanelOfNormals \
        -R $REF_GENOME_b37 \
        --germline-resource $GERMLINE_RESOURCE \
        -V gendb://$PURECN_G/pon_db \
        -O $PURECN_G/pon.vcf.gz

#Step 3 - Run Mutect on unmatched tumor samples

while IFS='' read -r line
        do
                if [ ! -z $TUMOR_PFX ]  && [ ! -f $PURECN_T/${TUMOR_PFX}.vcf.gz ]; then
                        #Run Mutect2 for each tumor sample
                        gatk Mutect2 \
                        -R $REF_GENOME_b37 \
                        -I $SAVEROOT/$TUMOR_PFX/recal/bqsr_${TUMOR_PFX}*.bam \
                        -pon $PURECN_G/pon.vcf.gz \
                        --germline-resource $GERMLINE_RESOURCE\
                        --genotype-germline-sites true \
                        --genotype-pon-sites true \
                        --interval-padding 75 \
                        -L $BEDFILE_0B \
                        -O $PURECN_T/${TUMOR_PFX}.vcf.gz

                fi

       done<$LIST_OF_SAMPLES

#Step 4 - calculate gc normalized coverages for germlines

count=1
            while IFS='' read -r line
            do
            TUMOR_PFX=$(echo $line| awk -F'|' '{print $1}')
            GERMLINE_PFX=$(echo $line| awk -F'|' '{print $2}')
            if [ ! -z $GERMLINE_PFX ] && [ $count -le 10 ] && [ ! -f $PURECN_G/${GERMLINE_PFX}_coverage_loess.txt.gz ] ; then
                    count=$((count+1))
                    echo $GERMLINE_PFX

                    Rscript $PURECN/Coverage.R \
                    --out-dir $PURECN_G \
                    --bam $SAVEROOT/$TUMOR_PFX/recal/bqsr_${GERMLINE_PFX}*.bam \
                    --intervals $OUT_REF/baits_hg19_intervals.txt
            fi

           done<$LIST_OF_SAMPLES

#Step 5 - calculate gc normalized coverages for tumors

while IFS='' read -r line
        do
        TUMOR_PFX=$(echo $line| awk -F'|' '{print $1}')
        if [ ! -z $TUMOR_PFX ] && [ ! -f $PURECN_T/${TUMOR_PFX}_coverage_loess.txt.gz ]; then
                echo $TUMOR_PFX
                Rscript $PURECN/Coverage.R \
                --out-dir $PURECN_T \
                --bam $SAVEROOT/$TUMOR_PFX/recal/bqsr_${TUMOR_PFX}*.bam \
                --intervals $OUT_REF/baits_hg19_intervals.txt
        fi

        done<$LIST_OF_SAMPLES

#Step 6 - create normal db

ls -a $PURECN_G/*_loess.txt.gz | cat > $PURECN_G/normal_coverages.list
        ulimit -n 5000
        #For a Mutect2/GATK4 normal panel GenomicsDB (beta)
        Rscript $PURECN/NormalDB.R --out-dir $PURECN_G \
        --coverage-files $PURECN_G/normal_coverages.list \
        --normal-panel $PURECN_G/pon_db \
        --genome hg19

#Step 7 - calculate tumor cellularity and ploidy

while IFS='' read -r line
        do
                TUMOR_PFX=$(echo $line| awk -F'|' '{print $1}')

        if [ ! -z $TUMOR_PFX ] && ! grep $TUMOR_PFX $PURECN_RESULTS/purecn_unpaired_tumor_cellularity_ploidy.txt; then
                echo $TUMOR_PFX
                #For each tumor bam files
                # Production pipeline run (dont need stats file if Mutect2 was run)
                Rscript $PURECN/PureCN.R --out $PURECN_T \
                --tumor $PURECN_T/bqsr_${TUMOR_PFX}_T_coverage_loess.txt.gz \
                --sampleid $TUMOR_PFX \
                --vcf $PURECN_T/${TUMOR_PFX}.vcf.gz \
                --fun-segmentation PSCBS \
                --normaldb $PURECN_G/normalDB_hg19.rds \
                --mapping-bias-file $PURECN_G/mapping_bias_hg19.rds \
                --intervals $OUT_REF/baits_hg19_intervals.txt \
                --genome hg19 \
                --model betabin \
                --force --post-optimize --seed 123

My assay is targeted capture with an off target read of 30% and but if I specify --off-target in Step 1, in Step 6, I see an error:

INFO [2022-07-12 03:14:00] Removing 15793 intervals with low coverage in normalDB.
FATAL [2022-07-12 03:14:00] No intervals left after coverage checks.
FATAL [2022-07-12 03:14:00]
FATAL [2022-07-12 03:14:00] This is most likely a user error due to invalid input data or
FATAL [2022-07-12 03:14:00] parameters (PureCN 2.2.0).
Error: No intervals left after coverage checks.
This is most likely a user error due to invalid input data or
parameters (PureCN 2.2.0).
Execution halted

If I exclude --off-target in Step 1, the pipeline runs but produces different cellularity and ploidy estimates from that produced by Sequenza.

Hoping to get some idea on whether I am following the right steps for my pipeline, thank you!