Reduced recall due to low MAPQ after surjecting long-read graph alignments to GRCh38

[matthewglasenapp]

Hello,

Cross-posting this issue after DeepVariant recommended following up with the vg team.

vg giraffe-mapped PacBio reads surjected to GRCh38 have very low MAPQ in the MHC region of chromosome 6, which leads to poor performance by DeepVariant. Do you have any tips for running vg giraffe + DeepVariant when working with long reads? My current workaround has been to reassign MAPQ based on the MAPQ scores from minimap2 alignments to the linear GRCh38 reference.

I've been using version v1.64.1 with the following command:

vg giraffe \
-b hifi \
-Z <reference_gbz> \
--ref-paths <ref_paths> \
-f <input_file> \
-p \
-P \
-o BAM \

[adamnovak]

We're still working on the preprint for good variant calling with long reads with Giraffe and DeepVariant; we have a HiFi model the Google team trained for us that works well that we are writing up.

But even there we haven't spot-checked MHC. In fact unless it's "high confidence" in the calling benchmarks we use, I don't think we would have been evaluating it at all.

It looks like you're likely to be using the HPRC v1.1 graph with GRCh38 as its backbone reference, with haplotype sampling?

My guess is that the MAPQ is low because there are still too many versions of MHC in the graph after haplotype sampling, or that the MHCs in the base graph are under-aligned somehow. Are your alignment scores also low, or are these well-aligned reads that aren't excessively softclipping but still get low MAPQ?

@glennhickey Are we confident in the assembly alignment there in v1.1? I just checked my hprc-v1.1-mc-grch38.gbz and it doesn't have the GRCh38 alt loci in there, which is one way MHC could get weird.

It could just be that the graph doesn't flatten all of MHC down onto the linear reference, and the graph splits apart for appreciable distances so the different versions get treated as distinct mapping locations, and then the lower MAPQ accurately reflects uncertainty about which version of MHC the read actually came from. But then vg surject would be expected to have trouble squishing things down onto one linear reference for you to see the low MAPQs there. (Do you also have a coverage deficit in your GRCh38 linear MHC region?)

Just in general I'd recommend trying the CHM13-based graphs, since CHM13 is a better reference than GRCh38 and using it as the backbone over GRCh38 produces better graphs. But I'm not sure that will really help.

[matthewglasenapp]

Hi Adam,

Thanks for the detailed response!

I believe I've been using the full HPRC v1.0 graph with GRCh38 as the backbone, with no haplotype sampling (hprc-v1.0-mc-grch38-minaf.0.1.gbz).

Although the reads had low MAPQ, their alignment scores look good.

For reads mapping to HLA-A after surjection (chr6:29941259–29949572), I found:
Mean Read Length: 3841
Mean MAPQ: 22
Mean AS: 3427

For comparison, the same reads aligned to linear GRCh38 with minimap2 give:
Mean Read Length: 3841
Mean MAPQ=59.7
Mean AS=3392

This seems consistent with your hypothesis that the low MAPQ might be due to redundancy from many similar haplotypes in the graph, rather than poor alignments. Is there a particular graph build you’d recommend for this kind of work? I am open to switching to CHM13 as the backbone.

[matthewglasenapp]

Update: I still observed the reduced MapQ despite normal AS when using the HPRC v2.0 graph hprc-v2.0-mc-grch38.gbz with haplotype sampling via vg giraffe.

vg giraffe \
-b hifi \
-Z hprc-v2.0-mc-grch38.gbz \
--haplotype-name hprc-v2.0-mc-grch38.hapl \
--kff-name HG002_kmc.kff \
--set-reference GRCh38 \
-f HG002.fastq.gz \
-p \
-P \
-o BAM \
--ref-paths hprc-v2.0-mc-grch38.dict 
-R <RG_ID> 
-N HG002 | samtools sort \
-o {output_file}

Read count: 1743
Mean length: 3841
Mean AS: 3427
Average MAPQ: 22

Here is the vg giraffe output:

Sampling haplotypes
Loading GBZ from hprc-v2.0-mc-grch38.gbz
Updating reference samples
Loading haplotype information from hprc-v2.0-mc-grch38.hapl
Assigning reference paths to GBWT construction jobs
Assigned 195 reference paths (0 unassigned) in 0.394079 seconds
Reading kmer counts
Read the kmer counts in 71.2694 seconds
Estimating kmer coverage
Estimated kmer coverage in 3.22431 seconds
Building GBWT
Running 13 GBWT construction jobs in parallel
Total: 2 haplotypes for 195 chains (148 full, 256622 subchains, 785 fragments); connection rate 0.0140557; included 195 reference paths; used 303642206 kmers with average score 0.222595
Finished the jobs in 447.32 seconds
Merging the partial indexes
Merged the indexes in 20.5 seconds
Built the GBWT in 467.821 seconds
Building GBWTGraph
Saving GBWT and GBWTGraph to hprc-v2.0-mc-grch38.HG002.gbz
Preparing Indexes
[IndexRegistry]: Constructing distance index for Giraffe.
[IndexRegistry]: Constructing minimizer index and associated zipcodes.
	use parameters -k 31 -w 50 -W 
[IndexRegistry]: Finding frequent kmers using the fast algorithm
[IndexRegistry]: Found 78672 kmers with more than 500 hits
Loading Minimizer Index
Loading Zipcodes
Loading GBZ
Loading Distance Index v2
Paging in Distance Index v2
Applying overlay
Initializing MinimizerMapper
Loading and initialization: 2196.15 seconds
Of which Distance Index v2 paging: 4.7011 seconds
Mapping reads to "-" (BAM)
--watchdog-timeout 30
--batch-size 10
--match 1
--mismatch 4
--gap-open 6
--gap-extend 1
--full-l-bonus 5
--max-multimaps 1
--hit-cap 0
--hard-hit-cap 13614
--score-fraction 1
--max-min 79
--min-coverage-flank 250
--num-bp-per-min 152
--downsample-window-length 227
--downsample-window-count 15
--distance-limit 200
--max-extensions 800
--max-alignments 3
--cluster-score 50
--pad-cluster-score 20
--cluster-coverage 0.3
--max-extension-mismatches 4
--extension-score 1
--extension-set 20
--rescue-attempts 15
--max-fragment-length 2000
--paired-distance-limit 2
--rescue-subgraph-size 4
--rescue-seed-limit 100
--mapq-score-window 0
--mapq-score-scale 0.001
--align-from-chains
--zipcode-tree-score-threshold 100
--pad-zipcode-tree-score-threshold 50
--zipcode-tree-coverage-threshold 0.5
--zipcode-tree-scale 2
--min-to-fragment 2
--max-to-fragment 15
--max-direct-chain 0
--gapless-extension-limit 0
--fragment-max-graph-lookback-bases 500
--fragment-max-graph-lookback-bases-per-base 0.025
--fragment-max-read-lookback-bases 18446744073709551615
--fragment-max-read-lookback-bases-per-base 1
--max-fragments 15000
--fragment-max-indel-bases 15000
--fragment-max-indel-bases-per-base 0.1
--fragment-gap-scale 1.44952
--fragment-points-per-possible-match 0
--fragment-score-fraction 0
--fragment-max-min-score 50000
--fragment-min-score 2
--fragment-set-score-threshold 70
--min-chaining-problems 6
--max-chaining-problems 2147483647
--max-graph-lookback-bases 20000
--max-graph-lookback-bases-per-base 0.10501
--max-read-lookback-bases 18446744073709551615
--max-read-lookback-bases-per-base 1
--max-indel-bases 5000
--max-indel-bases-per-base 2.45
--item-bonus 20
--item-scale 1
--gap-scale 0.2
--points-per-possible-match 0
--chain-score-threshold 200
--min-chains 2
--min-chain-score-per-base 0.1
--max-min-chain-score 1100
--max-skipped-bases 1000
--max-chains-per-tree 3
--max-chain-connection 233
--max-tail-length 68
--max-dp-cells 8000000000
--max-tail-gap 150
--max-middle-gap 500
--max-tail-dp-length 30000
--max-middle-dp-length 2147483647
--wfa-max-mismatches 2
--wfa-max-mismatches-per-base 0.05
--wfa-max-max-mismatches 15
--wfa-distance 33
--wfa-distance-per-base 0.195722
--wfa-max-distance 240
--min-unique-node-fraction 0
--prune-low-cplx
T10:	warning[MinimizerMapper::map_from_chains]: No seeds found for m84039_240622_113450_s1/23725535/ccs!
warning[MinimizerMapper::align_sequence_between]: Trimmed back tips 41 times on graph between 95136553+132 as 197+ and 95136718+12 as 1+ leaving 58 nodes and 2 tips for read m84039_240622_113450_s1/102436949/ccs
Mapped 417532 reads across 13 threads in 116.749 seconds with 0 additional single-threaded seconds.
Mapping speed: 275.101 reads per second per thread
Used 1613.45 CPU-seconds (including output).](url)
Achieved 258.782 reads per CPU-second (including output)
Used 12320383559993 CPU instructions (not including output).
Mapping slowness: 29.5076 M instructions per read at 7636.04 M mapping instructions per inclusive CPU-second
Memory footprint: 88.443 GB