Difference between revisions of "Relationship between Ploidy, Alleles and Genotypes"
(Created page with "= Introduction = Ploidy, alleles are independent characteristics that give rise to genotypes. = Table = { class="wikitable"  ! scope="col" Case ! scope="col" Alleles ...") 
(→Motivation) 

(101 intermediate revisions by the same user not shown)  
Line 1:  Line 1:  
= Introduction =  = Introduction =  
−  +  The VCF format encodes genotypes by the index of the enumeration of genotypes given ploidy number and alleles.  
+  This allows for direct access to a value associated with a genotype within an array when one works with genotype likelihoods.  
−  =  +  = Motivation = 
+  Plants species exhibit a diverse number of ploidy, for example, the strawberry is an octoploid and the pear is a triploid.  
+  
+  Copy number variations in somatic variant calling also leads to variable ploidy to consider when genotyping a locus.  
+  
+  While there are explicit functions that could be googled for handling haploid and diploid cases. It seems difficult to find the closed forms for the general case.  
+  This wiki fills in that need.  
+  
+  = The number of genotypes given a ploidy and alleles =  
+  
+  <math>  
+  \begin{align}  
+  F(P,A) = \binom{P+A1}{A1} \\  
+  \end{align}  
+  </math>  
+  
+  where P is the ploidy number and A is the number of alleles.  
+  
+  = Getting the index of a genotype in an enumerated list given a ploidy and alleles =  
+  
+  <math>  
+  \begin{align}  
+  G(a_1,.. , a_P) = \sum_{k=1}^P \binom{k+a_k1}{a_k1}  
+  \end{align}  
+  </math>  
+  
+  where P is the number of ploidy, <math>a_1</math>, <math>a_2</math> .. <math>a_P</math> are the alleles in numeric encoding (0 to A1) and are ordered (e.g. AB and ABCCCC are ordered but ACB is not ordered).  
+  
+  This is well defined because:  
+  
+  <math>  
+  \begin{align}  
+  \binom{n}{r} = \begin{cases}  
+  \frac{n!}{(nr)!r!} &, r \le n, r\ge0, n\ge0 \\  
+  0 & \text{otherwise}  
+  \end{cases}  
+  \end{align}  
+  </math>  
+  
+  Because <math>a_k</math> may be 0, we will see cases of <math>\binom{k1}{1}</math> when <math>a_k=0</math>. This is alright because of the definition of <math>\binom{n}{r}</math> which defines this case as 0.  
+  But to make it more sensible, we can define the function equivalently as:  
+  
+  <math>  
+  \begin{align}  
+  G(a_1,.. , a_P) = \sum_{k=1}^P \binom{k+a_k1}{k}  
+  \end{align}  
+  </math>  
+  
+  So when <math>a_k=0</math>, the binomial coefficient reads as <math>\binom{k1}{k}</math> which equals 0 since there are 0 ways to choose k items from k1 items.  
+  
+  = Getting the genotypes from a genotype index and a given ploidy =  
+  
+  The genotype index is computed by a summation of a series which is  
+  monotonically decreasing. This allows you to compute the inverse function  
+  from index to the ordered genotypes by using a "water rapids algorithm" with  
+  decreasing height of each mini water fall.  
+  
+  std::vector<int32_t> bcf_ip2g(int32_t genotype_index, uint32_t no_ploidy)  
+  {  
+  std::vector<int32_t> genotype(no_ploidy, 0);  
+  int32_t pth = no_ploidy;  
+  int32_t max_allele_index = genotype_index;  
+  int32_t leftover_genotype_index = genotype_index;  
+  while (pth>0)  
+  {  
+  for (int32_t allele_index=0; allele_index <= max_allele_index; ++allele_index)  
+  {  
+  int32_t i = choose(pth+allele_index1, pth);  
+  if (i>=leftover_genotype_index  allele_index==max_allele_index)  
+  {  
+  if (i>leftover_genotype_index) allele_index;  
+  leftover_genotype_index = choose(pth+allele_index1, pth);  
+  pth;  
+  max_allele_index = allele_index;  
+  genotype[pth] = allele_index;  
+  break;  
+  }  
+  }  
+  }  
+  return genotype;  
+  }  
+  
+  todo:: describe in a human understandable fashion.  
+  
+  = Simple cases =  
{ class="wikitable"  { class="wikitable"  
    
−  ! scope="col"  +  ! scope="col" Ploidy 
! scope="col" Alleles  ! scope="col" Alleles  
−  ! scope="col" Genotypes  +  ! scope="col" No. of Genotypes 
! scope="col" Index  ! scope="col" Index  
−  
    
−    +   1 
−    +   A 
+   <math>  
+  F(1, A) = \binom{1+A1}{A1} = A  
+  </math>  
+   <math>  
+  G(a_1) = F(1, a_1) = a_1  
+  </math>  
+    
+   2  
+   A  
+   <math>  
+  F(2,A) = \binom{2+A1}{A1} = \binom{A+1}{2}  
+  </math>  
+   <math>  
+  G(a_1,a_2) = F(1, a_1) + F(2, a_2) = a_1 + \binom{a_2+1}{2}  
+  </math>  
+  }  
+  
+  = Derivation for counting the number of genotypes =  
+  
+  There must always be P observed alleles and there can only be at most A alleles. This can be modeled by P+A1 points where you choose A1 points to be dividers that separate the alleles.  
+  Thus the number of ways you can observe this is <math> \binom{P+A1}{A1} </math> which is equivalent to <math> \binom{P+A1}{P} </math>.  
+  
+  = Derivation for getting the index of a genotype in an enumerated list =  
+  
+  == Observation of nested patterns ==  
+  An important observation here is that for the enumeration of A alleles for a given P ploidy, the enumeration of A1 alleles for P ploidy is a subsequence.  
+  
+  { class="wikitable"  
+    
+  ! scope="col" Index  
+  ! scope="col" A=4,P=3  
+  ! scope="col" A=3,P=3  
+  ! scope="col" A=2,P=3  
+  ! scope="col" A=1,P=3  
+    
+  0 <br>  
+  1 <br>  
+  2 <br>  
+  3 <br>  
+  4 <br>  
+  5 <br>  
+  6 <br>  
+  7 <br>  
+  8 <br>  
+  9 <br>  
+  10 <br>  
+  11 <br>  
+  12 <br>  
+  13 <br>  
+  14 <br>  
+  15 <br>  
+  16 <br>  
+  17 <br>  
+  18 <br>  
+  19 <br>  
+   AAA <br>  
+  AAB <br>  
+  ABB <br>  
+  BBB <br>  
+  AAC <br>  
+  ABC <br>  
+  BBC <br>  
+  ACC <br>  
+  BCC <br>  
+  CCC <br>  
+  AAD <br>  
+  ABD <br>  
+  BBD <br>  
+  ACD <br>  
+  BCD <br>  
+  CCD <br>  
+  ADD <br>  
+  BDD <br>  
+  CDD <br>  
+  DDD <br>  
+   AAA <br>  
+  AAB <br>  
+  ABB <br>  
+  BBB <br>  
+  AAC <br>  
+  ABC <br>  
+  BBC <br>  
+  ACC <br>  
+  BCC <br>  
+  CCC <br><br><br><br><br><br><br><br><br><br><br>  
+   AAA <br>  
+  AAB <br>  
+  ABB <br>  
+  BBB <br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>  
+   AAA <br>  
+  <br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>  
+    
+  }  
+  
+  Another important observation here is that for a genotype <math>(a_1, a_2, .. a_P)</math>, The <math>F(P, a_p)</math>th genotype to <math>(a_1, a_2, .. a_P)</math> all end with <math>a_p</math>, this does not help distinguish the order, so we need to only examine the genotype <math>(a_1,..., a_{P1})</math>. The sub genotype <math>(a_1,..., a_{P1})</math> is also ordered else <math>(a_1,..., a_{P})</math> is not ordered.  
+  
+  The nested genotype sequence is in red. The blue sequence shows the sequence of genotypes enumerated without involving <math>a_P</math>.  
+  
+  { class="wikitable"  
+    
+  ! scope="col" Index  
+  ! scope="col" A=4,P=3  
+  ! scope="col" A=4,P=2  
+  ! scope="col" A=4,P=1  
+    
+  0 <br>  
+  1 <br>  
+  2 <br>  
+  3 <br>  
+  4 <br>  
+  5 <br>  
+  6 <br>  
+  7 <br>  
+  8 <br>  
+  9 <br>  
+  10 <br>  
+  11 <br>  
+  12 <br>  
+  13 <br>  
+  14 <br>  
+  15 <br>  
+  16 <br>  
+  17 <br>  
+  18 <br>  
+  19 <br>  
+   <span style="color:#FF0000">AAA</span><br>  
+  <span style="color:#FF0000">AAB</span> <br>  
+  <span style="color:#FF0000">ABB</span> <br>  
+  <span style="color:#FF0000">BBB</span> <br>  
+  <span style="color:#FF0000">AAC</span> <br>  
+  <span style="color:#FF0000">ABC</span> <br>  
+  <span style="color:#FF0000">BBC</span> <br>  
+  <span style="color:#FF0000">ACC</span> <br>  
+  <span style="color:#FF0000">BCC</span> <br>  
+  <span style="color:#FF0000">CCC</span> <br>  
+  <span style="color:#FF0000">AAD</span> <br>  
+  <span style="color:#FF0000">ABD</span> <br>  
+  <span style="color:#FF0000">BBD</span> <br>  
+  <span style="color:#FF0000">ACD</span> <br>  
+  <span style="color:#FF0000">BCD</span> <br>  
+  <span style="color:#FF0000">CCD</span> <br>  
+  <span style="color:#FF0000">ADD</span> <br>  
+  <span style="color:#FF0000">BDD</span> <br>  
+  <span style="color:#FF0000">CDD</span> <br>  
+  <span style="color:#FF0000">DDD</span> <br>  
+  <span style="color:#0000FF">AAA</span><br>  
+  <span style="color:#0000FF">AAB</span> <br>  
+  <span style="color:#0000FF">ABB</span> <br>  
+  <span style="color:#0000FF">BBB</span> <br>  
+  <span style="color:#0000FF">AAC</span> <br>  
+  <span style="color:#0000FF">ABC</span> <br>  
+  <span style="color:#0000FF">BBC</span> <br>  
+  <span style="color:#0000FF">ACC</span> <br>  
+  <span style="color:#0000FF">BCC</span> <br>  
+  <span style="color:#0000FF">CCC</span> <br>  
+  <span style="color:#FF0000">AA</span>D <br>  
+  <span style="color:#FF0000">AB</span>D <br>  
+  <span style="color:#FF0000">BB</span>D <br>  
+  <span style="color:#FF0000">AC</span>D <br>  
+  <span style="color:#FF0000">BC</span>D <br>  
+  <span style="color:#FF0000">CC</span>D <br>  
+  <span style="color:#FF0000">AD</span>D <br>  
+  <span style="color:#FF0000">BD</span>D <br>  
+  <span style="color:#FF0000">CD</span>D <br>  
+  <span style="color:#FF0000">DD</span>D <br>  
+   AAA<br>  
+  AAB <br>  
+  ABB <br>  
+  BBB <br>  
+  AAC <br>  
+  ABC <br>  
+  BBC <br>  
+  ACC <br>  
+  BCC <br>  
+  CCC <br>  
+  <span style="color:#0000FF">AA</span>D <br>  
+  <span style="color:#0000FF">AB</span>D <br>  
+  <span style="color:#0000FF">BB</span>D <br>  
+  <span style="color:#0000FF">AC</span>D <br>  
+  <span style="color:#0000FF">BC</span>D <br>  
+  <span style="color:#0000FF">CC</span>D <br>  
+  <span style="color:#FF0000">A</span>DD <br>  
+  <span style="color:#FF0000">B</span>DD <br>  
+  <span style="color:#FF0000">C</span>DD <br>  
+  <span style="color:#FF0000">D</span>DD <br>  
+    
+  }  
+  
+  The above 2 observations are the key to breaking down the enumeration recursively.  
+  
+  == Derivation ==  
+  
+  <math>a_1, ... a_P</math> is ordered and indexed 0 to A1. <br> <br>  
<math>  <math>  
\begin{align}  \begin{align}  
−  +  G(a_1,.. , a_P) &= \ \{\text{genotypes of } a_P \text{ alleles for ploidy P}\} \ + G(a_1,.. , a_{P1})\\  
+  &= F(P, a_P) + G(a_1,..,a_{P1}) \\  
+  &= F(P, a_P) + F(P1, a_{P1}) + G(a_1,..,a_{P2}) \\  
+  &= F(P, a_P) + F(P1, a_{P1}) + ... + F(1,a_1) \\  
+  &= \sum_{k=1}^P F(k, a_k) \\  
+  &= \sum_{k=1}^P \binom{k+a_k1}{a_k1}  
\end{align}  \end{align}  
</math>  </math>  
−  +  
−  +  <br>  
−  +  
+  This algorithm is demonstrated in the following table to obtain the index of the genotype C/C/D.  
+  
+  <br>  
+  
+  { class="wikitable"  
    
−    +  ! scope="col" Index 
−    +  ! scope="col" Iteration 0 
−    +  ! scope="col" 
−    +  ! scope="col" Iteration 1 
−    +  ! scope="col" 
+  ! scope="col" Iteration 2  
    
−   #  +  0 <br> 
−     +  1 <br> 
−    +  2 <br> 
−    +  3 <br> 
−    +  4 <br> 
+  5 <br>  
+  6 <br>  
+  7 <br>  
+  8 <br>  
+  9 <br>  
+  10 <br>  
+  11 <br>  
+  12 <br>  
+  13 <br>  
+  14 <br>  
+  <span style="color:#FF0000">15</span><br>  
+   AAA <br>  
+  AAB <br>  
+  ABB <br>  
+  BBB <br>  
+  AAC <br>  
+  ABC <br>  
+  BBC <br>  
+  ACC <br>  
+  BCC <br>  
+  CCC <br>  
+  AAD <br>  
+  ABD <br>  
+  BBD <br>  
+  ACD <br>  
+  BCD <br>  
+  CCD  
+   AAA <br>  
+  AAB <br>  
+  ABB <br>  
+  BBB <br>  
+  AAC <br>  
+  ABC <br>  
+  BBC <br>  
+  ACC <br>  
+  BCC <br>  
+  CCC <br><br><br><br><br><br><br>  
+   <br><br><br><br><br><br> <br><br><br><br>  
+  AA <br>  
+  AB <br>  
+  BB <br>  
+  AC <br>  
+  BC <br>  
+  CC <br>  
+   <br><br><br><br><br><br> <br><br><br><br>  
+  AA <br>  
+  AB <br>  
+  BB <br>  
+  <br><br><br>  
+   <br><br><br><br><br><br> <br><br><br><br><br><br><br>  
+  A <br>  
+  B <br>  
+  C <br>  
+    
+  Function call  
+  G(CCD)  
+  F(3,3)  
+  G(CC)  
+  F(2,2)  
+  G(C) = F(1, 2)  
+    
+  value returned  
+    
+  10  
+    
+  3  
+  2  
+  index=10+3+2=<span style="color:#FF0000">15</span> (QED!)  
    
−  
−  
−  
−  
−  
}  }  
+  
+  = Algorithm for enumerating the genotypes given ploidy and alleles =  
+  
+  The below code is for enumerating genotypes and can be used to test the above equations.  
+  
+  uint32_t no = 0 // some global variable  
+  void print_genotypes(uint32_t A, uint32_t P, std::string genotype)  
+  {  
+  if (genotype.size()==P)  
+  {  
+  std::cerr << no << ") " << genotype << "\n";  
+  ++no;  
+  }  
+  else  
+  {  
+  for (uint32_t a=0; a<A; ++a)  
+  {  
+  std::string s(1,(char)(a+65));  
+  s.append(genotype);  
+  print_genotypes(a+1, P, s);  
+  }  
+  }  
+  }  
+  
+  = Acknowledgement =  
+  
+  To [mailto:pd3@sanger.ac.uk Petr Danecek] for double checking this.  
+  
+  = Maintained by =  
+  
+  This page is maintained by [mailto:atks@umich.edu Adrian]. 
Latest revision as of 21:07, 18 October 2017
Contents
 1 Introduction
 2 Motivation
 3 The number of genotypes given a ploidy and alleles
 4 Getting the index of a genotype in an enumerated list given a ploidy and alleles
 5 Getting the genotypes from a genotype index and a given ploidy
 6 Simple cases
 7 Derivation for counting the number of genotypes
 8 Derivation for getting the index of a genotype in an enumerated list
 9 Algorithm for enumerating the genotypes given ploidy and alleles
 10 Acknowledgement
 11 Maintained by
Introduction
The VCF format encodes genotypes by the index of the enumeration of genotypes given ploidy number and alleles. This allows for direct access to a value associated with a genotype within an array when one works with genotype likelihoods.
Motivation
Plants species exhibit a diverse number of ploidy, for example, the strawberry is an octoploid and the pear is a triploid.
Copy number variations in somatic variant calling also leads to variable ploidy to consider when genotyping a locus.
While there are explicit functions that could be googled for handling haploid and diploid cases. It seems difficult to find the closed forms for the general case. This wiki fills in that need.
The number of genotypes given a ploidy and alleles
where P is the ploidy number and A is the number of alleles.
Getting the index of a genotype in an enumerated list given a ploidy and alleles
where P is the number of ploidy, , .. are the alleles in numeric encoding (0 to A1) and are ordered (e.g. AB and ABCCCC are ordered but ACB is not ordered).
This is well defined because:
Because may be 0, we will see cases of when . This is alright because of the definition of which defines this case as 0. But to make it more sensible, we can define the function equivalently as:
So when , the binomial coefficient reads as which equals 0 since there are 0 ways to choose k items from k1 items.
Getting the genotypes from a genotype index and a given ploidy
The genotype index is computed by a summation of a series which is monotonically decreasing. This allows you to compute the inverse function from index to the ordered genotypes by using a "water rapids algorithm" with decreasing height of each mini water fall. std::vector<int32_t> bcf_ip2g(int32_t genotype_index, uint32_t no_ploidy) { std::vector<int32_t> genotype(no_ploidy, 0); int32_t pth = no_ploidy; int32_t max_allele_index = genotype_index; int32_t leftover_genotype_index = genotype_index; while (pth>0) { for (int32_t allele_index=0; allele_index <= max_allele_index; ++allele_index) { int32_t i = choose(pth+allele_index1, pth); if (i>=leftover_genotype_index  allele_index==max_allele_index) { if (i>leftover_genotype_index) allele_index; leftover_genotype_index = choose(pth+allele_index1, pth); pth; max_allele_index = allele_index; genotype[pth] = allele_index; break; } } } return genotype; }
todo:: describe in a human understandable fashion.
Simple cases
Ploidy  Alleles  No. of Genotypes  Index 

1  A  
2  A 
Derivation for counting the number of genotypes
There must always be P observed alleles and there can only be at most A alleles. This can be modeled by P+A1 points where you choose A1 points to be dividers that separate the alleles. Thus the number of ways you can observe this is which is equivalent to .
Derivation for getting the index of a genotype in an enumerated list
Observation of nested patterns
An important observation here is that for the enumeration of A alleles for a given P ploidy, the enumeration of A1 alleles for P ploidy is a subsequence.
Index  A=4,P=3  A=3,P=3  A=2,P=3  A=1,P=3 

0 1 
AAA AAB 
AAA AAB 
AAA AAB 
AAA

Another important observation here is that for a genotype , The th genotype to all end with , this does not help distinguish the order, so we need to only examine the genotype . The sub genotype is also ordered else is not ordered.
The nested genotype sequence is in red. The blue sequence shows the sequence of genotypes enumerated without involving .
Index  A=4,P=3  A=4,P=2  A=4,P=1 

0 1 
AAA AAB 
AAA AAB 
AAA AAB 
The above 2 observations are the key to breaking down the enumeration recursively.
Derivation
is ordered and indexed 0 to A1.
This algorithm is demonstrated in the following table to obtain the index of the genotype C/C/D.
Index  Iteration 0  Iteration 1  Iteration 2  

0 1 
AAA AAB 
AAA AAB 
AA 
AA 
A  
Function call  G(CCD)  F(3,3)  G(CC)  F(2,2)  G(C) = F(1, 2)  
value returned  10  3  2  index=10+3+2=15 (QED!) 
Algorithm for enumerating the genotypes given ploidy and alleles
The below code is for enumerating genotypes and can be used to test the above equations.
uint32_t no = 0 // some global variable void print_genotypes(uint32_t A, uint32_t P, std::string genotype) { if (genotype.size()==P) { std::cerr << no << ") " << genotype << "\n"; ++no; } else { for (uint32_t a=0; a<A; ++a) { std::string s(1,(char)(a+65)); s.append(genotype); print_genotypes(a+1, P, s); } } }
Acknowledgement
To Petr Danecek for double checking this.
Maintained by
This page is maintained by Adrian.