PEPPER-Margin-DeepVariant is a haplotype-aware variant calling pipeline for processing third-generation nanopore sequence data. It outperforms the short-read-based single-nucleotide-variant identification method at the whole-genome scale and produces high-quality single-nucleotide variants in segmental duplications and low-mappability regions where short-read-based genotyping fails. It also can provide highly contiguous phase blocks across the genome with nanopore reads and perform de novo assembly polishing to produce diploid assemblies with high accuracy. This pipeline is also applicable to PacBio HiFi data, providing an efficient solution with superior performance over the current WhatsHap-DeepVariant standard.
Allocate an interactive session and run the program. Sample session:
[user@biowulf]$ sinteractive --mem=8g -c64 --gres=gpu:v100x:1,lscratch:20 [user@cn3316 ~]$ module load PEPPER_deepvariant [+] Loading singularity 3.8.5-1 on cn4235 [+] Loading pepper_deepvariant 0.7 ...The basic usage command for PEPPER_deepvariant is:
[user@cn3123 user]$ run_pepper_margin_deepvariant -h
usage: run_pepper_margin_deepvariant [-h] [--version] {call_variant} ...
Run PEPPER-Margin-DeepVariant for variant calling.Example run: run_pepper_margin_deepvariant -b -f -o -t <--ont_r9_guppy5_sup/--ont_r10_q20/--hifi>
positional arguments:
{call_variant}
optional arguments:
-h, --help show this help message and exit
--version Show version.
[user@cn3123 user]$ run_pepper_margin_deepvariant call_variant -h
usage: run_pepper_margin_deepvariant call_variant -b BAM -f FASTA -o
OUTPUT_DIR -t THREADS
[-r REGION] [-g] [--gvcf]
[-s SAMPLE_NAME]
[--only_pepper]
[--only_haplotag]
[--phased_output]
[--skip_final_phased_bam]
[-p OUTPUT_PREFIX] [-k]
[--dry]
[--pepper_model PEPPER_MODEL]
[--pepper_quantized]
[--no_pepper_quantized]
[--pepper_downsample_rate PEPPER_DOWNSAMPLE_RATE]
[--pepper_region_size PEPPER_REGION_SIZE]
[--pepper_include_supplementary]
[--pepper_min_mapq PEPPER_MIN_MAPQ]
[--pepper_min_snp_baseq PEPPER_MIN_SNP_BASEQ]
[--pepper_min_indel_baseq PEPPER_MIN_INDEL_BASEQ]
[--pepper_snp_frequency PEPPER_SNP_FREQUENCY]
[--pepper_insert_frequency PEPPER_INSERT_FREQUENCY]
[--pepper_delete_frequency PEPPER_DELETE_FREQUENCY]
[--pepper_min_coverage_threshold PEPPER_MIN_COVERAGE_THRESHOLD]
[--pepper_candidate_support_threshold PEPPER_CANDIDATE_SUPPORT_THRESHOLD]
[--pepper_snp_candidate_frequency_threshold PEPPER_SNP_CANDIDATE_FREQUENCY_THRESHOLD]
[--pepper_indel_candidate_frequency_threshold PEPPER_INDEL_CANDIDATE_FREQUENCY_THRESHOLD]
[--pepper_skip_indels]
[--pepper_allowed_multiallelics PEPPER_ALLOWED_MULTIALLELICS]
[--pepper_snp_p_value PEPPER_SNP_P_VALUE]
[--pepper_insert_p_value PEPPER_INSERT_P_VALUE]
[--pepper_delete_p_value PEPPER_DELETE_P_VALUE]
[--pepper_snp_p_value_in_lc PEPPER_SNP_P_VALUE_IN_LC]
[--pepper_insert_p_value_in_lc PEPPER_INSERT_P_VALUE_IN_LC]
[--pepper_delete_p_value_in_lc PEPPER_DELETE_P_VALUE_IN_LC]
[--pepper_snp_q_cutoff PEPPER_SNP_Q_CUTOFF]
[--pepper_indel_q_cutoff PEPPER_INDEL_Q_CUTOFF]
[--pepper_snp_q_cutoff_in_lc PEPPER_SNP_Q_CUTOFF_IN_LC]
[--pepper_indel_q_cutoff_in_lc PEPPER_INDEL_Q_CUTOFF_IN_LC]
[--pepper_report_snp_above_freq PEPPER_REPORT_SNP_ABOVE_FREQ]
[--pepper_report_indel_above_freq PEPPER_REPORT_INDEL_ABOVE_FREQ]
[--margin_haplotag_model MARGIN_HAPLOTAG_MODEL]
[--margin_phase_model MARGIN_PHASE_MODEL]
[--dv_model DV_MODEL]
[--dv_model_snp DV_MODEL_SNP]
[--dv_model_indel DV_MODEL_INDEL]
[--dv_alt_aligned_pileup DV_ALT_ALIGNED_PILEUP]
[--dv_alt_aligned_pileup_snp DV_ALT_ALIGNED_PILEUP_SNP]
[--dv_alt_aligned_pileup_indel DV_ALT_ALIGNED_PILEUP_INDEL]
[--dv_pileup_image_width DV_PILEUP_IMAGE_WIDTH]
[--dv_realign_reads DV_REALIGN_READS]
[--dv_partition_size DV_PARTITION_SIZE]
[--dv_min_mapping_quality DV_MIN_MAPPING_QUALITY]
[--dv_min_base_quality DV_MIN_BASE_QUALITY]
[--dv_vsc_min_fraction_indels DV_VSC_MIN_FRACTION_INDELS]
[--dv_vsc_min_fraction_snps DV_VSC_MIN_FRACTION_SNPS]
[--dv_sort_by_haplotypes DV_SORT_BY_HAPLOTYPES]
[--dv_parse_sam_aux_fields DV_PARSE_SAM_AUX_FIELDS]
[--dv_add_hp_channel DV_ADD_HP_CHANNEL]
[--dv_use_hp_information DV_USE_HP_INFORMATION]
[--dv_use_multiallelic_mode DV_USE_MULTIALLELIC_MODE]
[--ont_r9_guppy5_sup]
[--ont_r10_q20] [--hifi]
[-h]
...
user@cn3123 user]$ cp $PEPPER_DATA/* .
user@cn3123 user]$ mkdir output
[user@cn3123 user]$ run_pepper_margin_deepvariant call_variant \
-b HG003_chr20_demo.bam \
-f GRCh38_no_alt.chr20.fa \
-o output \
-t 64 \
--hifi \
-g &
...
[04-19-2024 18:25:25] INFO: HiFi VARIANT CALLING MODE SELECTED.
[04-19-2024 18:25:25] INFO: MODE: PEPPER
[04-19-2024 18:25:25] INFO: THRESHOLDS ARE SET TO:
[04-19-2024 18:25:25] INFO: MIN MAPQ: 5
[04-19-2024 18:25:25] INFO: MIN SNP BASEQ: 10
[04-19-2024 18:25:25] INFO: MIN INDEL BASEQ: 10
[04-19-2024 18:25:25] INFO: MIN SNP FREQUENCY: 0.1
[04-19-2024 18:25:25] INFO: MIN INSERT FREQUENCY: 0.12
[04-19-2024 18:25:25] INFO: MIN DELETE FREQUENCY: 0.1
[04-19-2024 18:25:25] INFO: MIN COVERAGE THRESHOLD: 2
[04-19-2024 18:25:25] INFO: MIN CANDIDATE SUPPORT: 2
[04-19-2024 18:25:25] INFO: MIN SNP CANDIDATE FREQUENCY: 0.1
[04-19-2024 18:25:25] INFO: MIN INDEL CANDIDATE FREQUENCY: 0.1
[04-19-2024 18:25:25] INFO: SKIP INDEL CANDIDATES: False
[04-19-2024 18:25:25] INFO: MAX ALLOWED CANDIDATE IN ONE SITE: 4
[04-19-2024 18:25:25] INFO: MIN SNP PREDICTIVE VALUE: 0
[04-19-2024 18:25:25] INFO: MIN INSERT PREDICTIVE VALUE: 0
[04-19-2024 18:25:25] INFO: MIN DELETE PREDICTIVE VALUE: 0
[04-19-2024 18:25:25] INFO: SNP QV CUTOFF FOR RE-GENOTYPING: 15
[04-19-2024 18:25:25] INFO: INDEL QV CUTOFF FOR RE-GENOTYPING: 20
[04-19-2024 18:25:25] INFO: REPORT ALL SNPs ABOVE THRESHOLD: 0
[04-19-2024 18:25:25] INFO: REPORT ALL INDELs ABOVE THRESHOLD: 0
[04-19-2024 18:25:25] INFO: LOW COMPLEXITY REGION SETUP:
[04-19-2024 18:25:25] INFO: MIN SNP PREDICTIVE VALUE: 0
[04-19-2024 18:25:25] INFO: MIN INSERT PREDICTIVE VALUE: 0
[04-19-2024 18:25:25] INFO: MIN DELETE PREDICTIVE VALUE: 0
[04-19-2024 18:25:25] INFO: SNP QV CUTOFF FOR RE-GENOTYPING: 15
[04-19-2024 18:25:25] INFO: INDEL QV CUTOFF FOR RE-GENOTYPING: 20
[04-19-2024 18:25:25] INFO: CALL VARIANT MODULE SELECTED
[04-19-2024 18:25:26] INFO: RUN-ID: 04192024_182526
[04-19-2024 18:25:26] INFO: IMAGE OUTPUT: /data/denisovga/PEPPER_deepvariant_hifi/output/pepper/images_04192024_182526/
[04-19-2024 18:25:26] INFO: STEP 1/3 GENERATING IMAGES:
[04-19-2024 18:25:26] INFO: COMMON CONTIGS FOUND: ['chr20']
[04-19-2024 18:25:26] INFO: TOTAL CONTIGS: 1 TOTAL INTERVALS: 645 TOTAL BASES: 64444166
[04-19-2024 18:25:26] INFO: STARTING PROCESS: 0 FOR 21 INTERVALS
[04-19-2024 18:25:26] INFO: [THREAD 00] 10/21 COMPLETE (47%) [ELAPSED TIME: 0 Min 0 Sec]
[04-19-2024 18:25:26] INFO: [THREAD 00] 20/21 COMPLETE (95%) [ELAPSED TIME: 0 Min 0 Sec]
[04-19-2024 18:25:26] INFO: THREAD 0 FINISHED SUCCESSFULLY.
[04-19-2024 18:25:26] INFO: FINISHED IMAGE GENERATION
[04-19-2024 18:25:26] INFO: TOTAL ELAPSED TIME FOR GENERATING IMAGES: 0 Min 0 Sec
[04-19-2024 18:25:26] INFO: STEP 2/3 RUNNING INFERENCE
[04-19-2024 18:25:26] INFO: OUTPUT: /data/denisovga/PEPPER_deepvariant_hifi/output/pepper/predictions_04192024_182526/
[04-19-2024 18:25:26] INFO: TOTAL GPU AVAILABLE: 1
[04-19-2024 18:25:26] INFO: AVAILABLE GPU DEVICES: [0]
[04-19-2024 18:25:26] INFO: TOTAL CALLERS: 4
[04-19-2024 18:25:26] INFO: TOTAL THREADS PER CALLER: 8
[04-19-2024 18:25:36] INFO: TOTAL FILES: 64.
[04-19-2024 18:25:36] INFO: FILES COMPLETED: 1/64.
[04-19-2024 18:25:36] INFO: FILES COMPLETED: 2/64.
[04-19-2024 18:25:36] INFO: FILES COMPLETED: 3/64.
[04-19-2024 18:25:36] INFO: FILES COMPLETED: 4/64.
...
[user@cn3123 user]$ nvidia-smi
+---------------------------------------------------------------------------------------+
| NVIDIA-SMI 535.129.03 Driver Version: 535.129.03 CUDA Version: 12.2 |
|-----------------------------------------+----------------------+----------------------+
| GPU Name Persistence-M | Bus-Id Disp.A | Volatile Uncorr. ECC |
| Fan Temp Perf Pwr:Usage/Cap | Memory-Usage | GPU-Util Compute M. |
| | | MIG M. |
|=========================================+======================+======================|
| 0 NVIDIA A100-SXM4-80GB On | 00000000:C7:00.0 Off | 0 |
| N/A 33C P0 69W / 400W | 1662MiB / 81920MiB | 5% Default |
| | | Disabled |
+-----------------------------------------+----------------------+----------------------+
+---------------------------------------------------------------------------------------+
| Processes: |
| GPU GI CI PID Type Process name GPU Memory |
| ID ID Usage |
|=======================================================================================|
| 0 N/A N/A 2023948 C /usr/bin/python3 1648MiB |
+---------------------------------------------------------------------------------------+
...
[04-19-2024 18:26:32] INFO: FILES COMPLETED: 31/64.
[04-19-2024 18:26:32] INFO: FILES COMPLETED: 32/64.
[04-19-2024 18:26:32] INFO: PREDICTION GENERATED SUCCESSFULLY.
[04-19-2024 18:26:32] ELAPSED TIME: 0 Min 9 Sec
[04-19-2024 18:26:32] INFO: STEP 3/3 FINDING CANDIDATES
[04-19-2024 18:26:32] INFO: OUTPUT: /data/denisovga/PEPPER_deepvariant_hifi/output/pepper/
[04-19-2024 18:26:32] INFO: STARTING CANDIDATE FINDING.
[04-19-2024 18:26:33] INFO: FINISHED PROCESSING, TOTAL CANDIDATES FOUND: 659
[04-19-2024 18:26:33] INFO: FINISHED PROCESSING, TOTAL VARIANTS IN PEPPER: 312
[04-19-2024 18:26:33] INFO: FINISHED PROCESSING, TOTAL VARIANTS SELECTED FOR RE-GENOTYPING: 347
[04-19-2024 18:26:33] INFO: FINISHED PROCESSING, TOTAL SNP VARIANTS SELECTED FOR RE-GENOTYPING: 276
[04-19-2024 18:26:33] INFO: FINISHED PROCESSING, TOTAL INDEL VARIANTS SELECTED FOR RE-GENOTYPING: 71
[04-19-2024 18:26:33] INFO: TOTAL TIME SPENT ON CANDIDATE FINDING: 0 Min 1 Sec
[04-19-2024 18:26:33] INFO: TOTAL ELAPSED TIME FOR FINDING CANDIDATES: 0 Min 12 Sec
...
Running OpenMP with 8 threads.
> Parsing model parameters from file: /opt/margin_dir/params/phase/allParams.haplotag.pb-hifi.snp.json
> Parsed 659 total VCF entries from output/intermediate_files/PEPPER_VARIANT_FULL.vcf.gz; kept 422 HETs, skipped 0 for region, 121 for not being PASS, 116 for being homozygous, 0 for being INDEL
> Set up bam chunker in 0s with chunk size 100000 and overlap 10000 (for region=chr20), resulting in 3 total chunks
> Ordering chunks by estimated depth
> Setup complete, beginning run
> Polishing 66% complete (2/3, 0s). Estimated time remaining: 0s
> Starting merge
Expected three tokens in header line, got 2
This usually means you have multiple primary alignments with the same read ID.
You can identify whether this is the case with this command:
samtools view -F 0x904 YOUR.bam | cut -f 1 | sort | uniq -c | awk '$1 > 1'
Expected three tokens in header line, got 2
This usually means you have multiple primary alignments with the same read ID.
You can identify whether this is the case with this command:
samtools view -F 0x904 YOUR.bam | cut -f 1 | sort | uniq -c | awk '$1 > 1'
real 0m44.465s
user 1m26.239s
sys 0m0.131s
...
[user@cn3316 ~]$ exit salloc.exe: Relinquishing job allocation 46116226 [user@biowulf ~]$
Create a batch input file (e.g. fithichip.sh). For example:
#!/bin/bash module load PEPPER_deepvariant
Submit this job using the Slurm sbatch command.
sbatch [--cpus-per-task=#] [--mem=#] fithichip.sh