4
4
import lshmm as ls
5
5
import msprime
6
6
import numpy as np
7
+ import pytest
7
8
8
9
import tskit
9
10
@@ -1257,18 +1258,21 @@ def assertAllClose(self, A, B):
1257
1258
1258
1259
# Define a bunch of very small tree-sequences for testing a collection of
1259
1260
# parameters on
1261
+ @pytest .mark .skip (reason = "No plans to implement diploid LS HMM yet." )
1260
1262
def test_simple_n_10_no_recombination (self ):
1261
1263
ts = msprime .simulate (
1262
1264
10 , recombination_rate = 0 , mutation_rate = 0.5 , random_seed = 42
1263
1265
)
1264
1266
assert ts .num_sites > 3
1265
1267
self .verify (ts )
1266
1268
1269
+ @pytest .mark .skip (reason = "No plans to implement diploid LS HMM yet." )
1267
1270
def test_simple_n_6 (self ):
1268
1271
ts = msprime .simulate (6 , recombination_rate = 2 , mutation_rate = 7 , random_seed = 42 )
1269
1272
assert ts .num_sites > 5
1270
1273
self .verify (ts )
1271
1274
1275
+ @pytest .mark .skip (reason = "No plans to implement diploid LS HMM yet." )
1272
1276
def test_simple_n_8_high_recombination (self ):
1273
1277
ts = msprime .simulate (8 , recombination_rate = 20 , mutation_rate = 5 , random_seed = 42 )
1274
1278
assert ts .num_trees > 15
@@ -1326,11 +1330,10 @@ def verify(self, ts):
1326
1330
ts_check , mapping = ts .simplify (
1327
1331
range (1 , n + 1 ), filter_sites = False , map_nodes = True
1328
1332
)
1333
+ H_check = ts_check .genotype_matrix ()
1329
1334
G_check = np .zeros ((m , n , n ))
1330
1335
for i in range (m ):
1331
- G_check [i , :, :] = np .add .outer (
1332
- ts_check .genotype_matrix ()[i , :], ts_check .genotype_matrix ()[i , :]
1333
- )
1336
+ G_check [i , :, :] = np .add .outer (H_check [i , :], H_check [i , :])
1334
1337
1335
1338
cm_d = ls_forward_tree (s [0 , :], ts_check , r , mu )
1336
1339
ll_tree = np .sum (np .log10 (cm_d .normalisation_factor ))
@@ -1345,12 +1348,16 @@ def verify(self, ts):
1345
1348
self .assertAllClose (ll_tree , ll_mirror_tree_dict )
1346
1349
1347
1350
# Ensure that the decoded matrices are the same
1351
+ flipped_H_check = np .flip (H_check , axis = 0 )
1352
+ flipped_s = np .flip (s , axis = 1 )
1353
+
1348
1354
F_mirror_matrix , c , ll = ls .forwards (
1349
- np .flip (G_check , axis = 0 ),
1350
- np .flip (s , axis = 1 ),
1351
- r_flip ,
1352
- p_mutation = np .flip (mu ),
1353
- scale_mutation_based_on_n_alleles = False ,
1355
+ flipped_H_check ,
1356
+ flipped_s ,
1357
+ ploidy = 2 ,
1358
+ prob_recombination = r_flip ,
1359
+ prob_mutation = np .flip (mu ),
1360
+ scale_mutation_rate = False ,
1354
1361
)
1355
1362
1356
1363
self .assertAllClose (F_mirror_matrix , cm_mirror .decode ())
@@ -1367,14 +1374,18 @@ def verify(self, ts):
1367
1374
ts_check , mapping = ts .simplify (
1368
1375
range (1 , n + 1 ), filter_sites = False , map_nodes = True
1369
1376
)
1377
+ H_check = ts_check .genotype_matrix ()
1370
1378
G_check = np .zeros ((m , n , n ))
1371
1379
for i in range (m ):
1372
- G_check [i , :, :] = np .add .outer (
1373
- ts_check .genotype_matrix ()[i , :], ts_check .genotype_matrix ()[i , :]
1374
- )
1380
+ G_check [i , :, :] = np .add .outer (H_check [i , :], H_check [i , :])
1375
1381
1376
1382
F , c , ll = ls .forwards (
1377
- G_check , s , r , p_mutation = mu , scale_mutation_based_on_n_alleles = False
1383
+ reference_panel = H_check ,
1384
+ query = s ,
1385
+ ploidy = 2 ,
1386
+ prob_recombination = r ,
1387
+ prob_mutation = mu ,
1388
+ scale_mutation_rate = False ,
1378
1389
)
1379
1390
cm_d = ls_forward_tree (s [0 , :], ts_check , r , mu )
1380
1391
self .assertAllClose (cm_d .decode (), F )
@@ -1393,22 +1404,27 @@ def verify(self, ts):
1393
1404
ts_check , mapping = ts .simplify (
1394
1405
range (1 , n + 1 ), filter_sites = False , map_nodes = True
1395
1406
)
1407
+ H_check = ts_check .genotype_matrix ()
1396
1408
G_check = np .zeros ((m , n , n ))
1397
1409
for i in range (m ):
1398
- G_check [i , :, :] = np .add .outer (
1399
- ts_check .genotype_matrix ()[i , :], ts_check .genotype_matrix ()[i , :]
1400
- )
1410
+ G_check [i , :, :] = np .add .outer (H_check [i , :], H_check [i , :])
1401
1411
1402
1412
F , c , ll = ls .forwards (
1403
- G_check , s , r , p_mutation = mu , scale_mutation_based_on_n_alleles = False
1413
+ reference_panel = H_check ,
1414
+ query = s ,
1415
+ ploidy = 2 ,
1416
+ prob_recombination = r ,
1417
+ prob_mutation = mu ,
1418
+ scale_mutation_rate = False ,
1404
1419
)
1405
1420
B = ls .backwards (
1406
- G_check ,
1407
- s ,
1408
- c ,
1409
- r ,
1410
- p_mutation = mu ,
1411
- scale_mutation_based_on_n_alleles = False ,
1421
+ reference_panel = H_check ,
1422
+ query = s ,
1423
+ ploidy = 2 ,
1424
+ normalisation_factor_from_forward = c ,
1425
+ prob_recombination = r ,
1426
+ prob_mutation = mu ,
1427
+ scale_mutation_rate = False ,
1412
1428
)
1413
1429
1414
1430
# Note, need to remove the first sample from the ts, and ensure that
@@ -1447,22 +1463,28 @@ def verify(self, ts):
1447
1463
ts_check , mapping = ts .simplify (
1448
1464
range (1 , n + 1 ), filter_sites = False , map_nodes = True
1449
1465
)
1466
+ H_check = ts_check .genotype_matrix ()
1450
1467
G_check = np .zeros ((m , n , n ))
1451
1468
for i in range (m ):
1452
- G_check [i , :, :] = np .add .outer (
1453
- ts_check .genotype_matrix ()[i , :], ts_check .genotype_matrix ()[i , :]
1454
- )
1469
+ G_check [i , :, :] = np .add .outer (H_check [i , :], H_check [i , :])
1455
1470
ts_check = ts .simplify (range (1 , n + 1 ), filter_sites = False )
1471
+
1456
1472
phased_path , ll = ls .viterbi (
1457
- G_check , s , r , p_mutation = mu , scale_mutation_based_on_n_alleles = False
1473
+ reference_panel = H_check ,
1474
+ query = s ,
1475
+ ploidy = 2 ,
1476
+ prob_recombination = r ,
1477
+ prob_mutation = mu ,
1478
+ scale_mutation_rate = False ,
1458
1479
)
1459
- path_ll_matrix = ls .path_ll (
1460
- G_check ,
1461
- s ,
1462
- phased_path ,
1463
- r ,
1464
- p_mutation = mu ,
1465
- scale_mutation_based_on_n_alleles = False ,
1480
+ path_ll_matrix = ls .path_loglik (
1481
+ reference_panel = H_check ,
1482
+ query = s ,
1483
+ ploidy = 2 ,
1484
+ path = phased_path ,
1485
+ prob_recombination = r ,
1486
+ prob_mutation = mu ,
1487
+ scale_mutation_rate = False ,
1466
1488
)
1467
1489
1468
1490
c_v = ls_viterbi_tree (s [0 , :], ts_check , r , mu )
@@ -1471,13 +1493,14 @@ def verify(self, ts):
1471
1493
# Attempt to get the path
1472
1494
path_tree_dict = c_v .traceback ()
1473
1495
# Work out the likelihood of the proposed path
1474
- path_ll_tree = ls .path_ll (
1475
- G_check ,
1476
- s ,
1477
- np .transpose (path_tree_dict ),
1478
- r ,
1479
- p_mutation = mu ,
1480
- scale_mutation_based_on_n_alleles = False ,
1496
+ path_ll_tree = ls .path_loglik (
1497
+ reference_panel = H_check ,
1498
+ query = s ,
1499
+ ploidy = 2 ,
1500
+ path = np .transpose (path_tree_dict ),
1501
+ prob_recombination = r ,
1502
+ prob_mutation = mu ,
1503
+ scale_mutation_rate = False ,
1481
1504
)
1482
1505
1483
1506
self .assertAllClose (ll , ll_tree )
0 commit comments