Problem with PACBIO model

[ajsa-nukovic]

Hi everyone,

I tested DeepVariant 1.1.0 on PACBIO data (HG002, chr20), using PACBIO model and got following error in call_variants.py script:
I0519 14:15:32.843033 139847781275392 call_variants.py:338] Shape of input examples: [100, 221, 8] Traceback (most recent call last): File "/sbgenomics/workspaces/aa109ba2-5d46-4def-a398-4a7e1ee8806e/tasks/2555e949-a7d7-40a9-9a24-b002adf182c2/deepvariant-1-0-0/Bazel.runfiles_o79jsi96/runfiles/com_google_deepvariant/deepvariant/call_variants.py", line 502, in <module> tf.compat.v1.app.run() File "/usr/local/lib/python3.6/dist-packages/tensorflow/python/platform/app.py", line 40, in run _run(main=main, argv=argv, flags_parser=_parse_flags_tolerate_undef) File "/sbgenomics/workspaces/aa109ba2-5d46-4def-a398-4a7e1ee8806e/tasks/2555e949-a7d7-40a9-9a24-b002adf182c2/deepvariant-1-0-0/Bazel.runfiles_o79jsi96/runfiles/absl_py/absl/app.py", line 299, in run _run_main(main, args) File "/sbgenomics/workspaces/aa109ba2-5d46-4def-a398-4a7e1ee8806e/tasks/2555e949-a7d7-40a9-9a24-b002adf182c2/deepvariant-1-0-0/Bazel.runfiles_o79jsi96/runfiles/absl_py/absl/app.py", line 250, in _run_main sys.exit(main(argv)) File "/sbgenomics/workspaces/aa109ba2-5d46-4def-a398-4a7e1ee8806e/tasks/2555e949-a7d7-40a9-9a24-b002adf182c2/deepvariant-1-0-0/Bazel.runfiles_o79jsi96/runfiles/com_google_deepvariant/deepvariant/call_variants.py", line 492, in main use_tpu=FLAGS.use_tpu, File "/sbgenomics/workspaces/aa109ba2-5d46-4def-a398-4a7e1ee8806e/tasks/2555e949-a7d7-40a9-9a24-b002adf182c2/deepvariant-1-0-0/Bazel.runfiles_o79jsi96/runfiles/com_google_deepvariant/deepvariant/call_variants.py", line 351, in call_variants num_channels_in_checkpoint_model, example_shape[2])) **ValueError: The number of channels in examples and checkpoint should match, but the checkpoint has 9 channels while the examples have 8.**

My command line looks like this:

export HOME=/root && N_SHARDS=32 && LOGDIR=/opt/deepvariant/logs/ && mkdir -p "${LOGDIR}" && ( /usr/bin/time seq 0 $((N_SHARDS-1)) | parallel --eta --halt 2 --joblog "${LOGDIR}/log" --res "${LOGDIR}" python /opt/deepvariant/bin/make_examples.zip --mode calling --task {} --examples "./examples.tfrecord@${N_SHARDS}.gz" --alt_aligned_pileup diff_channels --reads /Projects/aa109ba2-5d46-4def-a398-4a7e1ee8806e/HG002.merged.bam --ref Projects/aa109ba2-5d46-4def-a398-4a7e1ee8806e/GRCh38ERCC.ensembl.fasta --norealign_reads --regions 20 --sample_name HG002 ) > ./make_examples.log 2>&1 && ( python /opt/deepvariant/bin/call_variants.zip --outfile ./call_variants_output.tfrecord.gz --examples ./examples.tfrecord@${N_SHARDS}.gz --checkpoint /opt/models/pacbio/model.ckpt --use_openvino --num_readers 32 ) > ./call_variants.log 2>&1 && ( python /opt/deepvariant/bin/postprocess_variants.zip --ref /Projects/aa109ba2-5d46-4def-a398-4a7e1ee8806e/GRCh38ERCC.ensembl.fasta --infile ./call_variants_output.tfrecord.gz --outfile ./HG002.vcf ) > ./postprocess_variants.log 2>&1

When I try and run with HYBRID model, everything goes smoothly. Do you have some input on this?

Thanks,
Ajsa

[pichuan]

Hi @ajsa-nukovic ,
Sorry for the confusion. Starting from v1.1.0, we added an additional channel to our PacBio model, and tried to simplify the flags in the one-step run_deepvariant by adding just one flag --use_hp_information, which you can set to false if you're BAM is not phased, and set to true if your BAM is phased.

Example:
https://github.com/google/deepvariant/blob/r1.1/docs/deepvariant-pacbio-model-case-study.md#run-deepvariant-on-haplotagged-chromosome-20-alignments

However, from your command, I see that you're running make_examples directly, which is a reasonable use. But you're specifying your own make_examples flag for 1.1.0 code and model, you'll need to add --add_hp_channel to make sure the last channel is added.

So, to summarize, if you're specifying your own flags for the make_examples step:

• If your BAM is NOT phased: use --sort_by_haplotypes=false --parse_sam_aux_fields=false --add_hp_channel=true.
• If your BAM is phased: use --sort_by_haplotypes=true --parse_sam_aux_fields=true --add_hp_channel=true.
  Basically, in v1.1.0, add_hp_channel needs to always be set to true when you're running make_examples for PacBio.

I'll see if I can update the r1.1 documentation to avoid future confusions.

[pichuan]

@ajsa-nukovic Hopefully adding the --add_hp_channel in make_examples work for you. If not, let me know.

[ajsa-nukovic]

Hi @pichuan ,

Yeah, it was definitely --add_hp_channel flag that made the difference. PACBIO model works great now.

Thanks a lot for your fast reply. :)