Revise probabilities in generateErrorMat‎()

[gregorgorjanc]

Looking at generateErrorMat‎() we have

https://github.com/AlphaGenes/tinyhouse./blob/65456ea92e58fc40886fa1d1654b3a66443f28b4/ProbMath.py#L216C30-L216C30

def generateErrorMat(error) :
    errorMat = np.array([[1-error, error/2, error/2, error/2], 
                            [error/2, 1-error, 1-error, error/2],
                            [error/2, error/2, error/2, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 1)[:,None]
    return errorMat

I am wondering if the probabilities should actually be [1-error, error/4, error/4, error/2] so that we spread the error equally across the 3 genotypes or give user an option to declare these genotype penetrance matrices. The current probabilities are for sure odd - we take error amount from one genotype and then redistribute error/2 to the other three genotypes - aha are we doing this effectively so that we give error/2 to hets and error/2 to the other homozygote? Technically it's correct due to normalisation (in errorMat/np.sum(errorMat, 1)[:,None]).

[gregorgorjanc]

I am wondering if the probabilities should actually be [1-error, error/4, error/4, error/2] so that we spread the error equally across the 3 genotypes or give user an option to declare these genotype penetrance matrices.

Dimensions here are observed genotype (with 3 states in rows) times true phased genotyped (with 4 states in columns). Penetrance function defines a probability distribution of observed genotype (with 3 states in row) for each true state. So, we need to look at columns of this matrix as one conditional distribution. Then these values make more sense, BUT, running the code recently with @AprilYUZhang we think the resulting probs are indeed odd (alternative observed homozygote had a bit larger probability than observed heterozygous genotype when true genotype was reference homozgyote - was this done on purpose or is this a mistake!?). We should review these functions and as mentioned in AlphaGenes/AlphaPeel#47 we should enable users to have control over these error probabilities.

[AprilYUZhang]

I agree, this change is clearer.
Here is the output.

       [0.005 , 0.495 , 0.495 , 0.005 ],
       [0.005 , 0.0025, 0.0025, 0.99  ]], dtype=float32)

It looks easy to explain that when we have genotype 1, the actual haplotype is AA Aa aA aa -> the probability should be [0.005 , 0.495 , 0.495 , 0.005 ]

def generateErrorMat(error) :
    errorMat = np.array([[1-error, error/4, error/4, error/2], 
                            [error/2, (1-error)/2, (1-error)/2, error/2],
                            [error/2, error/4, error/4, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 1)[:,None]
    #errorMat = errorMat / np.sum(errorMat, axis=1, keepdims=True)

    return errorMat

[gregorgorjanc]

@AprilYUZhang you are missing one row in the first output. Please correct so we can have fully informed discussion here. Ping me when you do that.

[AprilYUZhang]

I am afraid that I will change my idea. Now Alphapeel carried the right thing. I think if I change this matrix to two matrix multiplications, it will be easier to understand.
It is the process from the observed genotype mapping to the actual genotype and then projecting to actual haplotype value. The details are as follows.

The matrix projecting genotype to haplotype is key, if changing this mapping function to the following matrix, we can get the error matrix you arose.

$$ \begin{bmatrix} 1 & 0 &0 & 0 \\ 0 & 0.5 & 0.5 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} $$

here is my code

errorMat=np.matmul(np.array([[1-error, error/2, error/2], 
[error/2, 1-error, error/2], 
[error/2, error/2, 1-error]]), 
np.array([[1, 0, 0,0], 
[0, 1, 1,0], 
[0, 0, 0, 1]]))
errorMat = errorMat/np.sum(errorMat, 1)[:,None]
errorMat
# output 
array([[0.98507463, 0.00497512, 0.00497512, 0.00497512],
       [0.00251256, 0.49748744, 0.49748744, 0.00251256],
       [0.00497512, 0.00497512, 0.00497512, 0.98507463]])

errorMat=np.matmul(np.array([[1-error, error/2, error/2], 
[error/2, 1-error, error/2], 
[error/2, error/2, 1-error]]), 
np.array([[1, 0, 0,0], 
[0, 0.5, 0.5,0], 
[0, 0, 0, 1]]))
errorMat = errorMat/np.sum(errorMat, 1)[:,None]
errorMat
# output 
array([[0.99  , 0.0025, 0.0025, 0.005 ],
       [0.005 , 0.495 , 0.495 , 0.005 ],
       [0.005 , 0.0025, 0.0025, 0.99  ]])

[AprilYUZhang]

I think this Counter-intuitive thing (alternative observed homozygote had a bit larger probability than observed heterozygous genotype when true genotype was reference homozygote) results in normalization, the real error matrix used to invert the projection should be that before normalisation. The normalisation breaks the invertible process.

# real error matrix used to invert the projection
array([[0.99 , 0.005, 0.005, 0.005],
       [0.005, 0.99 , 0.99 , 0.005],
       [0.005, 0.005, 0.005, 0.99 ]])

It also got me thinking, is this normalization necessary? Although in the biologist's mind, the sum of the four haplotypes should be 1, but this is probability, and I think there should be no such constraints. And I'm still thinking about whether normalize affects our results, it seem to have no impact so far

[gregorgorjanc]

@AprilYUZhang let's see.

>>> import numpy as np
... def generateErrorMat(error) :
...     errorMat = np.array([[1-error, error/4, error/4, error/2], 
...                             [error/2, (1-error)/2, (1-error)/2, error/2],
...                             [error/2, error/4, error/4, 1-error]], dtype = np.float32)
...     errorMat = errorMat/np.sum(errorMat, 1)[:,None]
...     return errorMat
... 
... generateErrorMat(error=0.01)
array([[0.99  , 0.0025, 0.0025, 0.005 ],
       [0.005 , 0.495 , 0.495 , 0.005 ],
       [0.005 , 0.0025, 0.0025, 0.99  ]], dtype=float32)

The key point here is that rows are associated with observed genotypes [0, 1, 2], and columns are associated with true phased genotypes [aa, aA, Aa, AA].

My understanding is that penetrance function defines conditional distribution of observed genotypes [0, 1, 2] given the true genotypes (in our case true phased genotypes [aa, aA, Aa, AA]).

Your thinking was to first look at penetrance function of observed genotypes [0, 1, 2] given the true non-phased genotypes [aa, aA or Aa, AA]), which would indeed be:

>>> def generateErrorMatNonPhased(error) :
...     errorMat = np.array([[1-error, error/2, error/2], 
...                             [error/2, (1-error), error/2],
...                             [error/2, error/2, 1-error]], dtype = np.float32)
...     errorMat = errorMat/np.sum(errorMat, 1)[:,None]
...     return errorMat
... 
... generateErrorMatNonPhased(error=0.01)
array([[0.99 , 0.005, 0.005],
       [0.005, 0.99 , 0.005],
       [0.005, 0.005, 0.99 ]], dtype=float32)

I am too puzzled by why Andrew split the heterozygous case in two AND then divided the probability by 2 as well!? Maybe this is how the algorithm is worked out internally, but my first intuition tells me that penetrance function is saying "given this true phased genotype (say aA) we can expect these kind of genotypes". So I would expect code like this, where we have a column of probabilities that add up to 1 (since each column is a conditional distribution of an outcome (=observed genotype) given the input (=true phased genotype)):

>>> def generateErrorMatPhased(error) :
...     errorMat = np.array([[1-error, error/2, error/2, error/2], 
...                             [error/2, 1-error, 1-error, error/2],
...                             [error/2, error/2, error/2, 1-error]], dtype = np.float32)
...     errorMat = errorMat/np.sum(errorMat, 1)[:,None]
...     return errorMat
... 
>>> generateErrorMatPhased(error=0.01)
array([[0.98507464, 0.00497512, 0.00497512, 0.00497512],
       [0.00251256, 0.49748743, 0.49748743, 0.00251256],
       [0.00497512, 0.00497512, 0.00497512, 0.98507464]], dtype=float32)

Hmm, nope, I did not expect this! These columns do not sum up to 1 and as we started discussion above, if true phased genotype is aa, why are we assuming that we should expect observed genotypes [0, 1, 2] with probability of [0.985, 0.002, 0.005] - I am confused in particular why would we assume that aa is more likely to lead to 2 (with prob of 0.005) than to lead to 1 (with prob of 0.002)! What is going on here?

>>> error = 0.01
... errorMat = np.array([[1-error, error/2, error/2, error/2],
...                      [error/2, 1-error, 1-error, error/2],
...                      [error/2, error/2, error/2, 1-error]], dtype = np.float32)
... errorMat
... np.sum(errorMat, 1)
array([1.005, 1.99 , 1.005], dtype=float32)

Alright, this summation is not what we want (I think)!? I was expecting to get the sum for each column, not row. I was expecting this:

>>> error = 0.01
... errorMat = np.array([[1-error, error/2, error/2, error/2],
...                      [error/2, 1-error, 1-error, error/2],
...                      [error/2, error/2, error/2, 1-error]], dtype = np.float32)
... errorMat
... np.sum(errorMat, 0)
array([1., 1., 1., 1.], dtype=float32)

Which shows that we have already defined every column as a proper distribution (=the sum of probs is 1), but normalisation is OK to ensure we indeed get the sum of 1.

>>> errorMat/np.sum(errorMat, 0)
array([[0.99 , 0.005, 0.005, 0.005],
       [0.005, 0.99 , 0.99 , 0.005],
       [0.005, 0.005, 0.005, 0.99 ]], dtype=float32)

Note, I changed np.sum(errorMat, 1)[:,None] to np.sum(errorMat, 0).

So, I would expect then to see function like this:

>>> def generateErrorMatPhased(error) :
...     errorMat = np.array([[1-error, error/2, error/2, error/2], 
...                             [error/2, 1-error, 1-error, error/2],
...                             [error/2, error/2, error/2, 1-error]], dtype = np.float32)
...     errorMat = errorMat/np.sum(errorMat, 0)
...     return errorMat
>>> generateErrorMatPhased(error=0.01)
array([[0.99 , 0.005, 0.005, 0.005],
       [0.005, 0.99 , 0.99 , 0.005],
       [0.005, 0.005, 0.005, 0.99 ]], dtype=float32)

This now suggests that if true phased genotype is aa, we are assuming that we should expect observed genotypes [0, 1, 2] with probability of [0.990, 0.005, 0.005], which makes more sense in my head.

An even better would probably be:

>>> def generateErrorMatPhased2(error):
...     errorMat = np.array([[(1-error)**2,      0, 0, 0], 
...                          [2*(1-error)*error, 0, 0, 0],
...                          [error**2,          0, 0, 0]], dtype = np.float32)
...     errorMat = errorMat/np.sum(errorMat, 0)
...     return errorMat
... 
... generateErrorMatPhased2(error=0.01)
<positron-console-cell-84>:5: RuntimeWarning: invalid value encountered in divide
array([[9.801e-01,       nan,       nan,       nan],
       [1.980e-02,       nan,       nan,       nan],
       [1.000e-04,       nan,       nan,       nan]], dtype=float32)

I have done only one column as I need to move onwards. But this now says that if true phased genotype is aa, we are assuming that we should expect observed genotypes [0, 1, 2] with probability of [0.9801, 0.0198, 0.0001], which suggest that going from aa to 0 is the most likely outcome (with prob 0.9801), from aa to 1 is then possible (with prob 0.0198) and from aa to 2 is the least likely (with prob 0.0001).

@AprilYUZhang replace the 0 in generateErrorMatPhased2() above following the above logic and also see OneNote notes I shared with you.

I have two questions in my mind now:

1. Why has Andrew split the penetrance function in two columns (this part is OK since we want to separate the heterozygous case into aA and Aa phased genotypes) AND divided the probability of the heterozygous case by 2 - I think that speaking purely from the penetrance function we should look at input (and input is only either aA or Aa) and output and each column of this penetrance function should be a proper distribution (probs sum to 1). Has he done the above because there is something algorithmic in the code that requires the above confusing setting?

2. The normalisation with np.sum(errorMat, 1)[:,None] instead of with np.sum(errorMat, 0) looks like an error to me, but again, has he done the above because there is something algorithmic in the code that requires the above confusing setting?

I wonder if we have genuine 2 bugs in here and this will influence a lot of things down the road!?

@XingerTang do you have any thoughts on the above? You have been going through the code & algorithms quite a bit!

@RosCraddock can you please rerun one of your analysis with replacing the current function

def generateErrorMat(error) :
    errorMat = np.array([[1-error, error/4, error/4, error/2], 
                            [error/2, (1-error)/2, (1-error)/2, error/2],
                            [error/2, error/4, error/4, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 1)[:,None]
    return errorMat

with

def generateError(error):
    errorMat = np.array([[1-error, error/2, error/2, error/2], 
                            [error/2, 1-error, 1-error, error/2],
                            [error/2, error/2, error/2, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 0)
    return errorMat

Does this change any of the trends you see in your work? Please look trends in allele freq, geno freq, phased geno freq, and estimates of model parameters (base pop allele freq trends, error rate, anything else?).

We can later try to see if moving to generateErrorMatPhased2() would be even better.

[RosCraddock]

@gregorgorjanc I assume you meant change this:

def generateErrorMat(error) :
    errorMat = np.array([[1-error, error/2, error/2, error/2], 
                            [error/2, 1-error, 1-error, error/2],
                            [error/2, error/2, error/2, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 1)[:,None]
    return errorMat

to this:

def generateErrorMat(error) :
    errorMat = np.array([[1-error, error/4, error/4, error/2], 
                            [error/2, (1-error)/2, (1-error)/2, error/2],
                            [error/2, error/4, error/4, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 0)
    return errorMat

I did this in two stages: first changing to
errorMat = errorMat/np.sum(errorMat, 0),
then
errorMat = np.array([[1-error, error/4, error/4, error/2], [error/2, (1-error)/2, (1-error)/2, error/2], [error/2, error/4, error/4, 1-error]], dtype = np.float32)

Rerunning my analysis only provided a very small difference in the Miniature Bull Terrier.

• Genotype probability error marginally reduced from 0.004164 to 0.004160.
• The disease-causing allele frequency showed a marginally higher frequency in 2002 to 2008 (~0.001 dif each year) and lower frequency in 2018 to 2024 (<0.001 dif each year).

These changes were observed with errorMat = errorMat/np.sum(errorMat, 0) and maintained with

errorMat = np.array([[1-error, error/4, error/4, error/2], [error/2, (1-error)/2, (1-error)/2, error/2], [error/2, error/4, error/4, 1-error]], dtype = np.float32).

There was no difference in any of the phased genotype probability trends.

[AprilYUZhang]

Here is the phased genotype matrix.

def generateErrorMatPhased2(error):
    errorMat = np.array([[(1-error)**2,      (1-error)*error, (1-error)*error, error**2], 
                        [2*(1-error)*error, (1-error)**2+error**2, (1-error)**2+error**2, 2*(1-error)*error],
                        [error**2,          (1-error)*error, (1-error)*error, (1-error)**2]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 0)
    return errorMat

generateErrorMatPhased2(error=0.01)

# output
array([[9.8009998e-01, 9.9000009e-03, 9.9000009e-03, 9.9999997e-05],
       [1.9800000e-02, 9.8020005e-01, 9.8020005e-01, 1.9800000e-02],
       [9.9999997e-05, 9.9000009e-03, 9.9000009e-03, 9.8009998e-01]],
      dtype=float32)

[gregorgorjanc]

@AprilYUZhang your generateErrorMatPhased2() looks like it's giving back an array that I would expect.

[gregorgorjanc]

@RosCraddock good to see that you see improvement, even if small, with likely more correct normalisation.

@gregorgorjanc I assume you meant change this:

def generateErrorMat(error) :
    errorMat = np.array([[1-error, error/2, error/2, error/2], 
                            [error/2, 1-error, 1-error, error/2],
                            [error/2, error/2, error/2, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 1)[:,None]
    return errorMat

to this:

def generateErrorMat(error) :
    errorMat = np.array([[1-error, error/4, error/4, error/2], 
                            [error/2, (1-error)/2, (1-error)/2, error/2],
                            [error/2, error/4, error/4, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 0)
    return errorMat

I believe that the original function is

>>> import numpy as np
... def generateErrorMat(error) :
...     errorMat = np.array([[1-error, error/4, error/4, error/2], 
...                             [error/2, (1-error)/2, (1-error)/2, error/2],
...                             [error/2, error/4, error/4, 1-error]], dtype = np.float32)
...     errorMat = errorMat/np.sum(errorMat, 1)[:,None]
...     return errorMat
... 
... generateErrorMat(error=0.01)
array([[0.99  , 0.0025, 0.0025, 0.005 ],
       [0.005 , 0.495 , 0.495 , 0.005 ],
       [0.005 , 0.0025, 0.0025, 0.99  ]], dtype=float32)

and I am confused why the two intermediate columns don't give proper conditional distribution that sum to 1. If I sum both columns together, then we do get a sum of 1, but this conditional distribution should be proper for each phased genotype.

Hence, I would appreciate if you can test changing:

def generateErrorMat(error) :
    errorMat = np.array([[1-error, error/4, error/4, error/2], 
                            [error/2, (1-error)/2, (1-error)/2, error/2],
                            [error/2, error/4, error/4, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 0)
    return errorMat

to

def generateError(error):
    errorMat = np.array([[1-error, error/2, error/2, error/2], 
                            [error/2, 1-error, 1-error, error/2],
                            [error/2, error/2, error/2, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 0)
    return errorMat

The latter makes more sense to me.

Would be good to hear about this before we meet next week.

[RosCraddock]

To summarise, today's meeting (27/03/2025) focused on the generateErrorMat function:

def generateErrorMat(error) :
    errorMat = np.array([[1-error, error/2, error/2, error/2], 
                        [error/2, 1-error, 1-error, error/2],
                        [error/2, error/2, error/2, 1-error]], dtype = np.float32)
    errorMat = errorMat/np.sum(errorMat, 1)[:,None]
    return errorMat

When true (g) -> observed (og), which is equivalent to the probability of observed given the true probability (Pr(og|g)), then all the columns need to sum to 1. (observed by row, true by column)
But, when observed -> true, Pr(g|og), all the rows need to sum to 1. Currently, the generateErrorMat(error) function does the latter using normalisation. This normalisation is necessary since the setGenoProbsFromGenotypes() function grabs the row which matches the observed genotype to inform the genotype probabilities.

Genotype error Vs allele error

Currently in AlphaPeel, each genotype state is equally likely in error. For example, if the observed genotype is AA, Aa and aa are equally likely to be the true genotype. Although, you would expect Aa to be more likely than aa.

If an error per allele is used (as opposed to per genotype), this can be modelled as previously described by Yu:

    errorMat = np.array([[(1-error)**2,      (1-error)*error, (1-error)*error, error**2], 
                        [2*(1-error)*error, (1-error)**2+error**2, (1-error)**2+error**2, 2*(1-error)*error],
                        [error**2,          (1-error)*error, (1-error)*error, (1-error)**2]], dtype = np.float32)

This gives different weightings for each genotype state. However, the default error and estimation of genotype error is based on genotype rather than allele. If this is changed, we will need to update this, too.

Alternatively, we could provide different weightings for the genotype error across states. Thus, only the generateErrorMat() function would be changed.

Simulation of genotype error for testing

I have inserted the code below, but genotype error is added randomly based on a genotype error probability in the simulations used for testing. Each genotype state is weighted the same. So, we may need additional/adapt tests to consider the different error weights across genotype states.

generateGenoErr <- function(geno, error) {
  nLoci <- ncol(geno)
  for (ind in 1:nInd){
    for (locus in 1:nLoci){
      if (geno[ind, locus] != 9 && rbinom(n = 1, size = 1, prob = error) == 1) {
        if (geno[ind, locus] == 0) {
          geno[ind, locus] <- sample(c(1, 2), size = 1)
        } else if (geno[ind, locus] == 1) {
          geno[ind, locus] <- sample(c(0, 2), size = 1)
        } else {
          geno[ind, locus] <- sample(c(0, 1), size = 1)
        }
      }
    }
  }
  return(geno)
}