Merge pull request #1831 from dib-lab/khmer_saghi

Saghi cms unique kmers

Makefile

@@ -184,7 +184,7 @@ cppcheck-long: FORCE
 
 ## pep8        : check Python code style
 pep8: $(PYSOURCES) $(wildcard tests/*.py)
-	pep8 setup.py khmer/*.py scripts/*.py tests/*.py oxli/*.py
+	pep8 setup.py khmer/*.py scripts/*.py tests/*.py oxli/*.py examples/python-api/*.py
 
 pep8_report.txt: $(PYSOURCES) $(wildcard tests/*.py)
 	pep8 setup.py khmer/ scripts/ tests/ oxli/ \
@@ -389,6 +389,7 @@ py-demos: sharedobj
 	python examples/python-api/exact-counting.py
 	python examples/python-api/bloom.py
 	python examples/python-api/consume.py examples/c++-api/reads.fastq
+	python examples/python-api/mask.py
 
 COMMIT ?= $(shell git rev-parse HEAD)
 docker-container:

examples/python-api/bloom.py

@@ -11,7 +11,7 @@
 num_tables = 4
 
 bloomfilter = khmer.Nodetable(ksize, target_table_size, num_tables)
-bloomfilter.consume('GCTGCACCGATGTACGCAAAGCTATTTAAAACCATAACTATTCTCACTTA');
+bloomfilter.consume('GCTGCACCGATGTACGCAAAGCTATTTAAAACCATAACTATTCTCACTTA')
 
 print('count for "GCTGCACCGATGTACGCAAAG" is',
       bloomfilter.get('GCTGCACCGATGTACGCAAAG'))

examples/python-api/exact-counting.py

@@ -26,7 +26,7 @@
 cg.count('ATGGCA')
 cg.count('ACATGG')
 cg.count('AAAAAA')
-cg.count('TTTTTT') # this will be counted towards AAAAAA
+cg.count('TTTTTT')  # this will be counted towards AAAAAA
 
 # Show all >0 k-mer abundances from the table
 for i in range(nkmers):

examples/python-api/mask.py

@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+# The following program uses khmer to
+# find unique kmers to only one sequence.
+import khmer
+d1 = "ATGTACGGGCATTACGATTACCGATGTAG"
+d2 = "ATGACCAAACTCATTACGATTAGATATAG"
+ksize = 5
+target_table_size = 5e5
+num_tables = 4
+bf = khmer.Nodetable(ksize, target_table_size, num_tables)
+bf.consume(d1)
+cms = khmer.Counttable(ksize, target_table_size, num_tables)
+for kmer in cms.get_kmers(d2):
+    if bf.get(kmer) == 0:
+        cms.consume(kmer)
+
+# If kmer is in both sequences it should not be in cms but in bf
+assert cms.get('CATTA') == 0
+assert bf.get('CATTA') > 0
+# If kmer is in d1 but not d2 it should not be in cms but be in bf
+assert cms.get('ATGTA') == 0
+assert bf.get('ATGTA') > 0
+# If kmer is in d2 but not d1 it should be in cms and not in bf
+assert cms.get('TATAG') > 0
+assert bf.get('TATAG') == 0