fix IAF computation for all samples (#563) (#564)

* fix IAF computation for all samples * fix Mamba build Co-authored-by: Adam Wenocur <[email protected]>
TileDB-Inc · Aug 16, 2023 · 51c9bdb · 51c9bdb
1 parent a87964b
commit 51c9bdb
Show file tree

Hide file tree

Showing 11 changed files with 18 additions and 12 deletions.
diff --git a/apis/python/tests/test_tiledbvcf.py b/apis/python/tests/test_tiledbvcf.py
@@ -1166,7 +1166,7 @@ def test_ingest_with_stats(tmp_path):
     ).bool()
     assert (
         data_frame[data_frame["sample_name"] == "second"]["info_TILEDB_IAF"].iloc[0][0]
-        == 0.0625
+        == 0.9375
     )
     data_frame = ds.read(
         samples=sample_names,
@@ -1177,7 +1177,7 @@ def test_ingest_with_stats(tmp_path):
         data_frame[
             (data_frame["sample_name"] == "second") & (data_frame["pos_start"] == 4)
         ]["info_TILEDB_IAF"].iloc[0][0]
-        == 0.125
+        == 0.9375
     )
 
 

diff --git a/ci/azure-linux_mac.yml b/ci/azure-linux_mac.yml
@@ -152,6 +152,12 @@ steps:
           sudo chown -R $USER $CONDA
         fi
 
+        if [[ "$AGENT_OS" == "Darwin" ]]; then
+           conda install -y libarchive=3.6.2=h6d8d9f1_0 -n base -c conda-forge
+        else
+           conda install -y libarchive -n base -c conda-forge
+        fi
+
         # Install mamba for faster environment creation
         conda install -y mamba -n base -c conda-forge
 

diff --git a/libtiledbvcf/src/read/reader.cc b/libtiledbvcf/src/read/reader.cc
@@ -1298,8 +1298,8 @@ bool Reader::process_query_results_v4() {
         // If the allele is in GT, consider it in the pass computation
         // TODO: when supporting greater than diploid organisms, expand the
         // following boolean statement into a loop
-        if ((gt.size() > 0 && allele_index == gt[0]) ||
-            (gt.size() > 1 && allele_index == gt[1])) {
+        if (!is_ref && ((gt.size() > 0 && allele_index == gt[0]) ||
+                        (gt.size() > 1 && allele_index == gt[1]))) {
           pass = pass || allele_passes;
         } else {
           LOG_TRACE("  ignore allele {} not in GT", allele_index);

diff --git a/libtiledbvcf/src/stats/variant_stats_reader.h b/libtiledbvcf/src/stats/variant_stats_reader.h
@@ -137,7 +137,7 @@ class AFMap {
     if (next_allele == pos_map.second.end()) {
       return 0;
     }
-    return 1.0 * next_allele->second / num_samples;
+    return 1.0 * next_allele->second / num_samples / 2;
   }
 
   /**

diff --git a/libtiledbvcf/test/inputs/stats/eighth.vcf b/libtiledbvcf/test/inputs/stats/eighth.vcf
@@ -10,4 +10,4 @@
 chr1	1	.	T	C,<NON_REF>	829.77	.	ExcessHet=3.0103	GT	0/1
 chr1	2	.	G	GTTTA,T,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	2/2
 chr1	3	.	C	A,T	10.032	.	ExcessHet=8.34	GT	0/1
-chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	1/1
+chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	0/0
diff --git a/libtiledbvcf/test/inputs/stats/fifth.vcf b/libtiledbvcf/test/inputs/stats/fifth.vcf
@@ -10,4 +10,4 @@
 chr1	1	.	T	C,<NON_REF>	829.77	.	ExcessHet=3.0103	GT:AD:AB	0/1:1:1.03493
 chr1	2	.	G	GTTTA,T,<NON_REF>	343.73	.	ExcessHet=2.8532	GT:AD:AB	2/2:3:5
 chr1	3	.	C	A,T	10.032	.	ExcessHet=8.34	GT:AD:AB	2/2:20:3.5
-chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT:AD:AB	1/1:12:83.2
+chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT:AD:AB	0/0:12:83.2
diff --git a/libtiledbvcf/test/inputs/stats/first.vcf b/libtiledbvcf/test/inputs/stats/first.vcf
@@ -10,7 +10,7 @@
 chr1	1	.	T	C,<NON_REF>	829.77	.	ExcessHet=3.0103	GT	0/1
 chr1	2	.	G	GTTTA,T,<NON_REF>	1042.73	.	ExcessHet=4.8532	GT	0/1
 chr1	3	.	C	A,G	10.032	.	ExcessHet=2.0134	GT	2/2
-chr1	4	.	G	GTTTA,<NON_REF>	1042.73	.	ExcessHet=4.8532	GT	1/1
+chr1	4	.	G	GTTTA,<NON_REF>	1042.73	.	ExcessHet=4.8532	GT	0/0
 chr2	1	.	G	GTTTA	1042.73	.	ExcessHet=4.8532	GT	./1
 chr2	1	.	G	GTTTA	1042.73	.	ExcessHet=4.8532	GT	1/.
 chr2	3	.	G	GTTTA	1042.73	.	ExcessHet=4.8532	GT	1/1

diff --git a/libtiledbvcf/test/inputs/stats/fourth.vcf b/libtiledbvcf/test/inputs/stats/fourth.vcf
@@ -10,4 +10,4 @@
 chr1	1	.	T	C,<NON_REF>	829.77	.	ExcessHet=3.0103	GT	0/1
 chr1	2	.	G	GTTTA,T,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	0/1
 chr1	3	.	C	A,T	10.032	.	ExcessHet=8.34	GT	2/2
-chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	1/1
+chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	0/0
diff --git a/libtiledbvcf/test/inputs/stats/seventh.vcf b/libtiledbvcf/test/inputs/stats/seventh.vcf
@@ -10,4 +10,4 @@
 chr1	1	.	T	C,<NON_REF>	829.77	.	ExcessHet=3.0103	GT	0/1
 chr1	2	.	G	GTTTA,T,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	0/1
 chr1	3	.	C	A,T	10.032	.	ExcessHet=8.34	GT	0/1
-chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	1/1
+chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	0/0
diff --git a/libtiledbvcf/test/inputs/stats/sixth.vcf b/libtiledbvcf/test/inputs/stats/sixth.vcf
@@ -10,4 +10,4 @@
 chr1	1	.	T	C,<NON_REF>	829.77	.	ExcessHet=3.0103	GT	0/1
 chr1	2	.	G	GTTTA,T,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	2/2
 chr1	3	.	C	A,T	10.032	.	ExcessHet=8.34	GT	0/1
-chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	1/1
+chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	0/0
diff --git a/libtiledbvcf/test/inputs/stats/third.vcf b/libtiledbvcf/test/inputs/stats/third.vcf
@@ -10,4 +10,4 @@
 chr1	1	.	T	C,<NON_REF>	829.77	.	ExcessHet=3.0103	GT	0/1
 chr1	2	.	G	GTTTA,T,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	0/1
 chr1	3	.	C	A,T	10.032	.	ExcessHet=8.34	GT	0/1
-chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	1/1
+chr1	4	.	G	GTTTA,<NON_REF>	343.73	.	ExcessHet=2.8532	GT	0/0