deepvariant: Deepvariant fails to recognize PacBio BAM index file

Hi @Npaffen,

To expand a bit on my previous explanation, regarding the way you can view a variant coming from a neural network is through the idea of preserving $information$ $propagation$. You saw the previous description of how the different channels get encoded, but let’s start with a simpler version. Let say you have a one-line matrix with the following columns:

This could denote a simplified version of your read base representation, where you notice the middle denotes the variant, and the last column the channel it represents. Now I want to create new types of data from this, which will help me with identifying unique areas of patterns within it. This could be of the form of transformation of the values to ranges that are easier to detect differences among columns. One of these can be dividing all the values by 10, and labeling that channel 2. In this transformation, 10 represents a kernel I described previously and the output is a new feature map (a transformed matrix that helps with detecting unique features based on the numerical representation). Now the data would look like this:

Now imagine I create different transformations of these rows, to expand on specific areas among these values where intriguing patterns might emerge. Suppose I create 5 different transformations having then 5 channels with multiple copies of each row, in order to have a fuller dataset that mimics the number of reads. This data is multi-dimensional as it contains different values of X. This can be pretty hard to interpret, but we can collapse these differences to a 2D representation using t-SNE plots to visualize these differences. For example, if I do that to this dataset, I get the following plot:

This helps me understand which channels (feature maps) might be similar, or different. Now if we look at the heatmap for the original data it would look something like this (notice the variant in the middle):

If I generate a t-SNE plot of this, I see something like the following using colors to denote the different channels (feature maps):

Now comes the fun part! Let’s have different matrices (like the ones that generated the new channels above) that will identify interesting features that might expose one type of a genotype or another with confidence. Imagine that instead of dividing by 10, I find the values (matrices) that best helps separate the data in a BAM file I know should have specific variants at specific loci. I want to maximize that precision to be able to recall. Now you notice that the original data and transformations (feature maps) are linked preserving this $propagation$ $of$ $information$, as this flow of information is enabled through these sets of transformations. An interesting thing then begins to emerge as you move up the layers of transformation. For example, early on in the neural network’s set of transformations you will see patterns like this:

You might notice an explosion of features, with no specific patterns. These early steps are to generate a large variety of features to be able to have selection power for the later layers to use as input, for helping with the separation into distinct patterns for mapping to the different classes of genotypes confidently. For example, you can see distinct patterns forming as it reaches the later stages:

Once the pattern has been achieved like the following, then one can proceed with testing each genotype’s representation of the variant:

We want to see for which genotype the set of patterns (the feature map above) maximizes for, which will indicate the genotype present with a specific maximal probability.

First we test for $homozygous$ $reference$:

Next we test for $heterozygous$:

Finally we test for $homozygous$ $alternate$:

Now we can see there is a significant correlation with a heterozygous variant call.

So to call a variant site’s genotype you now have the power to traverse the whole neural network – through a set of transformations (you can tweak) to identify hidden patterns in the data – as it is all linked from the start (read encoding) to finish (genotype identification). This ability is enabled through the preservation of the propagation of information (which is deeply and directly linked) in order to help one optimize upon. That’s why you don’t want to throw anything away, as @AndrewCarroll and @pichuan mentioned as it helps you inspect and tweak the selection power and criteria as it is all interconnected.

Hope it helps, Paul

Yeah thanks too all here in this thread for the productive discussion. This helped me a lot to understand how deepvariant works! Good luck with the project!

Thank you @AndrewCarroll and @pichuan for the clarification. The calibration makes sense, and could be intriguing for inspecting DNN-resiliency. Having the same underlying Inception V3 network architecture for both PacBio and WGS, a point of natural comparability would be the logits kernel across all three genotypes:

Given visual similarity, these were confirmed via Euclidean distance (0.9931127, 0.8543731 and 1.052052, respectively). This indicates the feature set might exhibit strong similarity for interpretation.

Looking at one network (PacBio), it might be possible to confirm calibration by testing for network-resiliency. Via perturbation analysis it should be possible to get insight into a channel’s response under perturbation, and their binary interactions under such conditions. Keeping the variant unchanged within a window on each side for preserving the call, the inspection each channel vulnerability response to perturbation can be tested. This resulted in the following perturbation response ($c\_*$ denotes a channel, and $i\_*\_*$ represents a binary interaction between two channels):

The above can be mapped into a network of interactions among the channels:

Based on the above mapping, by testing well-interacting channels through a probabilistically value-update – within DeepVariant-acceptable values – it might be possible to check for shifts in genotype mimicking Mendelian violation. Selecting base_quality and staying within DeepVariant’s minimum acceptable value, random sampling with replacement was performed in the window outside the variant region. A shift in genotype was achieved giving a measure of network resiliency. Other channels being more strongly-connected could provide more aggressive shifts in genotype. This can be offset by restrictions in convolutional motifs within the network, or be shifted to the code via threshold limitations in make_examples.

Thank you and Happy 4th of July! Paul

Hi @Npaffen

Our typical recommendation for this would be to use DeepVariant WGS model to call variants with gVCF output in the Illumina samples, DeepVariant PACBIO model to call variants with gVCF output in the PacBio sample, and then use glnexus to jointly genotype the three samples together.

The variant call probabilities are well calibrated in each model individually, and so glnexus can operate on them.

Hi @pgrosu , The hybrid model you mentioned is combining PacBio and Illumina reads of the same individual. It is a DeepVariant model, not a DeepTrio model.

With the hybrid model, you can achieve better accuracy. This model was first developed as part of our submission to https://precision.fda.gov/challenges/10/results.

However, this is not a DeepTrio model and does not fit the use case that @Npaffen described unfortunately.

@AndrewCarroll and @pgrosu thanks for your view on things and @pgrosu especially for your fantastic explanation of the model!

From a data scientist perspective this makes totally sense. I fully agree that one should always gather as much data as possible since one can always remove data that seems to be not useful afterwards. My intention was not to say that one do not need those data or that you should not use the information which those HomRef sites emerged from. I talked to some of my colleagues yesterday and they agreed that they would expect a variant caller to produce variant calls and information/statistics related to those and all further information would be opt-in since they would have no need for this. But my team could just be the minority here and simply an opt-out option for this behavior would be highly recommended!

Hi @Npaffen @pgrosu

My sentence is written inexactly. I was thinking of the last word in the hyphen - reference in non-reference, and trying to clarify the reference calls. But this was more confusing rather than less.

Reference calls - 0/0 or ./. Variant calls - 0/1, 1/1, 1/2, 0/2, etc.

Hi @Npaffen

Yes, in phase variants, the reads are assigned a haplotag value. Briefly, in this process, a set of potential variants are scored with heuristics (no neural network) on the likelihood that they are heterozygous variants. A cluster of such variants forms a candidate seed for a haplotype. The evidence from multiple reads across multiple positions are used to identify the putative variants on that haplotype, and then reads are scored based on whether they fall into one of the haplotypes, the other, or cannot be phased. Because this haplotagging uses information from much longer stretches and more candidate variants than the individual process of variant calling, it has the advantage of a broader set of information.

This haplotagging is used to populate the information in the “haplotype channel” which is one of the inputs for DeepVariant long read data. We wrote a blog describing this channel and its impact.

Note that this process is only used to provide the information to the neural network for consider, the neural network will be able to learn when this channel is or is not reliable based on genome context, coverage, etc… the network’s call on the genotype is what finally goes into a variant.

As a result, haplotag is not used as input to generate the non-ref blocks of the gVCF, and as the final variants called are still from the neural network, the definition of a variant remains the same - a position with an ALT allele that receives a non-reference (0/0 or ./.) call.

We are currently working on a deeper description of the phasing logic used in DeepVariant, which may help understand or reproduce the haplotag method more easily.

Please let me know if anything in the explanation is unclear or can be elaborated further.