Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[Haplostats] EM is not ignoring missing data #41

Open
sjmack opened this issue Aug 24, 2017 · 1 comment
Open

[Haplostats] EM is not ignoring missing data #41

sjmack opened this issue Aug 24, 2017 · 1 comment
Assignees
@rsingle
Labels
bug Something isn't working, should not be used for new features: use "enhancement" for those
Milestone
Add Haplostats of...

Comments

@sjmack
Copy link
Collaborator

sjmack commented Aug 24, 2017

Using the same .ini and .pop files as in Issue #40 , it looks like [Haplostats] is including missing data (****) in its haploype estimation.

Here is an example of [Emhaplofreq] results:

Haplotype frequency est. for loci: A:B:DRB1
-------------------------------------------
Number of individuals: 47 (before-filtering)
Number of individuals: 45 (after-filtering)
Unique phenotypes: 45
Unique genotypes: 113
Number of haplotypes: 188
Loglikelihood under linkage equilibrium [ln(L_0)]: -359.216808
Loglikelihood obtained via the EM algorithm [ln(L_1)]: -340.336732
Number of iterations before convergence: 44

Haplotypes sorted by name             | Haplotypes sorted by frequency        
haplotype            frequency# copies| haplotype            frequency# copies
0101~1301~0402       0.02222  2.0     | 0201~1401~0402       0.04444  4.0     
0101~1301~1101       0.01111  1.0     | 03012~1401~0407      0.03333  3.0     
0101~1401~0901       0.01111  1.0     | 03012~1301~0402      0.03333  3.0     
0101~1520~1602       0.01111  1.0     | 0201~1401~0802       0.03333  3.0     
0101~18012~0407      0.01111  1.0     | 3204~1401~0802       0.02222  2.0     
0101~39021~0404      0.01111  1.0     | 0201~1401~1101       0.02222  2.0     
0101~39021~1602      0.01111  1.0     | 0101~4005~0802       0.02222  2.0     
0101~4005~0802       0.02222  2.0     | 03012~39021~0402     0.02222  2.0     
0101~67011~0802      0.01111  1.0     | 0201~1301~1602       0.02222  2.0     
0101~8101~1602       0.01111  1.0     | 0218~1401~0404       0.02222  2.0     
0201~1301~1602       0.02222  2.0     | 0210~51013~1602      0.02222  2.0     
0201~1401~0402       0.04444  4.0     | 0218~1401~0802       0.02222  2.0     
0201~1401~0404       0.01111  1.0     | 0101~1301~0402       0.02222  2.0     
0201~1401~0407       0.01111  1.0     | 2501~4005~0802       0.02222  2.0     
0201~1401~0802       0.03333  3.0     | 2501~1301~0802       0.02222  2.0     
0201~1401~0901       0.01111  1.0     | 2501~39021~0402      0.02222  2.0     
0201~1401~1101       0.02222  2.0     | 3204~39021~0802      0.02222  2.0     
0201~18012~0407      0.01111  1.0     | 03012~1520~0802      0.02222  2.0     
0201~39021~0901      0.01111  1.0     | 0201~1401~0901       0.01111  1.0     
0201~4005~0402       0.01111  1.0     | 0101~39021~1602      0.01111  1.0     
0201~51013~0802      0.01111  1.0     | 0210~1301~0407       0.01111  1.0     
0201~5703~1101       0.01111  1.0     | 0218~35091~0901      0.01111  1.0     
0210~1301~0402       0.01111  1.0     | 0210~78021~0901      0.01111  1.0     
0210~1301~0407       0.01111  1.0     | 3204~1301~0901       0.01111  1.0     
0210~1401~0404       0.01111  1.0     | 3204~51013~1101      0.01111  1.0     
0210~35091~0402      0.01111  1.0     | 2501~1401~0404       0.01111  1.0     
0210~39021~0802      0.01111  1.0     | 6814~1520~0407       0.01111  1.0     
0210~4005~0802       0.01111  1.0     | 0201~18012~0407      0.01111  1.0     
0210~51013~1602      0.02222  2.0     | 0201~4005~0402       0.01111  1.0     
0210~78021~0901      0.01111  1.0     | 0101~8101~1602       0.01111  1.0     
0218~1301~0402       0.01111  1.0     | 0210~1301~0402       0.01111  1.0     
0218~1401~0402       0.01111  1.0     | 6901~1301~0402       0.01111  1.0     
0218~1401~0404       0.02222  2.0     | 0210~4005~0802       0.01111  1.0     
0218~1401~0407       0.01111  1.0     | 03012~4005~0802      0.01111  1.0     
0218~1401~0802       0.02222  2.0     | 0101~1401~0901       0.01111  1.0     
0218~35091~0901      0.01111  1.0     | 0101~1301~1101       0.01111  1.0     
0218~4005~1101       0.01111  1.0     | 0101~1520~1602       0.01111  1.0     
03012~1301~0402      0.03333  3.0     | 0218~1301~0402       0.01111  1.0     
03012~1401~0407      0.03333  3.0     | 0201~51013~0802      0.01111  1.0     
03012~1401~0802      0.01111  1.0     | 0201~39021~0901      0.01111  1.0     
03012~1520~0802      0.02222  2.0     | 2501~1520~0404       0.01111  1.0     
03012~39021~0402     0.02222  2.0     | 6814~5703~1602       0.01111  1.0     
03012~39021~0802     0.01111  1.0     | 0201~1401~0404       0.01111  1.0     
03012~4005~0802      0.01111  1.0     | 0201~1401~0407       0.01111  1.0     
03012~51013~0901     0.01111  1.0     | 3204~1301~1602       0.01111  1.0     
2501~1301~0802       0.02222  2.0     | 0218~1401~0407       0.01111  1.0     
2501~1401~0404       0.01111  1.0     | 0218~4005~1101       0.01111  1.0     
2501~1520~0404       0.01111  1.0     | 03012~1401~0802      0.01111  1.0     
2501~39021~0402      0.02222  2.0     | 03012~39021~0802     0.01111  1.0     
2501~4005~0802       0.02222  2.0     | 0101~39021~0404      0.01111  1.0     
2501~51013~0901      0.01111  1.0     | 0210~1401~0404       0.01111  1.0     
2501~5703~0802       0.01111  1.0     | 6814~4005~0407       0.01111  1.0     
2501~8101~0802       0.01111  1.0     | 0101~18012~0407      0.01111  1.0     
3204~1301~0901       0.01111  1.0     | 6901~1520~0402       0.01111  1.0     
3204~1301~1602       0.01111  1.0     | 0210~35091~0402      0.01111  1.0     
3204~1401~0802       0.02222  2.0     | 0210~39021~0802      0.01111  1.0     
3204~1520~0802       0.01111  1.0     | 03012~51013~0901     0.01111  1.0     
3204~39021~0802      0.02222  2.0     | 6901~67011~1602      0.01111  1.0     
3204~51013~1101      0.01111  1.0     | 0218~1401~0402       0.01111  1.0     
6814~1520~0407       0.01111  1.0     | 0101~67011~0802      0.01111  1.0     
6814~4005~0407       0.01111  1.0     | 2501~8101~0802       0.01111  1.0     
6814~5703~1602       0.01111  1.0     | 7403~39021~0407      0.01111  1.0     
6901~1301~0402       0.01111  1.0     | 6901~1301~0404       0.01111  1.0     
6901~1301~0404       0.01111  1.0     | 2501~5703~0802       0.01111  1.0     
6901~1520~0402       0.01111  1.0     | 2501~51013~0901      0.01111  1.0     
6901~67011~1602      0.01111  1.0     | 0201~5703~1101       0.01111  1.0     
7403~39021~0407      0.01111  1.0     | 3204~1520~0802       0.01111  1.0

And here is an example of [Haplostats] results:

II. Multi-locus Analyses [haplo-stats]
======================================

Haplotype / linkage disequilibrium (LD) statistics
__________________________________________________

Pairwise LD estimates
---------------------
Locus pair            D      D'        Wn   ln(L_1)   ln(L_0)         S # permu p-value
A:B                   *0.516108  0.390722   -314.17                 NaN       - -      


Haplotype frequency est. for loci: A:B
--------------------------------------
Unique genotypes: 50
Number of haplotypes: 52
Loglikelihood under linkage equilibrium [ln(L_0)]: 
Loglikelihood obtained via the EM algorithm [ln(L_1)]: -314.173
Number of iterations before convergence: 

Haplotypes sorted by name             | Haplotypes sorted by frequency        
haplotype            frequency# copies| haplotype            frequency# copies
****~1301            0.0106382        | 0201~1401            0.1271851        
****~1401            0.0106382        | 0218~1401            0.0638297        
****~18012           0.0106382        | 03012~1401           0.0319148        
****~78021           0.0106382        | 0101~1301            0.0319148        
0101~1301            0.0319148        | 0210~1301            0.0319148        
0101~1401            0.0111127        | 2501~1301            0.0319148        
0101~1520            0.0212765        | 03012~1520           0.0319148        
0101~35091           0.0106382        | 03012~51013          0.0319148        
0101~39021           0.0208021        | 2501~39021           0.0319146        
0101~51013           0.0106382        | 3204~39021           0.0212768        
0101~8101            0.0212765        | 0101~1520            0.0212765        
0201~1301            0.0212765        | 0101~8101            0.0212765        
0201~1401            0.1271851        | 0201~1301            0.0212765        
0201~18012           0.0106382        | 0201~4005            0.0212765        
0201~39021           0.0111127        | 0210~4005            0.0212765        
0201~4005            0.0212765        | 0218~4005            0.0212765        
0201~5703            0.0106382        | 2501~51013           0.0212765        
0210~1301            0.0319148        | 03012~1301           0.0212765        
0210~1401            0.0212765        | 03012~39021          0.0212765        
0210~35091           0.0106382        | 3204~1401            0.0212765        
0210~39021           0.0106382        | 6901~1301            0.0212765        
0210~4005            0.0212765        | 0210~1401            0.0212765        
0210~51013           8.1821857        | 0101~39021           0.0208021        
0218~1301            0.0106382        | 0101~1401            0.0111127        
0218~1401            0.0638297        | 0201~39021           0.0111127        
0218~4005            0.0212765        | 0101~35091           0.0106382        
03012~1301           0.0212765        | 0101~51013           0.0106382        
03012~1401           0.0319148        | 0201~5703            0.0106382        
03012~1520           0.0319148        | 0201~18012           0.0106382        
03012~39021          0.0212765        | 0210~35091           0.0106382        
03012~4005           0.0106382        | 0210~39021           0.0106382        
03012~51013          0.0319148        | 0218~1301            0.0106382        
2501~1301            0.0319148        | 2501~1401            0.0106382        
2501~1401            0.0106382        | 2501~5703            0.0106382        
2501~39021           0.0319146        | 2501~67011           0.0106382        
2501~4005            2.9233526        | 03012~4005           0.0106382        
2501~51013           0.0212765        | 3204~1301            0.0106382        
2501~5703            0.0106382        | 3204~1520            0.0106382        
2501~67011           0.0106382        | 3204~78021           0.0106382        
3204~1301            0.0106382        | 6814~1520            0.0106382        
3204~1401            0.0212765        | 6814~4005            0.0106382        
3204~1520            0.0106382        | 6814~5703            0.0106382        
3204~39021           0.0212768        | 6901~18012           0.0106382        
3204~4005            0.0106380        | 6901~67011           0.0106382        
3204~78021           0.0106382        | 7403~39021           0.0106382        
6814~1520            0.0106382        | ****~1301            0.0106382        
6814~4005            0.0106382        | ****~1401            0.0106382        
6814~5703            0.0106382        | ****~18012           0.0106382        
6901~1301            0.0212765        | ****~78021           0.0106382        
6901~18012           0.0106382        | 3204~4005            0.0106380        
6901~67011           0.0106382        | 2501~4005            2.9233526        
7403~39021           0.0106382        | 0210~51013           8.1821857
@sjmack sjmack added the bug Something isn't working, should not be used for new features: use "enhancement" for those label Aug 24, 2017
@sjmack sjmack mentioned this issue Dec 5, 2017
Closed
@sjmack sjmack mentioned this issue Feb 3, 2018
Closed
@alexlancaster
Copy link
Owner

@rsingle not sure if this is still a problem with current haplo-stats in main branch? perhaps we don't apply the missing data filter in the same way we do with the emhaplofreq implementation? I'll take a look and see if I can still reproduce this.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@rsingle rsingle

Labels
bug Something isn't working, should not be used for new features: use "enhancement" for those
Projects
None yet
Milestone
Add Haplostats officially to release
Development

No branches or pull requests
3 participants
@alexlancaster @rsingle @sjmack