fix path and filter options in example code

example/read_mat.py

@@ -2,14 +2,15 @@
 
 import numpy as np
 from subprocess import run
+from pathlib import Path
 
 
 def read_ibdtools_mat(ibdtools_mat_fn):
     """read ibd matrix and turn it to a numpy array"""
 
     # decompress the matrix file
     flat_mat_fn = "tmp_ibdtools_mat_flat"
-    cmd = f"gzip -dc {ibdtools_mat_fn} > {flat_mat_fn}"
+    cmd = f"gzip -dc {str(ibdtools_mat_fn)} > {flat_mat_fn}"
     run(cmd, shell=True)
 
     # uint8: is_hap,                         offset 0
@@ -30,10 +31,10 @@ def read_ibdtools_mat(ibdtools_mat_fn):
     arr_sz = yy // 2  # each array item is uint16
     subpop_sz = np.uint64(xx // 4)  # each supop id is uint32
 
-    print(f"- is_hap      : {is_hap}")
-    print(f"- arr_sz      : {arr_sz}")
-    print(f"- d           : {d}")
-    print(f"- #subpop ids : {subpop_sz}")
+    # print(f"- is_hap      : {is_hap}")
+    # print(f"- arr_sz      : {arr_sz}")
+    # print(f"- d           : {d}")
+    # print(f"- #subpop ids : {subpop_sz}")
 
     # load subpop id
     if xx == 0:
@@ -63,5 +64,15 @@ def read_ibdtools_mat(ibdtools_mat_fn):
     return is_hap, subpop_ids, M
 
 
-res = read_ibdtools_mat("./gw.mat.mat")
-print(res, res[2].sum())
+path = list(Path().glob("gw.mat*mat"))[0]
+res = read_ibdtools_mat(path)
+M = res[2]
+print(
+    f"""
+    ibd matrix: dimension = {M.shape}, 
+    max element =  {M.max()}
+    min_element = {M.min()}
+    min_nonzero_element = {M[M>0].min()}
+    percentage of nonzero elements = {M[M>0].size / M.size}
+    """
+)

example/run_ibdtools.sh

@@ -2,19 +2,16 @@
 
 set -e
 
-# use msprime to simulate vcf and call ibd using hapibd. The output is used as input for ibdtools
-python3 ./simulate_data.py
-
 # use ibdtools to encode, split, deocde, sort, merge ibd and calculate chromosome-level total ibd matrix
 for i in {1..5}; do 
 	../build/ibdtools encode -i $i.ibd.gz -v $i.vcf.gz -g $i.map -c $i -o $i.eibd -m $i.meta -M 1
 	../build/ibdtools split -i $i.eibd -m $i.meta -W 2 -S 10 -o $i.split -M 1
 	../build/ibdtools decode -i $i.split1 -m $i.meta -o $i.ibd.gz -M 1
 	../build/ibdtools sort -i $i.split1 -o $i.sibd -M 1
 	../build/ibdtools merge -i $i.sibd -m $i.meta -o $i.mibd -M 1
-	../build/ibdtools matrix -i $i.mibd -m $i.meta -M 1 -T 2.5 -L 5 -o $i.mat
+	../build/ibdtools matrix -i $i.mibd -m $i.meta -M 1 -B 2.5 -o $i.mat
 done
 
 # collect all chromosome-wide total matrices and join them with comma
 MATRICES=`\ls {1..5}.mat*.mat | tr '\n' ','| sed s/,$//` 
-../build/ibdtools matrix -x $MATRICES -m 1.meta -T 2.5 -L 5 -o gw.mat -M 3
+../build/ibdtools matrix -x $MATRICES -m 1.meta -L 10 -o gw.mat -M 3

example/simulate_data.py

@@ -11,7 +11,7 @@
 chrlens = [(20 + i * 4) * bp_per_cm for i in range(1, 6)]
 
 for i, (chrno, chrlen) in enumerate(zip(chroms, chrlens)):
-    print(f"Simulating chr {chrno}")
+    print(f"Simulating chr {chrno} via msprime")
 
     # simulate ancestry
     ts = msprime.sim_ancestry(
@@ -32,11 +32,11 @@
             f, contig_id=f"{chrno}", individual_names=[f"ind{i:3d}" for i in range(50)]
         )
 
-    # add a low density region
+    # add a zero snp-density region
     cmd = f"""
         bgzip -f {chrno}.vcf; mv {chrno}.vcf.gz tmp_{chrno}.vcf.gz; bcftools index -f tmp_{chrno}.vcf.gz
         bcftools view -r {chrno}:1-5000000,{chrno}:8000000-{chrlen} tmp_{chrno}.vcf.gz -Oz -o {chrno}.vcf.gz
-        rm tmp_{chrno}.vcf.gz
+        rm tmp_{chrno}.vcf.gz*
     """
     run(cmd, shell=True)
 
@@ -60,8 +60,9 @@
             f.write(response.content)
 
     # run hapibd
+    print("\tCalling IBD via hap-ibd")
     cmd = f"java -Xmx1g -jar {hapibd_jar} gt={chrno}.vcf.gz map={chrno}.map out={chrno}"
-    run(cmd, shell=True)
+    _ = run(cmd, shell=True, capture_output=True)
 
     # clean up files
     pathlib.Path(f"{chrno}.log").unlink()