Answers.Body: <p>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.</p>
  
  <p>Mathematically I think that could be explained in terms of <a href="https://www.nature.com/articles/nrg1123" rel="nofollow noreferrer">linkage disequilibrium</a> (and <a href="https://en.wikipedia.org/wiki/Linkage_disequilibrium" rel="nofollow noreferrer">wikipedia page</a>), 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. </p>
  
  <p>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 <a href="https://en.wikipedia.org/wiki/Topologically_associating_domain" rel="nofollow noreferrer">Topologically Associating Domains</a>) 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.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-04T16:27:39.683Z
Answers.Text: Thanks. But, my problem is entirely different. I have looked into markov chain. That is how I was able to work out that algorithm and write a program. Now, I need a mathematical way to represent my algorithm. I think I can go ahead and read TAD paper in detail - for which I don't really have much time. I am sure there is a better way on working this out. Thanks
Answers.UserDisplayName: everestial
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-04T20:39:40.420Z
Answers.Text: Then could you explain how your problem is different to linkage disequilibrium? When you say you "need a mathematical way to represent my algorithm" do you refer to a mathematical formula like $y = a +b*x$?
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-05T12:01:18.470Z
Answers.Text: My problem is a more complicated than LD, though the very first principle of LD does apply in some sense. In LD you only have 2 genes, 2 alleles and trying to see the segregation pattern. And, even Markov Chains aren't quite involved. It is very easy to represent LD mathematically. On contrary my method has two blocks (length "n" and "m") and then two haplotypes within each block. Then the transitions counts from each "n" to each "m" which are summed. Yes. I just want a mathematical representation. Look for the links in my main question `I've read, this, and this, and this......`
Answers.UserDisplayName: everestial
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-03T20:28:57.417Z
Answers.LastActivityDate: 2018-04-03T20:28:57.417Z
Answers.UserDisplayName: llrs
Body: <p>I have crossposted this question in <a href="https://math.stackexchange.com/q/2716192">maths StackExchange</a>. 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.</p>
  
  <p>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. </p>
  
  <p>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).</p>
  
  <hr>
  
  <p>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:</p>
  
  <p>Fig 01: Example of the haplotype file</p>
  
  <pre><code>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 
  </code></pre>
  
  <p>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 <code>PI (phased index)</code> and  <code>PG_al  (phased genotype)</code>. For any sample, the site where two consecutive haplotype blocks are not joined represents the break point (and has different <code>PI</code> values). In the above example only sample <code>ms02g</code> 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 <code>sample ms02g</code>. </p>
  
  <p>Fig 02: Representing a break-point in the sample: ms02g</p>
  
  <pre><code>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
                     ×——————————×—————&gt; Break Point
  2     15882451     4         T|C
  2     15882454     4         T|C 
  2     15882493     4         T|C
  2     15882505     4         T|A
  </code></pre>
  
  <p>So, in the above haplotype file there is a breakpoint in <code>sample ms02g</code> at the position <em>15882091-15882451</em>. For sample  ms02g(PI-6) the haplotypes are  <code>C-T-A-G</code> and  <code>T-C-C-T</code>. Similarly, at PI-4 the haplotypes are <code>T-T-T-T</code> and <code>C-C-C-A</code>. 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.</p>
  
  <p>Using markov-chain transition to extend the haplotype block: </p>
  
  <p>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 <code>sample ms02g</code> . 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  <code>C-T-A-G-C-C-C-A</code> and <code>T-C-C-T-T-T-T-T</code>.</p>
  
  <p>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.</p>
  
  <p>Calculation of likelyhood estimates using markov transition:</p>
  
  <h3>Step 01:   prepare required haplotype configuration</h3>
  
  <ul>
  <li>The top <code>PI-6 is Block 01</code> and the bottom <code>PI-4 is Block 02</code>.</li>
  <li>The phased haplotype in the left within each block is <code>Hap-A</code> and on the right is <code>Hap-B</code>.</li>
  </ul>
  
  <p>Fig 03: representation of the haplotype breakpoint and block assignment</p>
  
  <pre><code>ms02g_PI    ms02g_PG_al
    6         C|T \
    6         T|C | Block 1
    6         A|C |
    6         G|T /
    ×——————————×—————&gt; Break Point
    4         T|C \
    4         T|C |
    4         T|C | Block 2
    4         T|A /
             ↓   ↓
         Hap-A   Hap-B
  </code></pre>
  
  <p>So, the two consecutive blocks can be extended in one of the two possible haplotype configurations:</p>
  
  <pre><code>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 
  </code></pre>
  
  <h3>Step 02: Compute transition matrix and estimate likelihood for each configuration.</h3>
  
  <p>Fig 04: Representation of the allele transition matrix (from alleles of former block-01 to alleles of later block-02).</p>
  
  <pre><code>Possible
  transitions     ms02g_PG_al
   │              ┌┬┬┬  C|T \
   │              ││││  T|C | Block 1
   │              ││││  A|C |
   │              ││││  G|T /
   └────────────&gt; ││││   ×—————&gt; Break Point
                  │││└&gt; T|C \
                  ││└─&gt; T|C | Block 2
                  │└──&gt; T|C |
                  └───&gt; 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. 
  </code></pre>
  
  <ul>
  <li>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. </li>
  <li>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.</li>
  <li>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.</li>
  </ul>
  
  <p>Fig 05 : Representation of nucleotide counts (emission counts) at positions 15881764 and 15882451.</p>
  
  <pre><code>pos\allele     A    T    G    C
  15881764       0    8    0    4
  15882451       1    7    0    4
  </code></pre>
  
  <p>Fig 06 : Representation of transition matrix counts (from pos 15881764) to (pos 15882451).</p>
  
  <p>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.</p>
  
  <pre><code>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
  </code></pre>
  
  <p>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.</p>
  
  <h3>Step 03: Compute the maximul likelihood for each configuration</h3>
  
  <p>Fig 06 : Likelihood estimate for parallel configuration using transition counts (probabilities)</p>
  
  <pre><code>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
  </code></pre>
  
  <p>note: 
        - "AtC" -> represent "A" to "C" transition
        - "+" represents the summation of the likelyhoods</p>
  
  <p>Fig 07 : Likelihood estimate for alternate configuration using transition counts (probabilities)</p>
  
  <pre><code>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
  </code></pre>
  
  <p>note: </p>
  
  <pre><code>  - the sum of the likelihoods &gt; 1, in the above example.
    - So, we can rather use the product of the likelyhoods.
  </code></pre>
  
  <p>Fig 08 : Likelhood estimate of Parallel vs. Alternate configuration.</p>
  
  <pre><code>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"
  </code></pre>
  
  <h3>Final Output data:</h3>
  
  <pre><code>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
  </code></pre>
  
  <p>So, how do I translate my overall computation process in mathematical and statistical language. </p>
  
  <p>I am looking to represent computation of:</p>
  
  <ul>
  <li>emission counts (and probabilities)</li>
  <li>transition counts (and probabilities)</li>
  <li>summation of the transition counts for each configuration</li>
  <li>computation of likelihood ratio and it's conversion to log2Odds</li>
  </ul>
  
  <p>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, <a href="https://www.probabilitycourse.com/chapter11/11_2_3_probability_distributions.php" rel="nofollow noreferrer">this</a>, and <a href="https://stats.stackexchange.com/questions/25540/summation-of-conditional-probability-in-bayesian-network">this</a>, and <a href="https://math.stackexchange.com/questions/389062/how-to-write-an-algorithm-as-an-equation">this</a>.</p>
  
  <p>FYI - I am biologist but have managed to learn python and write a python program to <a href="https://codereview.stackexchange.com/questions/186396/solve-the-phase-state-between-two-haplotype-blocks-using-markov-transition-proba">work it out</a>. The code isn't relevant in this question though.</p>
  
  <p>Now, it's my time to work out some maths and represent my algorithm. </p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-31T22:06:02.247Z
Comments.Text: Hi, I have edited a bit the question, I think that now is easier to read, feel free to edit further or roll back if you disagree. what is encoded in the haplotype columns (Fig 01) ? It seems that you have like several samples in the same file, and you are only interested in one sample. Have you looked at how do the mathematics represent markov chain  models?
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-31T23:57:32.400Z
Comments.Text: @Llopis: Yes, I am trying to solve one sample at a time; in this case `ms02g` - because this is the one that has a breakpoint (i.e different PI blocks). Other samples have full length haplotype, so they actually serve has a observation data for preparing the markov transition matrix.
Comments.UserDisplayName: everestial
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-04-01T00:04:25.450Z
Comments.Text: @Llopis: Yes, I looked into how to represent markov model. I can comprehend at several places but all. So, I cannot translate. Also, my data has `PI` values too which needs to be considered while computing transition counts (and probabilities).
Comments.UserDisplayName: everestial
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-04-04T16:24:15.223Z
Comments.Text: Are you asking how to take what you want to do and write a program to do it, or are you asking, given an algorithm, how to do a write-up of your algorithm to put in a paper?
Comments.UserDisplayName: Acccumulation
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-04-04T16:25:27.663Z
Comments.Text: I have already written a program. I want to be able to write it in a paper. I have left maths for more than a decade and seems like I will be quite slower if I tried to re-read the whole thing again. So, any way you can help is appreciable.
Comments.UserDisplayName: everestial
Tags: <python><bioconductor><hidden-markov-models><haplotypes>
Title: How to convert the given mathematical computation (on biological problem) to mathematical fomula, equation?
ThreadText: # How to convert the given mathematical computation (on biological problem) to mathematical fomula, equation?
  
  (Asked by: username_0 on 2018-03-31)
  
  I have crossposted this question in [maths StackExchange](https://math.stackexchange.com/q/2716192). 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](https://www.probabilitycourse.com/chapter11/11_2_3_probability_distributions.php), and [this](https://stats.stackexchange.com/questions/25540/summation-of-conditional-probability-in-bayesian-network), and [this](https://math.stackexchange.com/questions/389062/how-to-write-an-algorithm-as-an-equation).
  
  FYI - I am biologist but have managed to learn python and write a python program to [work it out](https://codereview.stackexchange.com/questions/186396/solve-the-phase-state-between-two-haplotype-blocks-using-markov-transition-proba). The code isn't relevant in this question though.
  
  Now, it's my time to work out some maths and represent my algorithm.
  
  **Tags:** python, bioconductor, hidden-markov-models, haplotypes
  
  > **Comments:**
  > **username_1**: Hi, I have edited a bit the question, I think that now is easier to read, feel free to edit further or roll back if you disagree. what is encoded in the haplotype columns (Fig 01) ? It seems that you have like several samples in the same file, and you are only interested in one sample. Have you looked at how do the mathematics represent markov chain  models?  
  > **username_0**: @Llopis: Yes, I am trying to solve one sample at a time; in this case `ms02g` - because this is the one that has a breakpoint (i.e different PI blocks). Other samples have full length haplotype, so they actually serve has a observation data for preparing the markov transition matrix.  
  > **username_0**: @Llopis: Yes, I looked into how to represent markov model. I can comprehend at several places but all. So, I cannot translate. Also, my data has `PI` values too which needs to be considered while computing transition counts (and probabilities).  
  > **username_2**: Are you asking how to take what you want to do and write a program to do it, or are you asking, given an algorithm, how to do a write-up of your algorithm to put in a paper?  
  > **username_0**: I have already written a program. I want to be able to write it in a paper. I have left maths for more than a decade and seems like I will be quite slower if I tried to re-read the whole thing again. So, any way you can help is appreciable.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2018-04-03. Score: 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](https://www.nature.com/articles/nrg1123) (and [wikipedia page](https://en.wikipedia.org/wiki/Linkage_disequilibrium)), 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](https://en.wikipedia.org/wiki/Topologically_associating_domain)) 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.
  
  > **Comments:**
  > **username_0**: Thanks. But, my problem is entirely different. I have looked into markov chain. That is how I was able to work out that algorithm and write a program. Now, I need a mathematical way to represent my algorithm. I think I can go ahead and read TAD paper in detail - for which I don't really have much time. I am sure there is a better way on working this out. Thanks  
  > **username_1**: Then could you explain how your problem is different to linkage disequilibrium? When you say you "need a mathematical way to represent my algorithm" do you refer to a mathematical formula like $y = a +b*x$?  
  > **username_0**: My problem is a more complicated than LD, though the very first principle of LD does apply in some sense. In LD you only have 2 genes, 2 alleles and trying to see the segregation pattern. And, even Markov Chains aren't quite involved. It is very easy to represent LD mathematically. On contrary my method has two blocks (length "n" and "m") and then two haplotypes within each block. Then the transitions counts from each "n" to each "m" which are summed. Yes. I just want a mathematical representation. Look for the links in my main question `I've read, this, and this, and this......`  
  
  ---
  
  