How to convert the given mathematical computation (on biological problem) to mathematical fomula, equation?

Question

I have crossposted this question in maths StackExchange. The problem is dominantly mathematical (this question) but the application of the problem is mainly biological. Hoping that people in this forum have faced similar problems, I am posting this in bioinformatics forum to get some ideas and probably solutions to it.

I am working on a problem that allows haplotype phasing. I have developed this computation method (expressed in detail below) and have developed a python code (not shown here) to solve the issue.

Now, I am having a hard time trying to translate this computation in mathematical and/or statistical language. I have revisited maths and I am able to comprehend the process at several places, but I cannot translate those ideas onto my computation (or algorithm).

---

Here is my computation (step by step) with a very workable example. I have tried to make it as comprehensive as possible, let me know if there is any confustion:

Fig 01: Example of the haplotype file

contig   pos     ref   all-alleles   ms01e_PI   ms01e_PG_al   ms02g_PI   ms02g_PG_al   ms03g_PI   ms03g_PG_al   ms04h_PI   ms04h_PG_al   ms05h_PI   ms05h_PG_al   ms06h_PI   ms06h_PG_al
2      15881764   .      .           4           C|T           6           C|T           7           T|T           7           T|T           7           C|T           7           C|T
2      15881767   .      .           4           C|C           6           T|C           7           C|C           7           C|C           7           T|C           7           C|C
2      15881989   .      .           4           C|C           6           A|C           7           C|C           7           C|C           7           A|T           7           A|C
2      15882091   .      .           4           G|T           6           G|T           7           T|A           7           A|A           7           A|T           7           A|C
2      15882451   .      .           4           C|T           4           T|C           7           T|T           7           T|T           7           C|T           7           C|A
2      15882454   .      .           4           C|T           4           T|C           7           T|T           7           T|T           7           C|T           7           C|T
2      15882493   .      .           4           C|T           4           T|C           7           T|T           7           T|T           7           C|T           7           C|T
2      15882505   .      .           4           A|T           4           T|A           7           T|T           7           T|T           7           A|C           7           A|T 

This (fig 01) is my main data file. The main idea is to take a pool of several samples which have short phased haplotype blocks represented as PI (phased index) and PG_al (phased genotype). For any sample, the site where two consecutive haplotype blocks are not joined represents the break point (and has different PI values). In the above example only sample ms02g has break point and needs to be phased. All other samples have haplotype block that bridges these two consecutive haplotype blocks and contains the data required to extend the phase state of sample ms02g.

Fig 02: Representing a break-point in the sample: ms02g

contig    pos    ms02g_PI    ms02g_PG_al
2     15881764     6         C|T
2     15881767     6         T|C
2     15881989     6         A|C
2     15882091     6         G|T
                   ×——————————×—————> Break Point
2     15882451     4         T|C
2     15882454     4         T|C 
2     15882493     4         T|C
2     15882505     4         T|A

So, in the above haplotype file there is a breakpoint in sample ms02g at the position 15882091-15882451. For sample ms02g(PI-6) the haplotypes are C-T-A-G and T-C-C-T. Similarly, at PI-4 the haplotypes are T-T-T-T and C-C-C-A. Since the haplotype is broken in two levels, we don’t know which phase from level-6 goes with which phase of level-4. But, all other samples have full haplotype intact that bridges this position. We can therefore use this information from other samples to join the two consecutive haplotype in sample ms02g.

Using markov-chain transition to extend the haplotype block:

Since, all other samples are completely phased bridging that breakpoint, I can run a markov-chain transition probabilities to solve the phase state in the sample ms02g . To the human eye/mind you can clearly see and say that left part of ms02g (PI-6, i.e C-T-A-G is more likely to go with right block of PI-4 C-C-C-A), thereby creating the extended haplotype block as C-T-A-G-C-C-C-A and T-C-C-T-T-T-T-T.

Below, I show how I can apply the first order markov transition matrix to compute the likelyhood estimates, calculate the log2Odds and then assign and extend the haplotype in proper configuration. And, I feed this logic to the computer using python.

Calculation of likelyhood estimates using markov transition:

Step 01: prepare required haplotype configuration

• The top PI-6 is Block 01 and the bottom PI-4 is Block 02.
• The phased haplotype in the left within each block is Hap-A and on the right is Hap-B.

Fig 03: representation of the haplotype breakpoint and block assignment

ms02g_PI    ms02g_PG_al
  6         C|T \
  6         T|C | Block 1
  6         A|C |
  6         G|T /
  ×——————————×—————> Break Point
  4         T|C \
  4         T|C |
  4         T|C | Block 2
  4         T|A /
           ↓   ↓
       Hap-A   Hap-B

So, the two consecutive blocks can be extended in one of the two possible haplotype configurations:

Parallel Configuration:
 Block01-HapA with Block02-HapA, so B01-HapB with B02-HapB 
Vs. Alternate Configuration:
 Block01-HapA with Block02-HapB, so B01-HapB with B02-HapA 

Step 02: Compute transition matrix and estimate likelihood for each configuration.

Fig 04: Representation of the allele transition matrix (from alleles of former block-01 to alleles of later block-02).

Possible
transitions     ms02g_PG_al
 │              ┌┬┬┬  C|T \
 │              ││││  T|C | Block 1
 │              ││││  A|C |
 │              ││││  G|T /
 └────────────> ││││   ×—————> Break Point
                │││└> T|C \
                ││└─> T|C | Block 2
                │└──> T|C |
                └───> T|A /
                     ↓   ↓ 
                 Hap-A   Hap-B

**note: In this example I am showing transition only from first 
  nucleotide of Block01 to first nucleotide of Block02. Actually we 
  prepare transition from all nucleotides of B01 to all nucleotide of 
  B02 in both parallel and alternate configuration. 

• I count the number of transitions from each nucleotide of PI-6 to each nucleotide of PI-4 for each haplotype configuration across all the samples and convert them to transition probabilities.
• And multiply the transition probabilities from the first nucleotide in PI-6 to all nucleotides of PI-4. Then similarly multiply the transition probability from 2nd nucleotide of PI-6 to all nucleotides in PI-4, and so on.
• When transition probabilities are calculated for all possible combination (from each position of PI-6 to each position of PI-4), then I compute the cumulative transition probabilities for each possible haplotype configuration.

Fig 05 : Representation of nucleotide counts (emission counts) at positions 15881764 and 15882451.

pos\allele     A    T    G    C
15881764       0    8    0    4
15882451       1    7    0    4

Fig 06 : Representation of transition matrix counts (from pos 15881764) to (pos 15882451).

This transition matrix is computed from nucleotides (A,T,G,C) at block01 to nucleotides (A,T,G,C) at block02 for all the positions. Transition counts are then converted to transition probabilities.

from     to
          A    T    G    C
  A       0    0    0    0
  T       1    6.5  0    0.5
  G       0    0    0    0
  C       0    0.5  0    3.5

Note: if the PI matches between two blocks the transition are counted as 1, else 0.5. Sample ms02g itself is also taken as an observation.

Step 03: Compute the maximul likelihood for each configuration

Fig 06 : Likelihood estimate for parallel configuration using transition counts (probabilities)

Parallel configuration:
Block-1-Hap-A (C-T-A-G) with Block-2-Hap-A (T-T-T-T)
  CtT × CtT × CtT × CtT = (0.5/4)*(0.5/4)*(0.5/4)*(0.5/4) = 0.000244
+ TtT × TtT × TtT × TtT = (0.5/2)*(0.5/2)*(0.5/2)*(0.5/2) = 0.003906
+ AtT × AtT × AtT × AtT = (0.5/3)*(0.5/3)*(0.5/3)*(0.5/3) = 0.0007716
+ GtT × GtT × GtT × GtT = (0.5/2)*(0.5/2)*(0.5/2)*(0.5/2) = 0.003906
——————— ————————— ——————— ————————— Max Sum (likelihoods) = 0.008828
                                    Average (likelihoods) = 0.002207

Block-1-Hap-B (T-C-C-T) with Block-2-Hap-B (C-C-C-A)
  TtC × TtC × TtC × TtA = (0.5/8)*(0.5/8)*(0.5/8)*(0.5/8) = 0.00001526
+ CtC × CtC × CtC × CtA = (2.5/10)*(2.5/10)*(2.5/10)*(2.5/10) = 0.003906
+ CtC × CtC × CtC × CtA = (1.5/8)*(1.5/8)*(1.5/8)*(1.5/8) = 0.001236
+ TtC × TtC × TtC × TtA = (0.5/4)*(0.5/4)*(0.5/4)*(0.5/4) = 0.000244
——————— ————————— ——————— ————————— Max Sum (likelihoods) = 0.0054016 
                                    Average (likelihoods) = 0.0013504

note: - "AtC" -> represent "A" to "C" transition - "+" represents the summation of the likelyhoods

Fig 07 : Likelihood estimate for alternate configuration using transition counts (probabilities)

Alternate configuration:
Block-1-Hap-A (C-T-A-G) with Block-2-Hap-B (C-C-C-A)
  CtC × CtC × CtC × CtA = (3.5/4)*(3.5/4)*(3.5/4)*(3.5/4) = 0.5861 
+ TtC × TtC × TtC × TtA = (1.5/2)*(1.5/2)*(1.5/2)*(1.5/2) = 0.3164
+ AtC × AtC × AtC × AtA = (2.5/3)*(2.5/3)*(2.5/3)*(2.5/3) = 0.4823
+ GtC × GtC × GtC × GtA = (1.5/2)*(1.5/2)*(1.5/2)*(1.5/2) = 0.3164
——————— ————————— ——————— ————————— Max Sum (likelyhoods) = 1.7012 
                                    Average (likelihoods) = 0.425311

Block-1-Hap-B (T-C-C-T) with Block-2-Hap-A (T-T-T-T) 
  TtC × TtC × TtC × TtA = (6.5/8)*(7.5/8)*(7.5/8)*(6.5/8) = 0.5802
+ CtC × CtC × CtC × CtA = (6.5/10)*(7.5/10)*(7.5/10)*(6.5/10) = 0.237
+ CtC × CtC × CtC × CtA = (5.5/8)*(6.5/8)*(6.5/8)*(6.5/8) = 0.36875
+ TtC × TtC × TtC × TtA = (3.5/4)*(3.5/4)*(3.5/4)*(2.5/4) = 0.4187
——————— ————————— ——————— ————————— Max Sum (likelyhoods) = 1.60465
                                    Average (likelihoods) = 0.4011625

note:

  - the sum of the likelihoods > 1, in the above example.
  - So, we can rather use the product of the likelyhoods.

Fig 08 : Likelhood estimate of Parallel vs. Alternate configuration.

Likelihood of Parallel vs. Alternate configuration

= likelihood of Parallel config / likelihood of Alternate config

= (0.002207 + 0.0013504)/ (0.425311 + 0.4011625)
= 0.0043043  (i.e 1/232)
Therefore, haplotype is 232 times more likely to be phased 
in "alternate configuration".

log2Odds = log2(0.0043043) = -7.860006
So, with our default cutoff threshold of "|5|", we will be extending 
the haplotype in "alternate configuration"

Final Output data:

contig   pos   ref   all-alleles   ms02g_PI   ms02g_PG_al
2   15881764   .    .                  6      C|T
2   15881767   .    .                  6      T|C
2   15881989   .    .                  6      A|C
2   15882091   .    .                  6      G|T
2   15882451   .    .                  6      C|T
2   15882454   .    .                  6      C|T
2   15882493   .    .                  6      C|T
2   15882505   .    .                  6      A|T

So, how do I translate my overall computation process in mathematical and statistical language.

I am looking to represent computation of:

• emission counts (and probabilities)
• transition counts (and probabilities)
• summation of the transition counts for each configuration
• computation of likelihood ratio and it's conversion to log2Odds

I am trying to represent my computation (algorithm) in the way it's shown in these links, but have a hard time doing it. I've read, this, and this, and this.

FYI - I am biologist but have managed to learn python and write a python program to work it out. The code isn't relevant in this question though.

Now, it's my time to work out some maths and represent my algorithm.

Answer 1

I am not sure I understood all, but what you do is like doing a Markov chain but instead of saying that the previous position(s) decides the fate of the next one(s) is the block what determines which other block follows it.

Mathematically I think that could be explained in terms of linkage disequilibrium (and wikipedia page), but instead of using a single position you use a block. I am not familiarly enough with the field to see if a similar (or this method) has been previously described, but I think you could find inspiration about how to describe it.

A related field of work that also deals with this kind of data (and problems) is structural chromatin, where they work with TADs (their blocks of Topologically Associating Domains) and study which ones are near (spatially instead of probabilistic) with others regions, that field could be also interesting to look for formulas and how do they present this kind of analysis.