-
Notifications
You must be signed in to change notification settings - Fork 0
/
filtercells.py
70 lines (54 loc) · 2.17 KB
/
filtercells.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#!/usr/bin/env python3
"""
Simple cell filtering (keeps only cells with a certain minimum UMI count)
"""
from argparse import ArgumentParser
from pathlib import Path
import sys
import shutil
import scipy
from xopen import xopen
def main():
parser = ArgumentParser()
parser.add_argument("--umis", type=int, default=200, help="Minimum UMI count")
parser.add_argument("raw", type=Path, help="Input directory with matrix.mtx, features.tsv and barcodes.tsv file")
parser.add_argument("target", type=Path, help="Output directory (must not exist)")
args = parser.parse_args()
log("Reading barcodes.tsv(.gz) ...")
try:
with xopen(args.raw / "barcodes.tsv") as f:
barcodes = list(f)
barcodes_name = "barcodes.tsv"
except FileNotFoundError:
with xopen(args.raw / "barcodes.tsv.gz") as f:
barcodes = list(f)
barcodes_name = "barcodes.tsv.gz"
log("Reading matrix.mtx(.gz) ...")
try:
matrix = scipy.io.mmread(args.raw / "matrix.mtx.gz")
matrix_name = "matrix.mtx.gz"
except FileNotFoundError:
matrix = scipy.io.mmread(args.raw / "matrix.mtx")
matrix_name = "matrix.mtx"
log(f"Input matrix: {matrix.shape[1]} cell barcodes; {matrix.shape[0]} features")
log("Filtering ...")
matrix = matrix.tocsr()
umi_counts = matrix.sum(axis=0).A1
mask = umi_counts >= args.umis
filtered_matrix = matrix[:, mask]
filtered_barcodes = [barcode for barcode, use_it in zip(barcodes, mask) if use_it]
log(f"Filtered matrix: {filtered_matrix.shape[1]} cell barcodes; {filtered_matrix.shape[0]} features")
log("Writing output ...")
args.target.mkdir(exist_ok=True)
with xopen(args.target / barcodes_name, mode="w") as f:
for barcode in filtered_barcodes:
f.write(barcode)
try:
shutil.copy(args.raw / "features.tsv.gz", args.target / "features.tsv.gz")
except FileNotFoundError:
shutil.copy(args.raw / "features.tsv", args.target / "features.tsv")
scipy.io.mmwrite(args.target / matrix_name, filtered_matrix)
def log(*args, **kwargs):
print(*args, **kwargs, file=sys.stderr)
if __name__ == "__main__":
main()