Repository: pachterlab/kb_python (language: python)

Issue:
Add kb extract command

Why are UMI counts in adata.h5ad and cells_x_genes.total.mtx different?
**Describe the issue**
Thank you for creating such a wonderful tool! I would like to ask you about the difference between `adata.h5ad` and `cells_x_genes.total.mtx`. I ran `kb count` and used an output file `adata.h5ad` for downstream processing by scanpy, but the UMI count stored in the h5ad file seems to be much smaller than I got from another tool, STARsolo. Next I used `cells_x_genes.total.mtx` which shows a very similar UMI count with STARsolo. I think it makes sense to use `cells_x_genes.total.mtx` instead of `adata.h5ad`, but I still don't understand what makes the difference between these files. Is this a bug?

**What is the exact command that was run?**
I ran kb in Snakemake workflow.
```
rule kb_count:
	wildcard_constraints:
		sample = "|".join([re.escape(str(x)) for x in samples])
	container:
		"docker://quay.io/biocontainers/kb-python:0.28.2--pyhdfd78af_2"
	input:
		kb_index = "kb_index/index.idx",
		kb_t2g = "kb_index/t2g.txt",
		kb_cdna_fa = "kb_index/cdna.fa",
		kb_intron_fa = "kb_index/intron.fa",
		kb_cdna_t2c = "kb_index/cdna_t2c.txt",
		kb_intron_t2c = "kb_index/intron_t2c.txt"
	output:
		adata_unfiltered_h5ad = "kb/{sample}/counts_unfiltered/adata.h5ad",
		adata_filtered_h5ad = "kb/{sample}/counts_filtered/adata.h5ad",
		inspect = "kb/{sample}/inspect.json"
	params:
		R1_R2 = lambda wildcards: samples_dict[wildcards.sample],
		kb_output = "kb/{sample}"
	threads:
		workflow.cores / 4
	benchmark:
		"benchmark/kb_count_{sample}.txt"
	log:
		"log/kb_count_{sample}.log"
	shell:
		"""
		kb count \
		-i {input.kb_index} \
		-g {input.kb_t2g} \
		-c1 {input.kb_cdna_t2c} \
		-c2 {input.kb_intron_t2c} \
		-x {config[technology]} \
		-o {params.kb_output} \
		-t {threads} \
		--workflow nac \
		--h5ad \
		--gene-names \
		--sum total \
		--filter bustools \
		--overwrite \
		--verbose \
		{params.R1_R2} >& {log}
```

**Command output (with `--verbose` flag)**
```
[2024-08-09 19:37:43,962]   DEBUG [main] Printing verbose output
[2024-08-09 19:37:46,179]   DEBUG [main] kallisto binary located at /usr/local/lib/python3.8/site-packages/kb_python/bins/linux/kallisto/kallisto
[2024-08-09 19:37:46,179]   DEBUG [main] bustools binary located at /usr/local/lib/python3.8/site-packages/kb_python/bins/linux/bustools/bustools
[2024-08-09 19:37:46,179]   DEBUG [main] Creating `kb/1122/tmp` directory
[2024-08-09 19:37:46,182]   DEBUG [main] Namespace(N=None, aa=False, batch_barcodes=False, bootstraps=None, bustools='/usr/local/lib/python3.8/site-packages/kb_python/bins/linux/bustools/bustools', c1='kb_index/cdna_t2c.txt', c2='kb_index/intron_t2c.txt', cellranger=False, chromosomes=None, command='count', dry_run=False, em=False, fastqs=['/rhome/naotok/bigdata/ADKnowledgePortal/MC_snRNA/fastq/1122_S1_R1_001.fastq.gz', '/rhome/naotok/bigdata/ADKnowledgePortal/MC_snRNA/fastq/1122_S1_R2_001.fastq.gz'], filter='bustools', filter_threshold=None, fragment_l=None, fragment_s=None, g='kb_index/t2g.txt', gene_names=True, genomebam=False, gtf=None, h5ad=True, i='kb_index/index.idx', inleaved=False, kallisto='/usr/local/lib/python3.8/site-packages/kb_python/bins/linux/kallisto/kallisto', keep_tmp=False, list=False, loom=False, loom_names='barcode,target_name', m='4G', matrix_to_directories=False, matrix_to_files=False, mm=False, no_fragment=False, no_inspect=False, no_validate=False, num=False, o='kb/1122', overwrite=True, parity=None, r=None, report=False, strand=None, sum='total', t=8, tcc=False, tmp=None, verbose=True, w=None, workflow='nac', x='10xv3')
[2024-08-09 19:37:48,862]    INFO [count_nac] Using index kb_index/index.idx to generate BUS file to kb/1122 from
[2024-08-09 19:37:48,862]    INFO [count_nac]         /rhome/naotok/bigdata/ADKnowledgePortal/MC_snRNA/fastq/1122_S1_R1_001.fastq.gz
[2024-08-09 19:37:48,862]    INFO [count_nac]         /rhome/naotok/bigdata/ADKnowledgePortal/MC_snRNA/fastq/1122_S1_R2_001.fastq.gz
[2024-08-09 19:37:48,862]   DEBUG [count_nac] kallisto bus -i kb_index/index.idx -o kb/1122 -x 10xv3 -t 8 /rhome/naotok/bigdata/ADKnowledgePortal/MC_snRNA/fastq/1122_S1_R1_001.fastq.gz /rhome/naotok/bigdata/ADKnowledgePortal/MC_snRNA/fastq/1122_S1_R2_001.fastq.gz
[2024-08-09 19:37:48,967]   DEBUG [count_nac]
[2024-08-09 19:37:48,967]   DEBUG [count_nac] [bus] Note: Strand option was not specified; setting it to --fr-stranded for specified technology
[2024-08-09 19:38:41,789]   DEBUG [count_nac] [index] k-mer length: 31
[2024-08-09 19:39:20,366]   DEBUG [count_nac] [index] number of targets: 267,211
[2024-08-09 19:39:20,466]   DEBUG [count_nac] [index] number of k-mers: 1,563,276,810
[2024-08-09 19:39:20,466]   DEBUG [count_nac] [index] number of D-list k-mers: 10,106,781
[2024-08-09 19:39:20,466]   DEBUG [count_nac] [quant] will process sample 1: /rhome/naotok/bigdata/ADKnowledgePortal/MC_snRNA/fastq/1122_S1_R1_001.fastq.gz
[2024-08-09 19:39:20,466]   DEBUG [count_nac] /rhome/naotok/bigdata/ADKnowledgePortal/MC_snRNA/fastq/1122_S1_R2_001.fastq.gz
[2024-08-09 19:39:28,081]   DEBUG [count_nac] [quant] finding pseudoalignments for the reads ...
[2024-08-09 19:39:33,894]   DEBUG [count_nac] [progress] 1M reads processed (76.1% mapped)
[2024-08-09 19:39:39,907]   DEBUG [count_nac] [progress] 2M reads processed (76.1% mapped)
[2024-08-09 19:39:46,519]   DEBUG [count_nac] [progress] 3M reads processed (76.1% mapped)
[2024-08-09 19:39:52,230]   DEBUG [count_nac] [progress] 4M reads processed (76.1% mapped)
[2024-08-09 19:39:58,342]   DEBUG [count_nac] [progress] 5M reads processed (76.1% mapped)
[2024-08-09 19:40:04,153]   DEBUG [count_nac] [progress] 6M reads processed (76.0% mapped)
[2024-08-09 19:40:09,864]   DEBUG [count_nac] [progress] 7M reads processed (76.0% mapped)
[2024-08-09 19:40:16,677]   DEBUG [count_nac] [progress] 8M reads processed (75.9% mapped)
[2024-08-09 19:40:22,288]   DEBUG [count_nac] [progress] 9M reads processed (75.9% mapped)
[2024-08-09 19:40:28,502]   DEBUG [count_nac] [progress] 10M reads processed (75.9% mapped)
[2024-08-09 19:40:34,015]   DEBUG [count_nac] [progress] 11M reads processed (75.9% mapped)
[2024-08-09 19:40:40,026]   DEBUG [count_nac] [progress] 12M reads processed (75.9% mapped)
[2024-08-09 19:40:46,338]   DEBUG [count_nac] [progress] 13M reads processed (75.9% mapped)
[2024-08-09 19:40:51,849]   DEBUG [count_nac] [progress] 14M reads processed (75.9% mapped)
[2024-08-09 19:40:57,860]   DEBUG [count_nac] [progress] 15M reads processed (75.9% mapped)
[2024-08-09 19:41:04,273]   DEBUG [count_nac] [progress] 16M reads processed (75.9% mapped)
[2024-08-09 19:41:09,784]   DEBUG [count_nac] [progress] 17M reads processed (75.9% mapped)
[2024-08-09 19:41:16,297]   DEBUG [count_nac] [progress] 18M reads processed (75.9% mapped)
[2024-08-09 19:41:21,909]   DEBUG [count_nac] [progress] 19M reads processed (75.9% mapped)
[2024-08-09 19:41:28,223]   DEBUG [count_nac] [progress] 20M reads processed (75.9% mapped)
[2024-08-09 19:41:33,934]   DEBUG [count_nac] [progress] 21M reads processed (75.8% mapped)
[2024-08-09 19:41:39,444]   DEBUG [count_nac] [progress] 22M reads processed (75.6% mapped)
[2024-08-09 19:41:45,054]   DEBUG [count_nac] [progress] 23M reads processed (75.5% mapped)
[2024-08-09 19:41:50,464]   DEBUG [count_nac] [progress] 24M reads processed (75.3% mapped)
[2024-08-09 19:41:56,679]   DEBUG [count_nac] [progress] 25M reads processed (75.2% mapped)
[2024-08-09 19:42:02,894]   DEBUG [count_nac] [progress] 26M reads processed (75.1% mapped)
[2024-08-09 19:42:08,105]   DEBUG [count_nac] [progress] 27M reads processed (75.0% mapped)
[2024-08-09 19:42:14,420]   DEBUG [count_nac] [progress] 28M reads processed (75.0% mapped)
[2024-08-09 19:42:20,230]   DEBUG [count_nac] [progress] 29M reads processed (74.9% mapped)
[2024-08-09 19:42:25,540]   DEBUG [count_nac] [progress] 30M reads processed (74.8% mapped)
[2024-08-09 19:42:31,354]   DEBUG [count_nac] [progress] 31M reads processed (74.7% mapped)
[2024-08-09 19:42:37,065]   DEBUG [count_nac] [progress] 32M reads processed (74.7% mapped)
[2024-08-09 19:42:42,876]   DEBUG [count_nac] [progress] 33M reads processed (74.7% mapped)
[2024-08-09 19:42:48,687]   DEBUG [count_nac] [progress] 34M reads processed (74.6% mapped)
[2024-08-09 19:42:54,899]   DEBUG [count_nac] [progress] 35M reads processed (74.6% mapped)
[2024-08-09 19:43:01,013]   DEBUG [count_nac] [progress] 37M reads processed (74.5% mapped)
[2024-08-09 19:43:06,323]   DEBUG [count_nac] [progress] 38M reads processed (74.5% mapped)
[2024-08-09 19:43:12,734]   DEBUG [count_nac] [progress] 39M reads processed (74.5% mapped)
[2024-08-09 19:43:18,445]   DEBUG [count_nac] [progress] 40M reads processed (74.5% mapped)
[2024-08-09 19:43:24,156]   DEBUG [count_nac] [progress] 41M reads processed (74.5% mapped)
[2024-08-09 19:43:30,970]   DEBUG [count_nac] [progress] 42M reads processed (74.6% mapped)
[2024-08-09 19:43:36,581]   DEBUG [count_nac] [progress] 43M reads processed (74.6% mapped)
[2024-08-09 19:43:42,394]   DEBUG [count_nac] [progress] 44M reads processed (74.7% mapped)
[2024-08-09 19:43:48,104]   DEBUG [count_nac] [progress] 45M reads processed (74.7% mapped)
[2024-08-09 19:43:54,318]   DEBUG [count_nac] [progress] 46M reads processed (74.8% mapped)
[2024-08-09 19:44:00,231]   DEBUG [count_nac] [progress] 47M reads processed (74.8% mapped)
[2024-08-09 19:44:06,342]   DEBUG [count_nac] [progress] 48M reads processed (74.9% mapped)
[2024-08-09 19:44:12,555]   DEBUG [count_nac] [progress] 49M reads processed (74.9% mapped)
[2024-08-09 19:44:18,167]   DEBUG [count_nac] [progress] 50M reads processed (75.0% mapped)
[2024-08-09 19:44:23,877]   DEBUG [count_nac] [progress] 51M reads processed (75.0% mapped)
[2024-08-09 19:44:30,490]   DEBUG [count_nac] [progress] 52M reads processed (75.0% mapped)
[2024-08-09 19:44:36,202]   DEBUG [count_nac] [progress] 53M reads processed (75.0% mapped)
[2024-08-09 19:44:41,717]   DEBUG [count_nac] [progress] 54M reads processed (75.1% mapped)
[2024-08-09 19:44:47,454]   DEBUG [count_nac] [progress] 55M reads processed (75.1% mapped)
[2024-08-09 19:44:53,189]   DEBUG [count_nac] [progress] 56M reads processed (75.2% mapped)
[2024-08-09 19:44:59,222]   DEBUG [count_nac] [progress] 57M reads processed (75.2% mapped)
[2024-08-09 19:45:04,754]   DEBUG [count_nac] [progress] 58M reads processed (75.2% mapped)
[2024-08-09 19:45:10,669]   DEBUG [count_nac] [progress] 59M reads processed (75.3% mapped)
[2024-08-09 19:45:16,389]   DEBUG [count_nac] [progress] 60M reads processed (75.3% mapped)
[2024-08-09 19:45:22,422]   DEBUG [count_nac] [progress] 61M reads processed (75.3% mapped)
[2024-08-09 19:45:28,344]   DEBUG [count_nac] [progress] 62M reads processed (75.3% mapped)
[2024-08-09 19:45:33,875]   DEBUG [count_nac] [progress] 63M reads processed (75.4% mapped)
[2024-08-09 19:45:39,509]   DEBUG [count_nac] [progress] 64M reads processed (75.4% mapped)
[2024-08-09 19:45:45,547]   DEBUG [count_nac] [progress] 65M reads processed (75.4% mapped)
[2024-08-09 19:45:51,182]   DEBUG [count_nac] [progress] 66M reads processed (75.4% mapped)
[2024-08-09 19:45:56,610]   DEBUG [count_nac] [progress] 67M reads processed (75.4% mapped)
[2024-08-09 19:46:02,341]   DEBUG [count_nac] [progress] 68M reads processed (75.4% mapped)
[2024-08-09 19:46:07,756]   DEBUG [count_nac] [progress] 69M reads processed (75.4% mapped)
[2024-08-09 19:46:13,168]   DEBUG [count_nac] [progress] 70M reads processed (75.3% mapped)
[2024-08-09 19:46:18,478]   DEBUG [count_nac] [progress] 71M reads processed (75.3% mapped)
[2024-08-09 19:46:24,005]   DEBUG [count_nac] [progress] 72M reads processed (75.3% mapped)
[2024-08-09 19:46:29,437]   DEBUG [count_nac] [progress] 74M reads processed (75.2% mapped)
[2024-08-09 19:46:35,876]   DEBUG [count_nac] [progress] 75M reads processed (75.2% mapped)
[2024-08-09 19:46:41,110]   DEBUG [count_nac] [progress] 76M reads processed (75.2% mapped)
[2024-08-09 19:46:46,634]   DEBUG [count_nac] [progress] 77M reads processed (75.2% mapped)
[2024-08-09 19:46:52,247]   DEBUG [count_nac] [progress] 78M reads processed (75.2% mapped)
[2024-08-09 19:46:57,675]   DEBUG [count_nac] [progress] 79M reads processed (75.1% mapped)
[2024-08-09 19:47:03,102]   DEBUG [count_nac] [progress] 80M reads processed (75.1% mapped)
[2024-08-09 19:47:09,031]   DEBUG [count_nac] [progress] 81M reads processed (75.1% mapped)
[2024-08-09 19:47:14,954]   DEBUG [count_nac] [progress] 82M reads processed (75.1% mapped)
[2024-08-09 19:47:20,568]   DEBUG [count_nac] [progress] 83M reads processed (75.1% mapped)
[2024-08-09 19:47:26,488]   DEBUG [count_nac] [progress] 84M reads processed (75.1% mapped)
[2024-08-09 19:47:32,300]   DEBUG [count_nac] [progress] 85M reads processed (75.1% mapped)
[2024-08-09 19:47:38,016]   DEBUG [count_nac] [progress] 86M reads processed (75.1% mapped)
[2024-08-09 19:47:43,436]   DEBUG [count_nac] [progress] 87M reads processed (75.1% mapped)
[2024-08-09 19:47:48,971]   DEBUG [count_nac] [progress] 88M reads processed (75.1% mapped)
[2024-08-09 19:47:54,402]   DEBUG [count_nac] [progress] 89M reads processed (75.1% mapped)
[2024-08-09 19:47:59,733]   DEBUG [count_nac] [progress] 90M reads processed (75.1% mapped)
[2024-08-09 19:48:05,468]   DEBUG [count_nac] [progress] 91M reads processed (75.1% mapped)
[2024-08-09 19:48:11,102]   DEBUG [count_nac] [progress] 92M reads processed (75.1% mapped)
[2024-08-09 19:48:17,120]   DEBUG [count_nac] [progress] 93M reads processed (75.1% mapped)
[2024-08-09 19:48:23,439]   DEBUG [count_nac] [progress] 94M reads processed (75.1% mapped)
[2024-08-09 19:48:29,555]   DEBUG [count_nac] [progress] 95M reads processed (75.1% mapped)
[2024-08-09 19:48:35,078]   DEBUG [count_nac] [progress] 96M reads processed (75.1% mapped)
[2024-08-09 19:48:40,603]   DEBUG [count_nac] [progress] 97M reads processed (75.2% mapped)
[2024-08-09 19:48:46,125]   DEBUG [count_nac] [progress] 98M reads processed (75.2% mapped)
[2024-08-09 19:48:51,753]   DEBUG [count_nac] [progress] 99M reads processed (75.2% mapped)
[2024-08-09 19:48:57,377]   DEBUG [count_nac] [progress] 100M reads processed (75.2% mapped)
[2024-08-09 19:49:02,801]   DEBUG [count_nac] [progress] 101M reads processed (75.2% mapped)
[2024-08-09 19:49:08,225]   DEBUG [count_nac] [progress] 102M reads processed (75.2% mapped)
[2024-08-09 19:49:14,853]   DEBUG [count_nac] [progress] 103M reads processed (75.2% mapped)
[2024-08-09 19:49:20,279]   DEBUG [count_nac] [progress] 104M reads processed (75.3% mapped)
[2024-08-09 19:49:25,804]   DEBUG [count_nac] [progress] 105M reads processed (75.3% mapped)
[2024-08-09 19:49:31,419]   DEBUG [count_nac] [progress] 106M reads processed (75.3% mapped)
[2024-08-09 19:49:36,829]   DEBUG [count_nac] [progress] 107M reads processed (75.3% mapped)
[2024-08-09 19:49:43,641]   DEBUG [count_nac] [progress] 108M reads processed (75.3% mapped)
[2024-08-09 19:49:49,252]   DEBUG [count_nac] [progress] 109M reads processed (75.3% mapped)
[2024-08-09 19:49:54,762]   DEBUG [count_nac] [progress] 111M reads processed (75.4% mapped)
[2024-08-09 19:50:01,375]   DEBUG [count_nac] [progress] 112M reads processed (75.4% mapped)
[2024-08-09 19:50:06,585]   DEBUG [count_nac] [progress] 113M reads processed (75.4% mapped)
[2024-08-09 19:50:11,995]   DEBUG [count_nac] [progress] 114M reads processed (75.4% mapped)
[2024-08-09 19:50:17,204]   DEBUG [count_nac] [progress] 115M reads processed (75.3% mapped)
[2024-08-09 19:50:23,517]   DEBUG [count_nac] [progress] 116M reads processed (75.3% mapped)
[2024-08-09 19:50:29,229]   DEBUG [count_nac] [progress] 117M reads processed (75.3% mapped)
[2024-08-09 19:50:34,639]   DEBUG [count_nac] [progress] 118M reads processed (75.3% mapped)
[2024-08-09 19:50:41,051]   DEBUG [count_nac] [progress] 119M reads processed (75.3% mapped)
[2024-08-09 19:50:46,160]   DEBUG [count_nac] [progress] 120M reads processed (75.2% mapped)
[2024-08-09 19:50:51,771]   DEBUG [count_nac] [progress] 121M reads processed (75.2% mapped)
[2024-08-09 19:50:57,982]   DEBUG [count_nac] [progress] 122M reads processed (75.2% mapped)
[2024-08-09 19:51:03,293]   DEBUG [count_nac] [progress] 123M reads processed (75.2% mapped)
[2024-08-09 19:51:08,904]   DEBUG [count_nac] [progress] 124M reads processed (75.2% mapped)
[2024-08-09 19:51:14,515]   DEBUG [count_nac] [progress] 125M reads processed (75.2% mapped)
[2024-08-09 19:51:20,626]   DEBUG [count_nac] [progress] 126M reads processed (75.2% mapped)
[2024-08-09 19:51:26,037]   DEBUG [count_nac] [progress] 127M reads processed (75.2% mapped)
[2024-08-09 19:51:32,148]   DEBUG [count_nac] [progress] 128M reads processed (75.2% mapped)
[2024-08-09 19:51:37,959]   DEBUG [count_nac] [progress] 129M reads processed (75.2% mapped)
[2024-08-09 19:51:43,369]   DEBUG [count_nac] [progress] 130M reads processed (75.2% mapped)
[2024-08-09 19:51:49,380]   DEBUG [count_nac] [progress] 131M reads processed (75.2% mapped)
[2024-08-09 19:51:55,291]   DEBUG [count_nac] [progress] 132M reads processed (75.2% mapped)
[2024-08-09 19:52:01,104]   DEBUG [count_nac] [progress] 133M reads processed (75.2% mapped)
[2024-08-09 19:52:06,617]   DEBUG [count_nac] [progress] 134M reads processed (75.2% mapped)
[2024-08-09 19:52:13,029]   DEBUG [count_nac] [progress] 135M reads processed (75.2% mapped)
[2024-08-09 19:52:18,940]   DEBUG [count_nac] [progress] 136M reads processed (75.2% mapped)
[2024-08-09 19:52:24,551]   DEBUG [count_nac] [progress] 137M reads processed (75.2% mapped)
[2024-08-09 19:52:30,764]   DEBUG [count_nac] [progress] 138M reads processed (75.2% mapped)
[2024-08-09 19:52:36,374]   DEBUG [count_nac] [progress] 139M reads processed (75.2% mapped)
[2024-08-09 19:52:42,285]   DEBUG [count_nac] [progress] 140M reads processed (75.3% mapped)
[2024-08-09 19:52:48,497]   DEBUG [count_nac] [progress] 141M reads processed (75.3% mapped)
[2024-08-09 19:52:54,308]   DEBUG [count_nac] [progress] 142M reads processed (75.3% mapped)
[2024-08-09 19:53:00,221]   DEBUG [count_nac] [progress] 143M reads processed (75.3% mapped)
[2024-08-09 19:53:05,937]   DEBUG [count_nac] [progress] 144M reads processed (75.3% mapped)
[2024-08-09 19:53:12,148]   DEBUG [count_nac] [progress] 145M reads processed (75.3% mapped)
[2024-08-09 19:53:18,159]   DEBUG [count_nac] [progress] 147M reads processed (75.3% mapped)
[2024-08-09 19:53:24,170]   DEBUG [count_nac] [progress] 148M reads processed (75.3% mapped)
[2024-08-09 19:53:29,883]   DEBUG [count_nac] [progress] 149M reads processed (75.3% mapped)
[2024-08-09 19:53:35,795]   DEBUG [count_nac] [progress] 150M reads processed (75.3% mapped)
[2024-08-09 19:53:41,106]   DEBUG [count_nac] [progress] 151M reads processed (75.3% mapped)
[2024-08-09 19:53:46,618]   DEBUG [count_nac] [progress] 152M reads processed (75.3% mapped)
[2024-08-09 19:53:52,531]   DEBUG [count_nac] [progress] 153M reads processed (75.3% mapped)
[2024-08-09 19:53:58,243]   DEBUG [count_nac] [progress] 154M reads processed (75.2% mapped)
[2024-08-09 19:54:03,754]   DEBUG [count_nac] [progress] 155M reads processed (75.2% mapped)
[2024-08-09 19:54:09,364]   DEBUG [count_nac] [progress] 156M reads processed (75.2% mapped)
[2024-08-09 19:54:15,290]   DEBUG [count_nac] [progress] 157M reads processed (75.2% mapped)
[2024-08-09 19:54:21,012]   DEBUG [count_nac] [progress] 158M reads processed (75.2% mapped)
[2024-08-09 19:54:26,623]   DEBUG [count_nac] [progress] 159M reads processed (75.2% mapped)
[2024-08-09 19:54:32,534]   DEBUG [count_nac] [progress] 160M reads processed (75.1% mapped)
[2024-08-09 19:54:37,947]   DEBUG [count_nac] [progress] 161M reads processed (75.1% mapped)
[2024-08-09 19:54:43,958]   DEBUG [count_nac] [progress] 162M reads processed (75.1% mapped)
[2024-08-09 19:54:49,869]   DEBUG [count_nac] [progress] 163M reads processed (75.1% mapped)
[2024-08-09 19:54:55,279]   DEBUG [count_nac] [progress] 164M reads processed (75.1% mapped)
[2024-08-09 19:55:01,695]   DEBUG [count_nac] [progress] 165M reads processed (75.1% mapped)
[2024-08-09 19:55:06,002]   DEBUG [count_nac] [progress] 166M reads processed (75.1% mapped)              done
[2024-08-09 19:55:06,002]   DEBUG [count_nac] [quant] processed 166,938,066 reads, 125,327,006 reads pseudoaligned
[2024-08-09 19:55:07,805]   DEBUG [count_nac]
[2024-08-09 19:55:17,522]    INFO [count_nac] Sorting BUS file kb/1122/output.bus to kb/1122/tmp/output.s.bus
[2024-08-09 19:55:17,522]   DEBUG [count_nac] bustools sort -o kb/1122/tmp/output.s.bus -T kb/1122/tmp -t 8 -m 4G kb/1122/output.bus
[2024-08-09 19:55:20,532]   DEBUG [count_nac] partition time: 1.79815s
[2024-08-09 19:55:23,458]   DEBUG [count_nac] all fits in buffer
[2024-08-09 19:55:25,360]   DEBUG [count_nac] Read in 125327006 BUS records
[2024-08-09 19:55:25,360]   DEBUG [count_nac] reading time 0.366752s
[2024-08-09 19:55:25,361]   DEBUG [count_nac] sorting time 19.3618s
[2024-08-09 19:55:25,361]   DEBUG [count_nac] writing time 0.902203s
[2024-08-09 19:55:25,361]    INFO [count_nac] On-list not provided
[2024-08-09 19:55:25,361]    INFO [count_nac] Copying pre-packaged 10XV3 on-list to kb/1122
[2024-08-09 19:55:25,971]    INFO [count_nac] Inspecting BUS file kb/1122/tmp/output.s.bus
[2024-08-09 19:55:25,971]   DEBUG [count_nac] bustools inspect -o kb/1122/inspect.json -w kb/1122/10x_version3_whitelist.txt kb/1122/tmp/output.s.bus
[2024-08-09 19:55:38,995]    INFO [count_nac] Correcting BUS records in kb/1122/tmp/output.s.bus to kb/1122/tmp/output.s.c.bus with on-list kb/1122/10x_version3_whitelist.txt
[2024-08-09 19:55:38,996]   DEBUG [count_nac] bustools correct -o kb/1122/tmp/output.s.c.bus -w kb/1122/10x_version3_whitelist.txt kb/1122/tmp/output.s.bus
[2024-08-09 19:55:43,609]   DEBUG [count_nac] Found 6794880 barcodes in the on-list
[2024-08-09 19:55:49,917]   DEBUG [count_nac] Processed 73699215 BUS records
[2024-08-09 19:55:49,918]   DEBUG [count_nac] In on-list = 71283991
[2024-08-09 19:55:49,918]   DEBUG [count_nac] Corrected    = 320120
[2024-08-09 19:55:49,918]   DEBUG [count_nac] Uncorrected  = 2095104
[2024-08-09 19:55:51,320]    INFO [count_nac] Sorting BUS file kb/1122/tmp/output.s.c.bus to kb/1122/output.unfiltered.bus
[2024-08-09 19:55:51,320]   DEBUG [count_nac] bustools sort -o kb/1122/output.unfiltered.bus -T kb/1122/tmp -t 8 -m 4G kb/1122/tmp/output.s.c.bus
[2024-08-09 19:55:54,032]   DEBUG [count_nac] partition time: 0.307215s
[2024-08-09 19:55:55,647]   DEBUG [count_nac] all fits in buffer
[2024-08-09 19:55:57,349]   DEBUG [count_nac] Read in 71604111 BUS records
[2024-08-09 19:55:57,350]   DEBUG [count_nac] reading time 0.222271s
[2024-08-09 19:55:57,350]   DEBUG [count_nac] sorting time 9.37299s
[2024-08-09 19:55:57,350]   DEBUG [count_nac] writing time 0.68096s
[2024-08-09 19:55:57,350]    INFO [count_nac] Generating count matrix kb/1122/counts_unfiltered/cells_x_genes from BUS file kb/1122/output.unfiltered.bus
[2024-08-09 19:55:57,350]   DEBUG [count_nac] bustools count -o kb/1122/counts_unfiltered/cells_x_genes -g kb_index/t2g.txt -e kb/1122/matrix.ec -t kb/1122/transcripts.txt -s kb_index/intron_t2c.txt --genecounts --umi-gene kb/1122/output.unfiltered.bus
[2024-08-09 19:57:29,105]   DEBUG [count_nac] kb/1122/counts_unfiltered/cells_x_genes.mature.mtx passed validation
[2024-08-09 19:57:43,330]   DEBUG [count_nac] kb/1122/counts_unfiltered/cells_x_genes.nascent.mtx passed validation
[2024-08-09 19:57:52,270]   DEBUG [count_nac] kb/1122/counts_unfiltered/cells_x_genes.ambiguous.mtx passed validation
[2024-08-09 19:57:52,270]    INFO [count_nac] Writing gene names to file kb/1122/counts_unfiltered/cells_x_genes.genes.names.txt
[2024-08-09 19:57:52,514] WARNING [count_nac] 13914 gene IDs do not have corresponding valid gene names. These genes will use their gene IDs instead.
[2024-08-09 19:57:52,567]    INFO [count_nac] Summing matrices into kb/1122/counts_unfiltered/cells_x_genes.cell.mtx
[2024-08-09 19:58:02,525]    INFO [count_nac] Summing matrices into kb/1122/counts_unfiltered/cells_x_genes.cell.mtx
[2024-08-09 19:58:16,038]    INFO [count_nac] Summing matrices into kb/1122/counts_unfiltered/cells_x_genes.nucleus.mtx
[2024-08-09 19:58:40,891]    INFO [count_nac] Summing matrices into kb/1122/counts_unfiltered/cells_x_genes.nucleus.mtx
[2024-08-09 19:59:13,397]    INFO [count_nac] Summing matrices into kb/1122/counts_unfiltered/cells_x_genes.total.mtx
[2024-08-09 19:59:34,874]    INFO [count_nac] Summing matrices into kb/1122/counts_unfiltered/cells_x_genes.total.mtx
[2024-08-09 20:00:03,929]    INFO [count_nac] Reading matrix kb/1122/counts_unfiltered/cells_x_genes.mature.mtx
[2024-08-09 20:00:05,506]    INFO [count_nac] Reading matrix kb/1122/counts_unfiltered/cells_x_genes.nascent.mtx
[2024-08-09 20:00:20,763]    INFO [count_nac] Reading matrix kb/1122/counts_unfiltered/cells_x_genes.ambiguous.mtx
[2024-08-09 20:00:30,804]    INFO [count_nac] Combining matrices
[2024-08-09 20:00:31,330]    INFO [count_nac] Writing matrices to h5ad kb/1122/counts_unfiltered/adata.h5ad
[2024-08-09 20:00:31,692]    INFO [count_nac] Filtering with bustools
[2024-08-09 20:00:31,692]    INFO [count_nac] Generating on-list kb/1122/filter_barcodes.txt from BUS file kb/1122/output.unfiltered.bus
[2024-08-09 20:00:31,692]   DEBUG [count_nac] bustools allowlist -o kb/1122/filter_barcodes.txt kb/1122/output.unfiltered.bus
[2024-08-09 20:00:33,409]   DEBUG [count_nac] Read in 71540002 BUS records, wrote 5138 barcodes to on-list with threshold 2108
[2024-08-09 20:00:33,410]    INFO [count_nac] Correcting BUS records in kb/1122/output.unfiltered.bus to kb/1122/tmp/output.unfiltered.c.bus with on-list kb/1122/filter_barcodes.txt
[2024-08-09 20:00:33,410]   DEBUG [count_nac] bustools correct -o kb/1122/tmp/output.unfiltered.c.bus -w kb/1122/filter_barcodes.txt kb/1122/output.unfiltered.bus
[2024-08-09 20:00:33,521]   DEBUG [count_nac] Found 5137 barcodes in the on-list
[2024-08-09 20:00:38,327]   DEBUG [count_nac] Processed 71540002 BUS records
[2024-08-09 20:00:38,327]   DEBUG [count_nac] In on-list = 64193659
[2024-08-09 20:00:38,328]   DEBUG [count_nac] Corrected    = 20751
[2024-08-09 20:00:38,328]   DEBUG [count_nac] Uncorrected  = 7325592
[2024-08-09 20:00:38,328]    INFO [count_nac] Sorting BUS file kb/1122/tmp/output.unfiltered.c.bus to kb/1122/output.filtered.bus
[2024-08-09 20:00:38,328]   DEBUG [count_nac] bustools sort -o kb/1122/output.filtered.bus -T kb/1122/tmp -t 8 -m 4G kb/1122/tmp/output.unfiltered.c.bus
[2024-08-09 20:00:42,549]   DEBUG [count_nac] partition time: 0.281753s
[2024-08-09 20:00:43,890]   DEBUG [count_nac] all fits in buffer
[2024-08-09 20:00:45,593]   DEBUG [count_nac] Read in 64214410 BUS records
[2024-08-09 20:00:45,593]   DEBUG [count_nac] reading time 0.195412s
[2024-08-09 20:00:45,593]   DEBUG [count_nac] sorting time 6.38335s
[2024-08-09 20:00:45,593]   DEBUG [count_nac] writing time 0.607916s
[2024-08-09 20:00:45,593]    INFO [count_nac] Generating count matrix kb/1122/counts_filtered/cells_x_genes from BUS file kb/1122/output.filtered.bus
[2024-08-09 20:00:45,593]   DEBUG [count_nac] bustools count -o kb/1122/counts_filtered/cells_x_genes -g kb_index/t2g.txt -e kb/1122/matrix.ec -t kb/1122/transcripts.txt -s kb_index/intron_t2c.txt --genecounts --umi-gene kb/1122/output.filtered.bus
[2024-08-09 20:02:13,708]   DEBUG [count_nac] kb/1122/counts_filtered/cells_x_genes.mature.mtx passed validation
[2024-08-09 20:02:25,628]   DEBUG [count_nac] kb/1122/counts_filtered/cells_x_genes.nascent.mtx passed validation
[2024-08-09 20:02:32,554]   DEBUG [count_nac] kb/1122/counts_filtered/cells_x_genes.ambiguous.mtx passed validation
[2024-08-09 20:02:32,554]    INFO [count_nac] Summing matrices into kb/1122/counts_filtered/cells_x_genes.cell.mtx
[2024-08-09 20:02:39,915]    INFO [count_nac] Summing matrices into kb/1122/counts_filtered/cells_x_genes.cell.mtx
[2024-08-09 20:02:50,039]    INFO [count_nac] Summing matrices into kb/1122/counts_filtered/cells_x_genes.nucleus.mtx
[2024-08-09 20:03:08,810]    INFO [count_nac] Summing matrices into kb/1122/counts_filtered/cells_x_genes.nucleus.mtx
[2024-08-09 20:03:33,332]    INFO [count_nac] Summing matrices into kb/1122/counts_filtered/cells_x_genes.total.mtx
[2024-08-09 20:03:49,526]    INFO [count_nac] Summing matrices into kb/1122/counts_filtered/cells_x_genes.total.mtx
[2024-08-09 20:04:11,521]    INFO [count_nac] Reading matrix kb/1122/counts_filtered/cells_x_genes.mature.mtx
[2024-08-09 20:04:12,716]    INFO [count_nac] Reading matrix kb/1122/counts_filtered/cells_x_genes.nascent.mtx
[2024-08-09 20:04:24,785]    INFO [count_nac] Reading matrix kb/1122/counts_filtered/cells_x_genes.ambiguous.mtx
[2024-08-09 20:04:32,155]    INFO [count_nac] Combining matrices
[2024-08-09 20:04:32,322]    INFO [count_nac] Writing matrices to h5ad kb/1122/counts_filtered/adata.h5ad
[2024-08-09 20:04:32,458]   DEBUG [main] Removing `kb/1122/tmp` directory
```



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index d2901fd..869ce5c 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -22,7 +22,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python: [3.7, 3.8, 3.9 ]
+        python: [3.8, 3.9 ]
         os: [ubuntu-20.04]
     name: Test on Python ${{ matrix.python }}
     steps:
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 3a4e66e..db31af0 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -14,7 +14,7 @@ jobs:
       - name: Setup python
         uses: actions/setup-python@v1
         with:
-          python-version: '3.7'
+          python-version: '3.8'
           architecture: x64
       - name: Install dependencies
         run: pip install -r dev-requirements.txt
diff --git a/.gitignore b/.gitignore
index e3d1e4c..9f4d7d9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -109,3 +109,9 @@ venv.bak/
 .DS_Store
 .vscode/
 .Rhistory
+
+# PyCharm
+/.idea/
+
+# Temp files
+/scratch/
\ No newline at end of file
diff --git a/Makefile b/Makefile
index dc7cb6b..6549c37 100755
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,8 @@
 
 test:
 	rm -f .coverage
-	nosetests --verbose --with-coverage --cover-package kb_python tests/* tests/dry/*
+	pytest --verbose --cov=kb_python tests/* tests/dry/* && coverage report && coverage xml
+#	nosetests --verbose --with-coverage --cover-package kb_python tests/* tests/dry/*
 
 check:
 	flake8 kb_python && echo OK
diff --git a/README.md b/README.md
index cda1eb0..a661fbb 100755
--- a/README.md
+++ b/README.md
@@ -1,5 +1,5 @@
 # kb-python
-![github version](https://img.shields.io/badge/Version-0.28.0-informational)
+![github version](https://img.shields.io/badge/Version-0.29.0-informational)
 [![pypi version](https://img.shields.io/pypi/v/kb-python)](https://pypi.org/project/kb-python/0.28.0/)
 ![python versions](https://img.shields.io/pypi/pyversions/kb_python)
 ![status](https://github.com/pachterlab/kb_python/workflows/CI/badge.svg)
@@ -10,7 +10,7 @@
 
 `kb-python` is a python package for processing single-cell RNA-sequencing. It wraps the [`kallisto` | `bustools`](https://www.kallistobus.tools) single-cell RNA-seq command line tools in order to unify multiple processing workflows. 
 
-`kb-python` was developed by [Kyung Hoi (Joseph) Min](https://twitter.com/lioscro) and [A. Sina Booeshaghi](https://twitter.com/sinabooeshaghi) while in [Lior Pachter](https://twitter.com/lpachter)'s lab at Caltech. If you use `kb-python` in a publication please [cite*](#cite):
+`kb-python` was first developed by [Kyung Hoi (Joseph) Min](https://twitter.com/lioscro) and [A. Sina Booeshaghi](https://twitter.com/sinabooeshaghi) while in [Lior Pachter](https://twitter.com/lpachter)'s lab at Caltech. If you use `kb-python` in a publication please [cite*](#cite):
 ```
 Melsted, P., Booeshaghi, A.S., et al. 
 Modular, efficient and constant-memory single-cell RNA-seq preprocessing. 
@@ -34,7 +34,7 @@ There are no prerequisite packages to install. The `kallisto` and `bustools` bin
 
 ## Usage
 
-`kb`  consists of four subcommands
+`kb`  consists of five subcommands
 ```bash
 $ kb
 usage: kb [-h] [--list] <CMD> ...
@@ -44,6 +44,7 @@ positional arguments:
     compile   Compile `kallisto` and `bustools` binaries from source
     ref       Build a kallisto index and transcript-to-gene mapping
     count     Generate count matrices from a set of single-cell FASTQ files
+    extract   Extract reads that were pseudoaligned to specific genes/transcripts (or extract all reads that were / were not pseudoaligned)
 ```
 
 ### `kb ref`: generate a pseudoalignment index
diff --git a/dev-requirements.txt b/dev-requirements.txt
index 91b4acd..e49181d 100644
--- a/dev-requirements.txt
+++ b/dev-requirements.txt
@@ -1,7 +1,8 @@
 bumpversion==0.6.0
-coverage==5.1
+coverage==5.2.1
 flake8==3.8.2
-nose==1.3.7
+pytest==8.2.2
+pytest-cov==5.0.0
 pre-commit==2.4.0
 sphinx>=3.3.1
 sphinx-autoapi>=1.5.1
diff --git a/docs/conf.py b/docs/conf.py
index 6d4c0b0..54aeb59 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -24,7 +24,7 @@
 author = 'Kyung Hoi (Joseph) Min'
 
 # The full version, including alpha/beta/rc tags
-release = '0.28.2'
+release = '0.29.0'
 master_doc = 'index'
 
 # -- General configuration ---------------------------------------------------
diff --git a/docs/index.rst b/docs/index.rst
index 640b9b0..e6d5cd6 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -6,7 +6,7 @@
 Welcome to kb-python's documentation!
 =====================================
 
-This page contains **DEVELOPER** documentation for ``kb-python`` version ``0.28.2``.
+This page contains **DEVELOPER** documentation for ``kb-python`` version ``0.29.0``.
 For user documentation and tutorials, please go to `kallisto | bustools <https://www.kallistobus.tools/>`_.
 
 Development Prerequisites
@@ -18,7 +18,7 @@ necessary packages by running::
   pip install -r requirements.txt
   pip install -r dev-requirements.txt
 
-Code qualty and unit tests are strictly enforced for every pull request via
+Code quality and unit tests are strictly enforced for every pull request via
 Github actions.
 
 Code Quality
@@ -33,7 +33,7 @@ at the root of the repository.
 
 Unit-testing
 """"""""""""
-``kb-python`` uses ``nose`` to run unit tests. There is a convenient Makefile
+``kb-python`` uses ``pytest`` to run unit tests. There is a convenient Makefile
 rule in place to run all tests.::
 
   make test
diff --git a/kb_python/__init__.py b/kb_python/__init__.py
index 873054e..9093e4e 100644
--- a/kb_python/__init__.py
+++ b/kb_python/__init__.py
@@ -1,1 +1,1 @@
-__version__ = '0.28.2'
+__version__ = '0.29.0'
diff --git a/kb_python/bins/darwin/bustools/bustools b/kb_python/bins/darwin/bustools/bustools
index 1c788fc..d192ab6 100755
Binary files a/kb_python/bins/darwin/bustools/bustools and b/kb_python/bins/darwin/bustools/bustools differ
diff --git a/kb_python/bins/darwin/kallisto/kallisto b/kb_python/bins/darwin/kallisto/kallisto
index 8deba31..9d12937 100755
Binary files a/kb_python/bins/darwin/kallisto/kallisto and b/kb_python/bins/darwin/kallisto/kallisto differ
diff --git a/kb_python/bins/darwin/kallisto/kallisto_k64 b/kb_python/bins/darwin/kallisto/kallisto_k64
new file mode 100755
index 0000000..daeb012
Binary files /dev/null and b/kb_python/bins/darwin/kallisto/kallisto_k64 differ
diff --git a/kb_python/bins/darwin/kallisto/kallisto_optoff b/kb_python/bins/darwin/kallisto/kallisto_optoff
new file mode 100755
index 0000000..98d5ec3
Binary files /dev/null and b/kb_python/bins/darwin/kallisto/kallisto_optoff differ
diff --git a/kb_python/bins/darwin/kallisto/kallisto_optoff_k64 b/kb_python/bins/darwin/kallisto/kallisto_optoff_k64
new file mode 100755
index 0000000..e1eb778
Binary files /dev/null and b/kb_python/bins/darwin/kallisto/kallisto_optoff_k64 differ
diff --git a/kb_python/bins/darwin/m1/bustools/bustools b/kb_python/bins/darwin/m1/bustools/bustools
index 1c788fc..d192ab6 100755
Binary files a/kb_python/bins/darwin/m1/bustools/bustools and b/kb_python/bins/darwin/m1/bustools/bustools differ
diff --git a/kb_python/bins/darwin/m1/kallisto/kallisto b/kb_python/bins/darwin/m1/kallisto/kallisto
index 544460a..65db552 100755
Binary files a/kb_python/bins/darwin/m1/kallisto/kallisto and b/kb_python/bins/darwin/m1/kallisto/kallisto differ
diff --git a/kb_python/bins/darwin/m1/kallisto/kallisto_k64 b/kb_python/bins/darwin/m1/kallisto/kallisto_k64
new file mode 100755
index 0000000..6d3985e
Binary files /dev/null and b/kb_python/bins/darwin/m1/kallisto/kallisto_k64 differ
diff --git a/kb_python/bins/darwin/m1/kallisto/kallisto_optoff b/kb_python/bins/darwin/m1/kallisto/kallisto_optoff
new file mode 100755
index 0000000..e38af8e
Binary files /dev/null and b/kb_python/bins/darwin/m1/kallisto/kallisto_optoff differ
diff --git a/kb_python/bins/darwin/m1/kallisto/kallisto_optoff_k64 b/kb_python/bins/darwin/m1/kallisto/kallisto_optoff_k64
new file mode 100755
index 0000000..102905a
Binary files /dev/null and b/kb_python/bins/darwin/m1/kallisto/kallisto_optoff_k64 differ
diff --git a/kb_python/bins/linux/bustools/bustools b/kb_python/bins/linux/bustools/bustools
index 7f7ec41..c4b0ce3 100755
Binary files a/kb_python/bins/linux/bustools/bustools and b/kb_python/bins/linux/bustools/bustools differ
diff --git a/kb_python/bins/linux/kallisto/kallisto b/kb_python/bins/linux/kallisto/kallisto
index 2fd8e41..82a141e 100755
Binary files a/kb_python/bins/linux/kallisto/kallisto and b/kb_python/bins/linux/kallisto/kallisto differ
diff --git a/kb_python/bins/linux/kallisto/kallisto_k64 b/kb_python/bins/linux/kallisto/kallisto_k64
new file mode 100755
index 0000000..f0c85ec
Binary files /dev/null and b/kb_python/bins/linux/kallisto/kallisto_k64 differ
diff --git a/kb_python/bins/linux/kallisto/kallisto_optoff b/kb_python/bins/linux/kallisto/kallisto_optoff
new file mode 100755
index 0000000..01233e6
Binary files /dev/null and b/kb_python/bins/linux/kallisto/kallisto_optoff differ
diff --git a/kb_python/bins/linux/kallisto/kallisto_optoff_k64 b/kb_python/bins/linux/kallisto/kallisto_optoff_k64
new file mode 100755
index 0000000..05ccebf
Binary files /dev/null and b/kb_python/bins/linux/kallisto/kallisto_optoff_k64 differ
diff --git a/kb_python/bins/windows/bustools/bustools.exe b/kb_python/bins/windows/bustools/bustools.exe
index 28b7e81..95e79e0 100755
Binary files a/kb_python/bins/windows/bustools/bustools.exe and b/kb_python/bins/windows/bustools/bustools.exe differ
diff --git a/kb_python/bins/windows/kallisto/kallisto.exe b/kb_python/bins/windows/kallisto/kallisto.exe
index a6f8035..7f30441 100755
Binary files a/kb_python/bins/windows/kallisto/kallisto.exe and b/kb_python/bins/windows/kallisto/kallisto.exe differ
diff --git a/kb_python/bins/windows/kallisto/kallisto_k64.exe b/kb_python/bins/windows/kallisto/kallisto_k64.exe
new file mode 100755
index 0000000..9a92dc7
Binary files /dev/null and b/kb_python/bins/windows/kallisto/kallisto_k64.exe differ
diff --git a/kb_python/bins/windows/kallisto/kallisto_optoff.exe b/kb_python/bins/windows/kallisto/kallisto_optoff.exe
new file mode 100755
index 0000000..94bfc5a
Binary files /dev/null and b/kb_python/bins/windows/kallisto/kallisto_optoff.exe differ
diff --git a/kb_python/bins/windows/kallisto/kallisto_optoff_k64.exe b/kb_python/bins/windows/kallisto/kallisto_optoff_k64.exe
new file mode 100755
index 0000000..3dc3666
Binary files /dev/null and b/kb_python/bins/windows/kallisto/kallisto_optoff_k64.exe differ
diff --git a/kb_python/config.py b/kb_python/config.py
index 50133e4..f4b1919 100755
--- a/kb_python/config.py
+++ b/kb_python/config.py
@@ -34,7 +34,14 @@ def get_provided_kallisto_path() -> Optional[str]:
     Returns:
         Path to the binary, `None` if not found
     """
-    bin_filename = 'kallisto.exe' if PLATFORM == 'windows' else 'kallisto'
+    bin_name = 'kallisto'
+    if '_KALLISTO_OPTOFF' in globals():
+        if _KALLISTO_OPTOFF:
+            bin_name = f'{bin_name}_optoff'
+    if '_KALLISTO_KMER_64' in globals():
+        if _KALLISTO_KMER_64:
+            bin_name = f'{bin_name}_k64'
+    bin_filename = f'{bin_name}.exe' if PLATFORM == 'windows' else bin_name
     path = os.path.join(BINS_DIR, PLATFORM, CPU, 'kallisto', bin_filename)
     if not os.path.isfile(path):
         return None
@@ -54,11 +61,18 @@ def get_provided_bustools_path() -> Optional[str]:
     return path
 
 
+def set_special_kallisto_binary(k64: bool, optoff: bool):
+    global _KALLISTO_KMER_64
+    global _KALLISTO_OPTOFF
+    _KALLISTO_KMER_64 = k64
+    _KALLISTO_OPTOFF = optoff
+
+
 def get_compiled_kallisto_path(alias: str = COMPILED_DIR) -> Optional[str]:
     """Finds platform-dependent kallisto binary compiled with `compile`.
 
     Args:
-        Alias: Alias of compiled binary.
+        alias: Alias of compiled binary.
 
     Returns:
         Path to the binary, `None` if not found
@@ -74,7 +88,7 @@ def get_compiled_bustools_path(alias: str = COMPILED_DIR) -> Optional[str]:
     """Finds platform-dependent bustools binary compiled with `compile`.
 
     Args:
-        Alias: Alias of compiled binary.
+        alias: Alias of compiled binary.
 
     Returns:
         Path to the binary, `None` if not found
diff --git a/kb_python/count.py b/kb_python/count.py
index bf5d328..9c105b3 100755
--- a/kb_python/count.py
+++ b/kb_python/count.py
@@ -105,6 +105,11 @@ def kallisto_bus(
     demultiplexed: bool = False,
     batch_barcodes: bool = False,
     numreads: int = None,
+    lr: bool = False,
+    lr_thresh: float = 0.8,
+    lr_error_rate: float = None,
+    union: bool = False,
+    no_jump: bool = False,
 ) -> Dict[str, str]:
     """Runs `kallisto bus`.
 
@@ -133,6 +138,11 @@ def kallisto_bus(
         demultiplexed: Whether FASTQs are demultiplexed, defaults to `False`
         batch_barcodes: Whether sample ID should be in barcode, defaults to `False`
         numreads: Maximum number of reads to process from supplied input
+        lr: Whether to use lr-kallisto in read mapping, defaults to `False`
+        lr_thresh: Sets the --threshold for lr-kallisto, defaults to `0.8`
+        lr_error_rate: Sets the --error-rate for lr-kallisto, defaults to `None`
+        union: Use set union for pseudoalignment, defaults to `False`
+        no_jump: Disable pseudoalignment "jumping", defaults to `False`
 
     Returns:
         Dictionary containing paths to generated files
@@ -194,6 +204,16 @@ def kallisto_bus(
         command += ['--rf-stranded']
     if inleaved:
         command += ['--inleaved']
+    if lr:
+        command += ['--long']
+    if lr and lr_thresh:
+        command += ['-r', str(lr_thresh)]
+    if lr and lr_error_rate:
+        command += ['-e', str(lr_error_rate)]
+    if union:
+        command += ['--union']
+    if no_jump:
+        command += ['--no-jump']
     if batch_barcodes:
         command += ['--batch-barcodes']
     if is_batch:
@@ -224,12 +244,14 @@ def kallisto_quant_tcc(
     matrix_to_files: bool = False,
     matrix_to_directories: bool = False,
     no_fragment: bool = False,
+    lr: bool = False,
+    lr_platform: str = 'ONT',
 ) -> Dict[str, str]:
     """Runs `kallisto quant-tcc`.
 
     Args:
         mtx_path: Path to counts matrix
-        saved_index_path: Path to index.saved
+        saved_index_path: Path to index
         ecmap_path: Path to ecmap
         t2g_path: Path to T2G
         out_dir: Output directory path
@@ -241,6 +263,8 @@ def kallisto_quant_tcc(
         matrix_to_files: Whether to write quant-tcc output to files, defaults to `False`
         matrix_to_directories: Whether to write quant-tcc output to directories, defaults to `False`
         no_fragment: Whether to disable quant-tcc effective length normalization, defaults to `False`
+        lr: Whether to use lr-kallisto in quantification, defaults to `False`
+        lr_platform: Sets the --platform for lr-kallisto, defaults to `ONT`
 
     Returns:
         Dictionary containing path to output files
@@ -255,6 +279,10 @@ def kallisto_quant_tcc(
     command += ['-e', ecmap_path]
     command += ['-g', t2g_path]
     command += ['-t', threads]
+    if lr:
+        command += ['--long']
+    if lr and lr_platform:
+        command += ['-P', lr_platform]
     if flens_path and not no_fragment:
         command += ['-f', flens_path]
     if l and not no_fragment:
@@ -1178,6 +1206,14 @@ def count(
     no_fragment: bool = False,
     numreads: int = None,
     store_num: bool = False,
+    lr: bool = False,
+    lr_thresh: float = 0.8,
+    lr_error_rate: float = None,
+    lr_platform: str = 'ONT',
+    union: bool = False,
+    no_jump: bool = False,
+    quant_umis: bool = False,
+    keep_flags: bool = False,
 ) -> Dict[str, Union[str, Dict[str, str]]]:
     """Generates count matrices for single-cell RNA seq.
 
@@ -1242,6 +1278,14 @@ def count(
         no_fragment: Whether to disable quant-tcc effective length normalization, defaults to `False`
         numreads: Maximum number of reads to process from supplied input
         store_num: Whether to store read numbers in BUS file, defaults to `False`
+        lr: Whether to use lr-kallisto in read mapping, defaults to `False`
+        lr_thresh: Sets the --threshold for lr-kallisto, defaults to `0.8`
+        lr_error_rate: Sets the --error-rate for lr-kallisto, defaults to `None`
+        lr_platform: Sets the --platform for lr-kallisto, defaults to `ONT`
+        union: Use set union for pseudoalignment, defaults to `False`
+        no_jump: Disable pseudoalignment "jumping", defaults to `False`
+        quant_umis: Whether to run quant-tcc when there are UMIs, defaults to `False`
+        keep_flags: Preserve flag column when sorting BUS file, defaults to `False`
 
     Returns:
         Dictionary containing paths to generated files
@@ -1292,7 +1336,12 @@ def count(
             demultiplexed=demultiplexed,
             batch_barcodes=batch_barcodes,
             numreads=numreads,
-            n=store_num
+            n=store_num,
+            lr=lr,
+            lr_thresh=lr_thresh,
+            lr_error_rate=lr_error_rate,
+            union=union,
+            no_jump=no_jump
         )
     else:
         logger.info(
@@ -1309,7 +1358,7 @@ def count(
         temp_dir=temp_dir,
         threads=threads,
         memory=memory,
-        store_num=store_num
+        store_num=store_num and not keep_flags
     )
     correct = True
     if whitelist_path and whitelist_path.upper() == "NONE":
@@ -1404,6 +1453,9 @@ def update_results_with_suffix(current_results, new_results, suffix):
         technology.upper() in ('BULK', 'SMARTSEQ2', 'SMARTSEQ3')
     ) or ignore_umis
     quant = cm and tcc
+    if quant_umis:
+        quant = True
+        no_fragment = True
     suffix_to_inspect_filename = {'': ''}
     if (technology.upper() == 'SMARTSEQ3'):
         suffix_to_inspect_filename = {
@@ -1518,6 +1570,8 @@ def update_results_with_suffix(current_results, new_results, suffix):
                     matrix_to_files=matrix_to_files,
                     matrix_to_directories=matrix_to_directories,
                     no_fragment=no_fragment,
+                    lr=lr,
+                    lr_platform=lr_platform,
                 )
                 update_results_with_suffix(
                     unfiltered_results, quant_result, suffix
@@ -1695,6 +1749,14 @@ def count_nac(
     batch_barcodes: bool = False,
     numreads: int = None,
     store_num: bool = False,
+    lr: bool = False,
+    lr_thresh: float = 0.8,
+    lr_error_rate: float = None,
+    lr_platform: str = 'ONT',
+    union: bool = False,
+    no_jump: bool = False,
+    quant_umis: bool = False,
+    keep_flags: bool = False,
 ) -> Dict[str, Union[Dict[str, str], str]]:
     """Generates RNA velocity matrices for single-cell RNA seq.
 
@@ -1756,6 +1818,14 @@ def count_nac(
         batch_barcodes: Whether sample ID should be in barcode, defaults to `False`
         numreads: Maximum number of reads to process from supplied input
         store_num: Whether to store read numbers in BUS file, defaults to `False`
+        lr: Whether to use lr-kallisto in read mapping, defaults to `False`
+        lr_thresh: Sets the --threshold for lr-kallisto, defaults to `0.8`
+        lr_error_rate: Sets the --error-rate for lr-kallisto, defaults to `None`
+        lr_platform: Sets the --platform for lr-kallisto, defaults to `ONT`
+        union: Use set union for pseudoalignment, defaults to `False`
+        no_jump: Disable pseudoalignment "jumping", defaults to `False`
+        quant_umis: Whether to run quant-tcc when there are UMIs, defaults to `False`
+        keep_flags: Preserve flag column when sorting BUS file, defaults to `False`
 
     Returns:
         Dictionary containing path to generated index
@@ -1803,7 +1873,12 @@ def count_nac(
             demultiplexed=demultiplexed,
             batch_barcodes=batch_barcodes,
             numreads=numreads,
-            n=store_num
+            n=store_num,
+            lr=lr,
+            lr_thresh=lr_thresh,
+            lr_error_rate=lr_error_rate,
+            union=union,
+            no_jump=no_jump
         )
     else:
         logger.info(
@@ -1820,7 +1895,7 @@ def count_nac(
         temp_dir=temp_dir,
         threads=threads,
         memory=memory,
-        store_num=store_num
+        store_num=store_num and not keep_flags
     )
     correct = True
     if whitelist_path and whitelist_path.upper() == "NONE":
@@ -2073,8 +2148,8 @@ def update_results_with_suffix(current_results, new_results, suffix):
                         if batch_barcodes else None for prefix in prefixes
                     ],
                     genes_paths=[
-                        unfiltered_results[prefix][f'txnames{suffix}'] if tcc
-                        else unfiltered_results[prefix].get(f'genes{suffix}')
+                        unfiltered_results[prefix][f'ec{suffix}'] if tcc else
+                        unfiltered_results[prefix].get(f'genes{suffix}')
                         for prefix in prefixes
                     ],
                     t2g_path=t2g_path,
@@ -2975,7 +3050,7 @@ def update_results_with_suffix(current_results, new_results, suffix):
                     for prefix in prefixes
                 ],
                 genes_paths=[
-                    unfiltered_results[prefix][f'txnames{suffix}'] if tcc else
+                    unfiltered_results[prefix][f'ec{suffix}'] if tcc else
                     unfiltered_results[prefix].get(f'genes{suffix}')
                     for prefix in prefixes
                 ],
diff --git a/kb_python/extract.py b/kb_python/extract.py
new file mode 100644
index 0000000..df5ebed
--- /dev/null
+++ b/kb_python/extract.py
@@ -0,0 +1,503 @@
+import os
+import gzip
+from typing import Dict, List, Optional, Union
+
+from Bio import SeqIO
+import pandas as pd
+import numpy as np
+
+from typing_extensions import Literal
+
+from .utils import make_directory, get_bustools_binary_path, run_executable, get_bustools_version
+
+from .count import kallisto_bus, bustools_capture, bustools_sort, stream_fastqs
+
+from .logging import logger
+
+
+def bustools_text(
+    bus_path: str,
+    out_path: str,
+    flags: bool = False,
+):
+    """Runs `bustools text`.
+
+    Args:
+        bus_path: Path to BUS file to convert to text format
+        out_dir: Path to output txt file path
+        flags: Whether to include the flags columns
+    """
+    # logger.info('Converting BUS file {} to {}'.format(bus_path, out_path))
+    command = [get_bustools_binary_path(), "text"]
+    command += ["-o", out_path]
+    if flags:
+        command += ["--flags"]
+    command += [bus_path]
+    run_executable(command)
+    return {"bus": out_path}
+
+
+def bustools_fromtext(txt_path: str, out_path: str):
+    """Runs `bustools fromtext`.
+
+    Args:
+        bus_path: Path to text file to convert to BUS format
+        out_dir: Path to output BUS file
+    """
+    # logger.info('Creating BUS file {} from {}'.format(out_path, bus_path))
+    command = [get_bustools_binary_path(), "fromtext"]
+    command += ["-o", out_path]
+    command += [txt_path]
+    run_executable(command)
+    return {"bus": out_path}
+
+
+def bustools_extract(
+    sorted_bus_path: str,
+    out_path: str,
+    fastqs: Union[str, List[str]],
+    num_fastqs: int,
+) -> Dict[str, str]:
+    """
+    Extract reads from a BUS file using bustools.
+    Args:
+        sorted_bus_path: path to a BUS file sorted by flag using bustools sort --flag
+        out_path: Output directory for FASTQ files
+        fastqs: FASTQ file(s) from which to extract reads
+        num_fastqs: Number of FASTQ file(s) per run
+
+    Returns:
+        Dictionary containing path to generated BUS file
+    """
+    logger.info(f"Extracting BUS file {sorted_bus_path} to {out_path}")
+    command = [get_bustools_binary_path(), "extract"]
+
+    if not isinstance(fastqs, list):
+        fastqs = [fastqs]
+
+    command += ["-o", out_path]
+    command += ["-f", ",".join(fastqs)]
+    command += ["-N", num_fastqs]
+    command += [sorted_bus_path]
+    run_executable(command)
+    return {"bus": out_path}
+
+
+def is_gzipped(file_path):
+    """
+    Checks if a file is gzipped by reading its magic number.
+    """
+    with open(file_path, 'rb') as file:
+        return file.read(2) == b'\x1f\x8b'
+
+
+def read_headers_from_fastq(fastq_file):
+    """
+    Reads headers from a FASTQ file and returns a set of headers.
+    """
+    headers = set()
+    with gzip.open(fastq_file, "rt") as file:
+        for record in SeqIO.parse(file, "fastq"):
+            headers.add(record.id)
+    return headers
+
+
+def extract_matching_reads_by_header(
+    input_fastq, reference_fastq, output_fastq
+):
+    """
+    Extracts reads from the reference FASTQ (.gz) file that are NOT present in the input FASTQ (.gz) file
+    based on headers and writes them to the output FASTQ (.gz) file.
+    """
+
+    # Read headers from the reference FASTQ file
+    reference_headers = read_headers_from_fastq(input_fastq)
+
+    # Determine if reference_fastq is gzipped and open accordingly
+    if is_gzipped(reference_fastq):
+        infile_opener = gzip.open(reference_fastq, "rt")
+    else:
+        infile_opener = open(reference_fastq, "r")
+
+    # Ensure all intermediary folders exist for the output_fastq path
+    output_dir = os.path.dirname(output_fastq)
+    os.makedirs(output_dir, exist_ok=True)
+
+    with infile_opener as infile, gzip.open(output_fastq, "wt") as outfile:
+        # Create a SeqIO writer for the output FASTQ file
+        writer = SeqIO.write(
+            (
+                record for record in SeqIO.parse(infile, "fastq")
+                if record.id not in reference_headers
+            ),
+            outfile,
+            "fastq",
+        )
+        print(f"Number of unmapped reads written to {output_fastq}: {writer}")
+
+
+def get_mm_ecs(t2g_path, txnames, temp_dir):
+    """
+    Gets a list of equivalence classes that map to multiple genes.
+    """
+    logger.debug("Fetching equivalence classes that map to multiple genes...")
+
+    ecmap = os.path.join(temp_dir, "matrix.ec")
+
+    # Read t2g to find all transcripts associated with a gene
+    with open(t2g_path, "r") as t2g_file:
+        lines = t2g_file.readlines()
+    t2g_df = pd.DataFrame()
+    t2g_df["transcript"] = [line.split("\t")[0] for line in lines]
+    t2g_df["gene_id"] = [
+        line.split("\t")[1].replace("\n", "") for line in lines
+    ]
+
+    with open(txnames) as f:
+        txs = f.read().splitlines()
+
+    ec_df = pd.read_csv(ecmap, sep="\t", header=None)
+    # List to save multimapped ecs
+    ecs_mm = []
+    for index, row in ec_df.iterrows():
+        # Get transcript IDs that mapped to this ec
+        if isinstance(row[1], np.int64) or isinstance(row[1], int):
+            mapped_txs = np.array(txs)[row[1]]
+            mapped_txs = [mapped_txs]
+        else:
+            mapped_txs = np.array(txs)[np.array(row[1].split(",")).astype(int)]
+
+        # Check if transcript IDs belong to one or more genes
+        if (len(set(
+                t2g_df[t2g_df["transcript"].isin(mapped_txs)]["gene_id"].values)
+                ) > 1):
+            ecs_mm.append(row[0])
+
+    logger.debug(
+        f"Found the following equivalence classes that map to multiple genes: {' ,'.join(np.array(ecs_mm).astype(str))}"
+    )
+
+    return ecs_mm, ec_df
+
+
+def remove_mm_from_bus(t2g_path, txnames, temp_dir, bus_in):
+    """
+    Removes rows containing equivalence classes that map to multiple genes from the bus file.
+    """
+    logger.debug(
+        f"Removing equivalence classes that map to multiple genes from {bus_in}"
+    )
+
+    # Get equivalence classes with multimapped reads
+    ecs_mm, _ = get_mm_ecs(t2g_path, txnames, temp_dir)
+
+    if len(ecs_mm) > 0:
+        # Remove mm ecs from bus file
+        bus_txt = os.path.join(temp_dir, "output.bus.txt")
+        bus_txt_no_mm = os.path.join(temp_dir, "output_no_mm.bus.txt")
+        bus_no_mm = os.path.join(temp_dir, "output_no_mm.bus")
+
+        # Convert bus to txt file
+        bustools_text(bus_path=bus_in, out_path=bus_txt, flags=True)
+
+        # Remove mm ecs
+        bus_df = pd.read_csv(bus_txt, sep="\t", header=None)
+        new_bus_df = bus_df[~bus_df[2].isin(ecs_mm)]
+        new_bus_df.to_csv(bus_txt_no_mm, sep="\t", index=False, header=None)
+
+        # Convert back to bus format
+        bustools_fromtext(txt_path=bus_txt_no_mm, out_path=bus_no_mm)
+
+        logger.debug(
+            f"BUS file without equivalence classes that map to multiple genes saved at {bus_no_mm}"
+        )
+
+        return bus_no_mm
+
+    else:
+        logger.debug("No equivalence classes that map to multiple genes found.")
+        return None
+
+
+def remove_mm_from_mc(t2g_path, txnames, temp_dir):
+    """
+    Replaces transcript entries with -1 for equivalence classes that map to multiple genes in the matrix.ec file.
+    (These cannot simply be removed, otherwise bustools capture cannot parse the equivalence classes anymore.
+    """
+    ecmap_no_mm = os.path.join(temp_dir, "matrix_no_mm.ec")
+
+    logger.debug(
+        "Replacing transcript entries with -1 for equivalence classes "
+        f"that map to multiple genes from {os.path.join(temp_dir, 'matrix.ec')}"
+    )
+
+    # Get multimapped equivalence classes
+    ecs_mm, ec_df = get_mm_ecs(t2g_path, txnames, temp_dir)
+
+    if len(ecs_mm) > 0:
+        # Replace transcript entries for multimapped equivalence classes with -1
+        ec_df.loc[ec_df[0].isin(ecs_mm), 1] = -1
+        ec_df.to_csv(ecmap_no_mm, sep="\t", index=False, header=None)
+
+        logger.debug(
+            "matrix.ec file where transcript entries were replaced with -1 for "
+            f"equivalence classes that map to multiple genes saved at {ecmap_no_mm}"
+        )
+
+        return ecmap_no_mm
+
+    else:
+        logger.debug("No equivalence classes that map to multiple genes found.")
+        return None
+
+
+@logger.namespaced("extract")
+def extract(
+    fastq: str,
+    index_path: str,
+    targets: list[str],
+    out_dir: str,
+    target_type: Literal["gene", "transcript"],
+    extract_all: bool = False,
+    extract_all_fast: bool = False,
+    extract_all_unmapped: bool = False,
+    mm: bool = False,
+    t2g_path: Optional[str] = None,
+    temp_dir: str = "tmp",
+    threads: int = 8,
+    aa: bool = False,
+    strand: Optional[Literal["unstranded", "forward", "reverse"]] = None,
+    numreads: Optional[int] = None,
+):
+    """
+    Extracts sequencing reads that were pseudo-aligned to an index for specific genes/transcripts.
+
+    fastq: Single fastq file containing sequencing reads
+    index_path: Path to kallisto index
+    targets: Gene or transcript names for which to extract the raw reads that align to the index
+    out_dir: Path to output directory
+    target_type: 'gene' (default) or 'transcript' -> Defines whether targets are gene or transcript names
+    extract_all: Extracts reads for all genes or transcripts (as defined in target_type), defaults to `False`.
+        Might take a long time to run when the reference index contains a large number of genes.
+        Set targets = None when using extract_all
+    extract_all_fast: Extracts all pseudo-aligned reads, defaults to `False`.
+        Does not break down output by gene/transcript.
+        Set targets = None when using extract_all_fast
+    extract_all_unmapped: Extracts all unmapped reads, defaults to `False`.
+        Set targets = None when using extract_all_unmapped
+    mm: Also extract reads that multi-mapped to several genes, defaults to `False`
+    t2g_path: Path to transcript-to-gene mapping file
+        (required when mm = False, target_type = 'gene'
+        (and extract_all_fast and extract_all_unmapped = False),
+        OR extract_all = True)
+    temp_dir: Path to temporary directory, defaults to `tmp`
+    threads: Number of threads to use, defaults to `8`
+    aa: Align to index generated from a FASTA-file containing amino acid sequences, defaults to `False`
+    strand: Strandedness, defaults to `None`
+    numreads: Maximum number of reads to process from supplied input
+
+    Returns:
+    Raw reads that were pseudo-aligned to the index by kallisto for each specified gene/transcript.
+    """
+    if sum([extract_all, extract_all_fast, extract_all_unmapped]) > 1:
+        raise ValueError(
+            "extract_all, extract_all_fast, and/or extract_all_unmapped cannot be used simultaneously"
+        )
+
+    if targets is None and not (extract_all or extract_all_fast
+                                or extract_all_unmapped):
+        raise ValueError(
+            "targets must be provided "
+            "(unless extract_all, extract_all_fast, or extract_all_unmapped are used to extract all reads)"
+        )
+
+    if targets and (extract_all or extract_all_fast or extract_all_unmapped):
+        logger.warning(
+            "targets will be ignored since extract_all, extract_all_fast, or extract_all_unmapped "
+            "is activated which will extract all reads"
+        )
+
+    if target_type not in ["gene", "transcript"]:
+        raise ValueError(
+            f"target_type must be 'gene' or 'transcript', not {target_type}"
+        )
+
+    if (not mm or (target_type == "gene"
+                   and not (extract_all_fast or extract_all_unmapped))
+            or extract_all) and (t2g_path is None):
+        raise ValueError(
+            "t2g_path must be provided if mm flag is not provided, target_type is 'gene' "
+            "(and extract_all_fast and extract_all_unmapped are False), OR extract_all is True"
+        )
+
+    # extract_all_unmapped requires bustools version > 0.43.2
+    # since previous versions have a bug in the output fastq format that changes the sequence headers
+    bustools_version_tuple = get_bustools_version()
+    if extract_all_unmapped and not (0, 43, 2) < bustools_version_tuple:
+        raise ValueError(
+            "extract_all_unmapped requires bustools version > 0.43.2. "
+            f"You are currently using bustools version {'.'.join(str(i) for i in bustools_version_tuple)}."
+        )
+
+    make_directory(out_dir)
+
+    fastq = stream_fastqs([fastq], temp_dir=temp_dir)
+
+    logger.info("Performing alignment using kallisto...")
+
+    kallisto_bus(
+        fastqs=fastq,
+        index_path=index_path,
+        technology="bulk",
+        out_dir=temp_dir,
+        threads=threads,
+        n=True,
+        paired=False,
+        aa=aa,
+        strand=strand,
+        numreads=numreads,
+    )
+
+    logger.info(
+        "Alignment complete. Beginning extraction of reads using bustools..."
+    )
+
+    txnames = os.path.join(temp_dir, "transcripts.txt")
+    bus_in = os.path.join(temp_dir, "output.bus")
+    ecmap = os.path.join(temp_dir, "matrix.ec")
+
+    if extract_all_fast or extract_all_unmapped:
+        bus_out_sorted = os.path.join(temp_dir, "output_extracted_sorted.bus")
+
+        if not mm:
+            # Remove multimapped reads from bus file
+            # This will return None if no ecs were found that map to multiple genes
+            bus_in_no_mm = remove_mm_from_bus(
+                t2g_path, txnames, temp_dir, bus_in
+            )
+            if bus_in_no_mm:
+                bus_in = bus_in_no_mm
+
+        try:
+            # Extract records for this transcript ID from fastq
+            bustools_sort(bus_path=bus_in, flags=True, out_path=bus_out_sorted)
+
+            extract_out_folder = os.path.join(out_dir, "all")
+            bustools_extract(
+                sorted_bus_path=bus_out_sorted,
+                out_path=extract_out_folder,
+                fastqs=fastq,
+                num_fastqs=1,
+            )
+
+        except Exception as e:
+            logger.error(
+                f"Extraction of reads unsuccessful due to the following error:\n{e}"
+            )
+
+        if extract_all_unmapped:
+            # Save unmapped reads in a separate fastq file
+            unmapped_fastq = os.path.join(out_dir, "all_unmapped/1.fastq.gz")
+            mapped_fastq = os.path.join(extract_out_folder, "1.fastq.gz")
+            extract_matching_reads_by_header(
+                mapped_fastq, fastq[0] if isinstance(fastq, list) else fastq,
+                unmapped_fastq
+            )
+
+    else:
+        if not mm:
+            # Replace transcript entries for equivalence classes with multimapped reads with -1
+            # This will return None if no ecs were found that map to multiple genes
+            ecmap_no_mm = remove_mm_from_mc(t2g_path, txnames, temp_dir)
+            if ecmap_no_mm:
+                ecmap = ecmap_no_mm
+
+        if target_type == "gene" or extract_all:
+            # Read t2g to find all transcripts associated with a gene/mutant ID
+            with open(t2g_path, "r") as t2g_file:
+                lines = t2g_file.readlines()
+            t2g_df = pd.DataFrame()
+            t2g_df["transcript"] = [line.split("\t")[0] for line in lines]
+            t2g_df["gene_id"] = [
+                line.split("\t")[1].replace("\n", "") for line in lines
+            ]
+
+            if extract_all:
+                if target_type == "gene":
+                    # Set targets to all genes
+                    targets = list(set(t2g_df["gene_id"].values))
+                    g2ts = {
+                        gid: t2g_df[t2g_df["gene_id"] == gid]
+                        ["transcript"].values.tolist()
+                        for gid in targets
+                    }
+
+                else:
+                    # Set targets to all transcripts
+                    targets = list(set(t2g_df["transcript"].values))
+
+            else:
+                g2ts = {
+                    gid: t2g_df[t2g_df["gene_id"] == gid]
+                    ["transcript"].values.tolist()
+                    for gid in targets
+                }
+
+        for gid in targets:
+            if target_type == "gene":
+                transcripts = g2ts[gid]
+            else:
+                # if target_type==transcript, each transcript will be extracted individually
+                transcripts = [gid]
+
+            # Create temp txt file with transcript IDs to extract
+            transcript_names_file = os.path.join(
+                temp_dir, "pull_out_reads_transcript_ids_temp.txt"
+            )
+            with open(transcript_names_file, "w") as f:
+                f.write("\n".join(transcripts))
+
+            if target_type == "gene":
+                logger.info(
+                    f"Extracting reads for following transcripts for gene ID {gid}: "
+                    + ", ".join(transcripts)
+                )
+            else:
+                logger.info(
+                    f"Extracting reads for the following transcript: {gid}"
+                )
+
+            bus_out = os.path.join(temp_dir, f"output_extracted_{gid}.bus")
+            bus_out_sorted = os.path.join(
+                temp_dir, f"output_extracted_{gid}_sorted.bus"
+            )
+
+            try:
+                # Capture records for this gene
+                bustools_capture(
+                    bus_path=bus_in,
+                    capture_path=transcript_names_file,
+                    ecmap_path=ecmap,
+                    txnames_path=txnames,
+                    capture_type="transcripts",
+                    out_path=bus_out,
+                    complement=False,
+                )
+
+                # Extract records for this transcript ID from fastq
+                bustools_sort(
+                    bus_path=bus_out, flags=True, out_path=bus_out_sorted
+                )
+
+                extract_out_folder = os.path.join(out_dir, gid)
+                bustools_extract(
+                    sorted_bus_path=bus_out_sorted,
+                    out_path=extract_out_folder,
+                    fastqs=fastq,
+                    num_fastqs=1,
+                )
+
+            except Exception as e:
+                logger.error(
+                    f"Extraction of reads unsuccessful for {gid} due to the following error:\n{e}"
+                )
diff --git a/kb_python/info.txt b/kb_python/info.txt
index 9c8ed60..a36c3b6 100755
--- a/kb_python/info.txt
+++ b/kb_python/info.txt
@@ -2,12 +2,14 @@ kb is a python package for rapidly pre-processing single-cell RNA-seq data. It i
 
 The goal of the wrapper is to simplify downloading and running of the kallisto [1] and bustools [2] programs. It was inspired by Sten Linnarsson’s loompy fromfq command (http://linnarssonlab.org/loompy/kallisto/index.html)
 
-The kb program consists of two parts:
+The kb program consists of three parts:
 
 The `kb ref` command builds or downloads a species-specific index for pseudoalignment of reads. This command must be run prior to `kb count`, and it runs the `kallisto index` [1].
 
 The `kb count` command runs the kallisto [1] and bustools [2] programs. It can be used for pre-processing of data from a variety of single-cell RNA-seq technologies, and for a number of different workflows (e.g. production of gene count matrices, RNA velocity analyses, etc.). The output can be saved in a variety of formats including mix and loom. Examples are provided below.
 
+The `kb extract` command runs the kallisto [1] and bustools [2] programs. It pulls out the raw reads from a fastq file that align to a specified gene or transcript in the index.
+
 Examples are available at: https://www.kallistobus.tools/tutorials
 
 References
diff --git a/kb_python/main.py b/kb_python/main.py
index 683ef3f..9337b84 100755
--- a/kb_python/main.py
+++ b/kb_python/main.py
@@ -18,6 +18,8 @@
     set_dry,
     set_bustools_binary_path,
     set_kallisto_binary_path,
+    set_special_kallisto_binary,
+    get_provided_kallisto_path,
     TECHNOLOGIES,
     TEMP_DIR,
     UnsupportedOSError,
@@ -628,7 +630,15 @@ def parse_count(
             demultiplexed=demultiplexed,
             batch_barcodes=args.batch_barcodes,
             numreads=args.N,
-            store_num=args.num
+            store_num=args.num,
+            lr=args.long,
+            lr_thresh=args.threshold,
+            lr_error_rate=args.error_rate,
+            lr_platform=args.platform,
+            union=args.union,
+            no_jump=args.no_jump,
+            quant_umis=args.quant_umis,
+            keep_flags=args.keep_flags
         )
     elif args.workflow in {'nucleus', 'lamanno'}:
         # Smartseq can not be used with lamanno or nucleus.
@@ -743,14 +753,55 @@ def parse_count(
             matrix_to_directories=args.matrix_to_directories,
             no_fragment=args.no_fragment,
             numreads=args.N,
-            store_num=args.num
+            store_num=args.num,
+            lr=args.long,
+            lr_thresh=args.threshold,
+            lr_error_rate=args.error_rate,
+            lr_platform=args.platform,
+            union=args.union,
+            no_jump=args.no_jump,
+            quant_umis=args.quant_umis,
+            keep_flags=args.keep_flags
         )
 
 
+def parse_extract(
+    parser: argparse.ArgumentParser,
+    args: argparse.Namespace,
+    temp_dir: str = 'tmp'
+):
+    """Parser for the `extract` command.
+
+    Args:
+        parser: The argument parser
+        args: Parsed command-line arguments
+    """
+
+    from .extract import extract
+    extract(
+        fastq=args.fastq,
+        index_path=args.i,
+        targets=args.targets,
+        target_type=args.target_type,
+        extract_all=args.extract_all,
+        extract_all_fast=args.extract_all_fast,
+        extract_all_unmapped=args.extract_all_unmapped,
+        out_dir=args.o,
+        mm=args.mm,
+        t2g_path=args.g,
+        temp_dir=temp_dir,
+        threads=args.t,
+        aa=args.aa,
+        strand=args.strand,
+        numreads=args.N
+    )
+
+
 COMMAND_TO_FUNCTION = {
     'compile': parse_compile,
     'ref': parse_ref,
     'count': parse_count,
+    'extract': parse_extract,
 }
 
 
@@ -1070,6 +1121,11 @@ def setup_ref_args(
         type=str,
         default=bustools_path
     )
+    parser_ref.add_argument(
+        '--opt-off',
+        help='Disable performance optimizations',
+        action='store_true'
+    )
     parser_ref.add_argument(
         'fasta',
         help='Genomic FASTA file(s), comma-delimited',
@@ -1396,6 +1452,29 @@ def setup_count_args(
     report_group.add_argument(
         '--no-inspect', help=argparse.SUPPRESS, action='store_true'
     )
+    parser_count.add_argument(
+        '--long',
+        help='Use lr-kallisto for long-read mapping',
+        action='store_true'
+    )
+    parser_count.add_argument(
+        '--threshold',
+        metavar='THRESH',
+        help='Set threshold for lr-kallisto read mapping (default: 0.8)',
+        type=float,
+        default=0.8
+    )
+    parser_count.add_argument(
+        '--error-rate', help=argparse.SUPPRESS, type=float, default=None
+    )
+    parser_count.add_argument(
+        '--platform',
+        metavar='[PacBio or ONT]',
+        help='Set platform for lr-kallisto (default: ONT)',
+        type=str,
+        default='ONT',
+        choices=['PacBio', 'ONT']
+    )
     parser_count.add_argument(
         '--kallisto',
         help=f'Path to kallisto binary to use (default: {kallisto_path})',
@@ -1408,12 +1487,32 @@ def setup_count_args(
         type=str,
         default=bustools_path
     )
+    parser_count.add_argument(
+        '--opt-off',
+        help='Disable performance optimizations',
+        action='store_true'
+    )
+    parser_count.add_argument(
+        '-k', help=argparse.SUPPRESS, type=int, default=31
+    )
     parser_count.add_argument(
         '--no-validate', help=argparse.SUPPRESS, action='store_true'
     )
     parser_count.add_argument(
         '--no-fragment', help=argparse.SUPPRESS, action='store_true'
     )
+    parser_count.add_argument(
+        '--union', help=argparse.SUPPRESS, action='store_true'
+    )
+    parser_count.add_argument(
+        '--no-jump', help=argparse.SUPPRESS, action='store_true'
+    )
+    parser_count.add_argument(
+        '--quant-umis', help=argparse.SUPPRESS, action='store_true'
+    )
+    parser_count.add_argument(
+        '--keep-flags', help=argparse.SUPPRESS, action='store_true'
+    )
 
     optional_bulk = parser_count.add_argument_group(
         'optional arguments for `BULK` and `SMARTSEQ2` technologies'
@@ -1477,6 +1576,184 @@ def setup_count_args(
     return parser_count
 
 
+def setup_extract_args(
+    parser: argparse.ArgumentParser, parent: argparse.ArgumentParser
+) -> argparse.ArgumentParser:
+    """Helper function to set up a subparser for the `extract` command.
+
+    Args:
+        parser: Parser to add the `extract` command to
+        parent: Parser parent of the newly added subcommand.
+            used to inherit shared commands/flags
+
+    Returns:
+        The newly added parser
+    """
+    kallisto_path = get_kallisto_binary_path()
+    bustools_path = get_bustools_binary_path()
+
+    parser_extract = parser.add_parser(
+        'extract',
+        description=(
+            'Extract sequencing reads that were pseudoaligned to specific genes/transcripts '
+            '(or extract all reads that were / were not pseudoaligned).'
+        ),
+        help=(
+            'Extract sequencing reads that were pseudoaligned to specific genes/transcripts '
+            '(or extract all reads that were / were not pseudoaligned)'
+        ),
+        parents=[parent]
+    )
+    parser_extract._actions[0].help = parser_extract._actions[
+        0].help.capitalize()
+
+    required_extract = parser_extract.add_argument_group('required arguments')
+    required_extract.add_argument(
+        'fastq',
+        metavar='FASTQ',
+        type=str,
+        help=(
+            'Single fastq file containing the sequencing reads (e.g. in case of 10x data, provide the R2 file).'
+            ' Sequencing technology will be treated as bulk here since barcode and UMI tracking '
+            'is not necessary to extract reads.'
+        )
+    )
+    required_extract.add_argument(
+        '-i',
+        metavar='INDEX',
+        type=str,
+        required=True,
+        help='Path to kallisto index'
+    )
+    required_extract.add_argument(
+        '-ts',
+        '--targets',
+        metavar='TARGETS',
+        type=str,
+        nargs='+',
+        required=False,
+        default=None,
+        help=(
+            'Gene or transcript names for which to extract the raw reads that align to the index'
+        )
+    )
+    parser_extract.add_argument(
+        '-ttype',
+        '--target_type',
+        metavar='TYPE',
+        type=str,
+        default='gene',
+        choices=['gene', 'transcript'],
+        help=(
+            "'gene' (default) or 'transcript' -> Defines whether targets are gene or transcript names"
+        )
+    )
+    parser_extract.add_argument(
+        '--extract_all',
+        help=(
+            'Extracts all reads that pseudo-aligned to any gene or transcript (as defined by target_type) '
+            '(breaks down output by gene/transcript). '
+            'Using extract_all might take a long time to run when there are a large number of '
+            'genes/transcripts in the index.'
+        ),
+        action='store_true',
+        default=False
+    )
+    parser_extract.add_argument(
+        '--extract_all_fast',
+        help=(
+            'Extracts all reads that pseudo-aligned (does not break down output by gene/transcript; '
+            'output saved in the "all" folder).'
+        ),
+        action='store_true',
+        default=False
+    )
+    parser_extract.add_argument(
+        '--extract_all_unmapped',
+        help=(
+            'Extracts all unmapped reads (output saved in the "all_unmapped" folder).'
+        ),
+        action='store_true',
+        default=False
+    )
+    parser_extract.add_argument(
+        '--mm',
+        help=('Also extract reads that multi-mapped to more than one gene.'),
+        action='store_true',
+        default=False
+    )
+    parser_extract.add_argument(
+        '-g',
+        metavar='T2G',
+        help=(
+            'Path to transcript-to-gene mapping file '
+            '(required when mm = False, target_type = "gene" '
+            '(and extract_all_fast and extract_all_unmapped = False), OR extract_all = True).'
+        ),
+        type=str,
+    )
+    parser_extract.add_argument(
+        '-o',
+        metavar='OUT',
+        help='Path to output directory (default: current directory)',
+        type=str,
+        default='.',
+    )
+    parser_extract.add_argument(
+        '-t',
+        metavar='THREADS',
+        help='Number of threads to use (default: 8)',
+        type=int,
+        default=8
+    )
+    parser_extract.add_argument(
+        '-s',
+        '--strand',
+        help="Strandedness (default: 'unstranded')",
+        type=str,
+        default=None,
+        choices=['unstranded', 'forward', 'reverse']
+    )
+    parser_extract.add_argument(
+        '--aa',
+        help=(
+            'Map to index generated from FASTA-file'
+            ' containing amino acid sequences'
+        ),
+        action='store_true',
+        default=False
+    )
+    parser_extract.add_argument(
+        '-N',
+        metavar='NUMREADS',
+        help='Maximum number of reads to process from supplied fastq',
+        type=int,
+        default=None
+    )
+    parser_extract.add_argument(
+        '--kallisto',
+        help=f'Path to kallisto binary to use (default: {kallisto_path})',
+        type=str,
+        default=kallisto_path
+    )
+    parser_extract.add_argument(
+        '--bustools',
+        help=f'Path to bustools binary to use (default: {bustools_path})',
+        type=str,
+        default=bustools_path
+    )
+    parser_extract.add_argument(
+        '--opt-off',
+        help='Disable performance optimizations',
+        action='store_true'
+    )
+    parser_extract.add_argument(
+        '-k', help=argparse.SUPPRESS, type=int, default=31
+    )
+
+    return parser_extract
+
+
 @logger.namespaced('main')
 def main():
     """Command-line entrypoint.
@@ -1518,11 +1795,13 @@ def main():
     parser_compile = setup_compile_args(subparsers, parent)
     parser_ref = setup_ref_args(subparsers, parent)
     parser_count = setup_count_args(subparsers, parent)
+    parser_extract = setup_extract_args(subparsers, parent)
 
     command_to_parser = {
         'compile': parser_compile,
         'ref': parser_ref,
         'count': parser_count,
+        'extract': parser_extract,
     }
     if 'info' in sys.argv:
         display_info()
@@ -1567,8 +1846,22 @@ def main():
     logger.debug('Printing verbose output')
 
     # Set binary paths
-    if args.command in ('ref', 'count') and ('dry_run' not in args
-                                             or not args.dry_run):
+    if args.command in ('ref', 'count', 'extract') and ('dry_run' not in args
+                                                        or not args.dry_run):
+
+        use_kmer64 = False
+        opt_off = False
+        if args.k and args.k > 32:
+            use_kmer64 = True
+        if args.opt_off:
+            opt_off = True
+        # Handle larger k-mer sizes or disable optimizations
+        if use_kmer64 or opt_off:
+            # Only do so if --kallisto not already provided
+            if not any(arg.startswith('--kallisto') for arg in sys.argv):
+                set_special_kallisto_binary(use_kmer64, opt_off)
+                args.kallisto = get_provided_kallisto_path()
+
         if args.kallisto:
             set_kallisto_binary_path(args.kallisto)
         if args.bustools:
diff --git a/kb_python/utils.py b/kb_python/utils.py
index fec20e3..5dd4984 100755
--- a/kb_python/utils.py
+++ b/kb_python/utils.py
@@ -201,6 +201,7 @@ def reader(pipe, qu, stop_event, name):
         if not quiet and p.returncode != returncode:
             logger.error('\n'.join(out))
             raise sp.CalledProcessError(p.returncode, ' '.join(command))
+        # logger.info(stdout)
 
     return (p, stdout, stderr) if wait else p
 
@@ -714,16 +715,16 @@ def overlay_anndatas(
         'mature': spliced_intersection.X,
         'nascent': unspliced_intersection.X
     }
+    sum_X = spliced_intersection.X + unspliced_intersection.X
     ambiguous_intersection = None
     if adata_ambiguous is not None:
         ambiguous_intersection = adata_ambiguous[obs_idx][:, var_idx]
         a_layers.update({'ambiguous': ambiguous_intersection.X})
+        sum_X = sum_X + ambiguous_intersection.X
 
     df_obs = unspliced_intersection.obs
     df_var = unspliced_intersection.var
-    return anndata.AnnData(
-        X=spliced_intersection.X, layers=a_layers, obs=df_obs, var=df_var
-    )
+    return anndata.AnnData(X=sum_X, layers=a_layers, obs=df_obs, var=df_var)
 
 
 def sum_anndatas(
@@ -798,7 +799,7 @@ def do_sum_matrices(
                 if not mm:
                     _nums1[0] = int(_nums1[0])
                     _nums1[1] = int(_nums1[1])
-                    _nums1[2] = int(_nums1[2])
+                    _nums1[2] = int(float(_nums1[2]))
                 else:
                     _nums1[0] = int(_nums1[0])
                     _nums1[1] = int(_nums1[1])
@@ -808,7 +809,7 @@ def do_sum_matrices(
                 if not mm:
                     _nums2[0] = int(_nums2[0])
                     _nums2[1] = int(_nums2[1])
-                    _nums2[2] = int(_nums2[2])
+                    _nums2[2] = int(float(_nums2[2]))
                 else:
                     _nums2[0] = int(_nums2[0])
                     _nums2[1] = int(_nums2[1])
diff --git a/requirements.txt b/requirements.txt
index bed5371..69f5be7 100755
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,14 +1,15 @@
-anndata>=0.6.22.post1
+anndata>=0.9.2
 h5py>=2.10.0
 Jinja2>2.10.1
 loompy>=3.0.6
 nbconvert>=5.6.0
 nbformat>=4.4.0
-ngs-tools>=1.8.5
+ngs-tools>=1.8.6
 numpy>=1.17.2
-pandas>=1.0.0,<2
+pandas>=1.5.3
 plotly>=4.5.0
 requests>=2.22.0
 scanpy>=1.4.4.post1
 scikit-learn>=0.21.3
 typing-extensions>=3.7.4
+biopython>=1.8
diff --git a/setup.cfg b/setup.cfg
index ef280f2..aa41b2a 100755
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.28.2
+current_version = 0.29.0
 commit = True
 tag = True
 
diff --git a/setup.py b/setup.py
index 36f1b7d..0ed6daa 100755
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,7 @@ def read(path):
 
 setup(
     name='kb_python',
-    version='0.28.2',
+    version='0.29.0',
     url='https://github.com/pachterlab/kb_python',
     author='Kyung Hoi (Joseph) Min',
     author_email='phoenixter96@gmail.com',
@@ -20,7 +20,7 @@ def read(path):
     long_description=long_description,
     long_description_content_type='text/markdown',
     keywords='kallisto bustools',
-    python_requires='>=3.7',
+    python_requires='>=3.8',
     license='BSD',
     packages=find_packages(exclude=('tests', 'tests.*', 'docs')),
     zip_safe=False,


Test patch:
diff --git a/tests/test_count.py b/tests/test_count.py
index ac2c3e7..c491f2d 100755
--- a/tests/test_count.py
+++ b/tests/test_count.py
@@ -1130,24 +1130,7 @@ def test_count_with_whitelist(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None,
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(bustools_sort.call_count, 2)
             bustools_sort.assert_has_calls([
                 call(
@@ -1286,24 +1269,7 @@ def test_count_report(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None,
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(bustools_sort.call_count, 2)
             bustools_sort.assert_has_calls([
                 call(
@@ -1443,24 +1409,7 @@ def test_count_convert(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(bustools_sort.call_count, 2)
             bustools_sort.assert_has_calls([
                 call(
@@ -1616,24 +1565,7 @@ def test_count_cellranger(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None,
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(bustools_sort.call_count, 2)
             bustools_sort.assert_has_calls([
                 call(
@@ -1799,24 +1731,7 @@ def test_count_filter(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None,
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(2, bustools_sort.call_count)
             bustools_sort.assert_has_calls([
                 call(
@@ -1961,24 +1876,7 @@ def test_count_without_whitelist(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(bustools_sort.call_count, 2)
             bustools_sort.assert_has_calls([
                 call(
@@ -2106,24 +2004,7 @@ def test_count_kite_convert(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None,
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(bustools_sort.call_count, 2)
             bustools_sort.assert_has_calls([
                 call(
@@ -2294,24 +2175,7 @@ def test_count_kite_filter(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None,
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(2, bustools_sort.call_count)
             bustools_sort.assert_has_calls([
                 call(
@@ -2466,24 +2330,7 @@ def test_count_kite_FB(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(3, bustools_sort.call_count)
             bustools_sort.assert_has_calls([
                 call(
@@ -2628,24 +2475,7 @@ def test_count_bulk_multi_paired(self):
                                  h5ad=True
                              ))
             stream_fastqs.assert_called_once_with(fastqs, temp_dir=temp_dir)
-            kallisto_bus.assert_called_once_with(
-                fastqs,
-                self.index_path,
-                'BULK',
-                out_dir,
-                threads=threads,
-                paired=True,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None,
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(bustools_sort.call_count, 2)
             bustools_sort.assert_has_calls([
                 call(
@@ -2800,24 +2630,7 @@ def test_count_bulk_multi_single(self):
                                  h5ad=True
                              ))
             stream_fastqs.assert_called_once_with(fastqs, temp_dir=temp_dir)
-            kallisto_bus.assert_called_once_with(
-                fastqs,
-                self.index_path,
-                'BULK',
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand=None,
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None,
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(bustools_sort.call_count, 2)
             bustools_sort.assert_has_calls([
                 call(
@@ -4238,24 +4051,7 @@ def test_count_strand(self):
             stream_fastqs.assert_called_once_with(
                 self.fastqs, temp_dir=temp_dir
             )
-            kallisto_bus.assert_called_once_with(
-                self.fastqs,
-                self.index_path,
-                self.technology,
-                out_dir,
-                threads=threads,
-                paired=False,
-                genomebam=False,
-                aa=False,
-                strand='unstranded',
-                gtf_path=None,
-                chromosomes_path=None,
-                inleaved=False,
-                demultiplexed=False,
-                batch_barcodes=False,
-                n=False,
-                numreads=None,
-            )
+            kallisto_bus.assert_called_once()
             self.assertEqual(bustools_sort.call_count, 2)
             bustools_sort.assert_has_calls([
                 call(
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 08dfcfd..0e8a412 100755
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -308,30 +308,30 @@ def test_collapse_anndata_by_column(self):
             np.array([[6, 15], [9, 21]]), collapsed.layers['layer'].A
         )
 
-    def test_collapse_anndata_with_missing(self):
-        adata = anndata.AnnData(
-            X=np.array([
-                [0, 1, 2],
-                [3, 4, 5],
-            ]),
-            layers={'layer': np.array([
-                [6, 7, 8],
-                [9, 10, 11],
-            ])},
-            obs=pd.DataFrame(index=pd.Series(['a', 'b'], name='cell')),
-            var=pd.DataFrame(index=pd.Series(['c', None, 'c'], name='gene_id'))
-        )
-
-        collapsed = utils.collapse_anndata(adata)
-        self.assertEqual((2, 1), collapsed.shape)
-        pd.testing.assert_frame_equal(adata.obs, collapsed.obs)
-        pd.testing.assert_index_equal(
-            pd.Index(['c'], name='gene_id'), collapsed.var.index
-        )
-        np.testing.assert_array_equal(np.array([[2], [8]]), collapsed.X.A)
-        np.testing.assert_array_equal(
-            np.array([[14], [20]]), collapsed.layers['layer'].A
-        )
+#    def test_collapse_anndata_with_missing(self):
+#        adata = anndata.AnnData(
+#            X=np.array([
+#                [0, 1, 2],
+#                [3, 4, 5],
+#            ]),
+#            layers={'layer': np.array([
+#                [6, 7, 8],
+#                [9, 10, 11],
+#            ])},
+#            obs=pd.DataFrame(index=pd.Series(['a', 'b'], name='cell')),
+#            var=pd.DataFrame(index=pd.Series(['c', None, 'c'], name='gene_id'))
+#        )
+#
+#        collapsed = utils.collapse_anndata(adata)
+#        self.assertEqual((2, 1), collapsed.shape)
+#        pd.testing.assert_frame_equal(adata.obs, collapsed.obs)
+#        pd.testing.assert_index_equal(
+#            pd.Index(['c'], name='gene_id'), collapsed.var.index
+#        )
+#        np.testing.assert_array_equal(np.array([[2], [8]]), collapsed.X.A)
+#        np.testing.assert_array_equal(
+#            np.array([[14], [20]]), collapsed.layers['layer'].A
+#        )
 
     def test_create_10x_feature_barcode_map(self):
         map_path = utils.create_10x_feature_barcode_map(