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 bottomPI-4 is Block 02
. - The phased haplotype in the left within each block is
Hap-A
and on the right isHap-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.
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. $\endgroup$PI
values too which needs to be considered while computing transition counts (and probabilities). $\endgroup$