diff --git a/.github/workflows/pytest-workflow.yml b/.github/workflows/pytest-workflow.yml index 91b8199..80e4cd7 100644 --- a/.github/workflows/pytest-workflow.yml +++ b/.github/workflows/pytest-workflow.yml @@ -19,7 +19,7 @@ jobs: name: test strategy: matrix: - python-version: ['3.9', '3.10', '3.11', '3.12', '3.13'] + python-version: ['3.9', '3.10', '3.11', '3.12', '3.13', '3.14'] steps: - name: Checkout uses: actions/checkout@v3.3.0 diff --git a/irescue/_version.py b/irescue/_version.py index 07e1f21..2aa5cc3 100644 --- a/irescue/_version.py +++ b/irescue/_version.py @@ -1 +1 @@ -__version__ = '1.2.0b2' +__version__ = "1.2.0b3" diff --git a/irescue/count.py b/irescue/count.py index 9a06f60..2f544f2 100644 --- a/irescue/count.py +++ b/irescue/count.py @@ -1,35 +1,41 @@ #!/usr/bin/env python +import gzip +import os from collections import Counter from itertools import combinations -import numpy as np + import networkx as nx -from irescue.misc import get_ranges, getlen, writerr, run_shell_cmd -from irescue.network import build_substr_idx, gen_ec_pairs +import numpy as np +from scipy.sparse import lil_matrix + from irescue.em import run_em -import gzip -import os +from irescue.misc import get_ranges, getlen, run_shell_cmd, writerr +from irescue.network import build_substr_idx, gen_ec_pairs + class EquivalenceClass: def __init__( - self, - index: int, - umi: bytes, - features: set, - count: int + self, index: int, umi: bytes, features: set, count: int ) -> None: self.index = index self.umi = umi self.features = features self.count = count + def to_tuple(self): return (self.umi, self.features, self.count) + def hdist(self, umi): return sum(1 for i, j in zip(self.umi, umi) if i != j) + def connect(self, eqc, threshold): - return (self.count >= (2 * eqc.count) - 1 - and self.features.intersection(eqc.features) - and self.hdist(eqc.umi) <= threshold) + return ( + self.count >= (2 * eqc.count) - 1 + and self.features.intersection(eqc.features) + and self.hdist(eqc.umi) <= threshold + ) + def pathfinder(graph, node, path=[], features=None): """ @@ -37,22 +43,26 @@ def pathfinder(graph, node, path=[], features=None): a starting node. Can be used iteratively to find all possible paths. """ if not features: - features = graph.nodes[node]['ft'] + features = graph.nodes[node]["ft"] path += [node] for next_node in graph.successors(node): - if (features.intersection(graph.nodes[next_node]['ft']) - and next_node not in path): + if ( + features.intersection(graph.nodes[next_node]["ft"]) + and next_node not in path + ): path = pathfinder(graph, next_node, path, features) return path + def index_features(features_file): idx = {} - with gzip.open(features_file, 'rb') as f: + with gzip.open(features_file, "rb") as f: for i, line in enumerate(f, start=1): - ft = line.strip().split(b'\t')[0] + ft = line.strip().split(b"\t")[0] idx[ft] = i return idx + def parse_maps(maps_file, feature_index): """ maps_file : str @@ -61,17 +71,17 @@ def parse_maps(maps_file, feature_index): CB, [(UMI , {FT , ...} , count ) , ...] """ - with gzip.open(maps_file, 'rb') as f: - cb, umi, feat, count = f.readline().strip().split(b'\t') + with gzip.open(maps_file, "rb") as f: + cb, umi, feat, count = f.readline().strip().split(b"\t") i = 0 it = cb count = int(count) - feat = {feature_index[ft] for ft in feat.split(b',')} + feat = {feature_index[ft] for ft in feat.split(b",")} eqcl = [EquivalenceClass(i, umi, feat, count)] for line in f: - cb, umi, feat, count = line.strip().split(b'\t') + cb, umi, feat, count = line.strip().split(b"\t") count = int(count) - feat = {feature_index[ft] for ft in feat.split(b',')} + feat = {feature_index[ft] for ft in feat.split(b",")} if cb == it: i += 1 eqcl.append(EquivalenceClass(i, umi, feat, count)) @@ -82,8 +92,10 @@ def parse_maps(maps_file, feature_index): eqcl = [EquivalenceClass(i, umi, feat, count)] yield it, eqcl -def compute_cell_counts(equivalence_classes, features_index, max_iters, - tolerance, dumpEC, no_umi): + +def compute_cell_counts( + equivalence_classes, features_index, max_iters, tolerance, dumpEC, no_umi +): """ Calculate TE counts of a single cell, given a list of equivalence classes. @@ -100,17 +112,20 @@ def compute_cell_counts(equivalence_classes, features_index, max_iters, """ # initialize TE counts and dedup log counts = Counter() - em_array = [] + em_array_rows = {} dump = {} if dumpEC else None number_of_features = len(features_index) + # With UMIs (10X-like datasets) if not no_umi: # build cell-wide UMI deduplication graph graph = nx.DiGraph() # add nodes with annotated features graph.add_nodes_from( - [(x.index, {'ft': x.features, 'c': x.count}) - for x in equivalence_classes] + [ + (x.index, {"ft": x.features, "c": x.count}) + for x in equivalence_classes + ] ) # make an iterator of umi pairs if len(equivalence_classes) > 25: @@ -129,9 +144,11 @@ def compute_cell_counts(equivalence_classes, features_index, max_iters, # collect graph metadata in a dictionary dump = {i: equivalence_classes[i].to_tuple() for i in graph.nodes} # split cell-wide graph into subgraphs of connected nodes - subgraphs = [graph.subgraph(x) for x in - nx.connected_components(graph.to_undirected())] - + subgraphs = [ + graph.subgraph(x) + for x in nx.connected_components(graph.to_undirected()) + ] + # solve UMI deduplication for each subgraph of connected nodes for subg in subgraphs: # find all parent nodes in graph @@ -140,7 +157,9 @@ def compute_cell_counts(equivalence_classes, features_index, max_iters, # if no parents are found due to bidirected edges, take all nodes # and the union of all features (i.e. all nodes are parents). parents = list(subg.nodes) - features = [list(set.union(*[subg.nodes[x]['ft'] for x in subg]))] + features = [ + tuple(set.union(*[subg.nodes[x]["ft"] for x in subg])) + ] else: # if parents node are found, features will be determined below. features = None @@ -150,7 +169,7 @@ def compute_cell_counts(equivalence_classes, features_index, max_iters, # find paths starting from each parent node for parent in parents: # populate this list with nodes utilized in paths - blacklist = [] + blacklist = set() # find paths in list of nodes starting from parent path = [] subg_copy = subg.copy() @@ -159,20 +178,23 @@ def compute_cell_counts(equivalence_classes, features_index, max_iters, # make a copy of subgraph and remove nodes already used # in a path if node not in blacklist: - path = pathfinder(subg_copy, node, path=[], features=None) + path = pathfinder( + subg_copy, node, path=[], features=None + ) for x in path: - blacklist.append(x) + blacklist.add(x) subg_copy.remove_node(x) paths[parent].append(path) # find the path configuration leading to the minimum number of # deduplicated UMIs -> list of lists of nodes path_config = [ - paths[k] for k, v in paths.items() + paths[k] + for k, v in paths.items() if len(v) == min([len(x) for x in paths.values()]) ][0] if not features: # take features from parent node of selected path configuration - features = [list(subg.nodes[x[0]]['ft']) for x in path_config] + features = [tuple(subg.nodes[x[0]]["ft"]) for x in path_config] else: # if features was already determined (i.e. no parent nodes), # multiplicate the feature's list by the number of paths @@ -183,65 +205,75 @@ def compute_cell_counts(equivalence_classes, features_index, max_iters, if len(feats) == 1: counts[feats[0]] += 1.0 elif len(feats) > 1: - row = [1 if x in feats else 0 - for x in range(1, number_of_features+1)] - em_array.append(row) + em_array_rows[len(em_array_rows)] = feats else: writerr(nx.to_dict_of_lists(subg)) - writerr([subg.nodes[x]['ft'] for x in subg.nodes]) - writerr([subg.nodes[x]['c'] for x in subg.nodes]) + writerr([subg.nodes[x]["ft"] for x in subg.nodes]) + writerr([subg.nodes[x]["c"] for x in subg.nodes]) writerr(path_config) writerr(path) writerr(features) writerr(feats) - writerr("Error: no common features detected in subgraph's" - " path.", error=True) + writerr( + "Error: no common features detected in subgraph's" + " path.", + error=True, + ) # add EC log to dump if dumpEC: for i, path_ in enumerate(path_config): # add empty fields to parent node parent_ = path_[0] path_.pop(0) - dump[parent_] += (b'', b'') + dump[parent_] += (b"", b"") # if child nodes are present, add parent node informations for x in path_: # add parent's UMI sequence and dedup features dump[x] += (dump[parent_][0], features[i]) - else: ## in case of UMI-less + # Without UMIs (Smart-seq-like datasets) + else: for eqc in equivalence_classes: feats = list(eqc.features) if len(feats) == 1: counts[feats[0]] += 1.0 else: - row = [1 if x in feats else 0 - for x in range(1, number_of_features+1)] - em_array.append(row) + em_array_rows[len(em_array_rows)] = feats if dumpEC: dump[eqc.index] = eqc.to_tuple() - + # EM stats placeholder in case of no multimapped UMIs em_stats = (None, None, None, None) - if em_array: + em_array = None + + if em_array_rows: # optimize the assignment of UMI from multimapping reads - em_array = np.array(em_array) + em_array = lil_matrix( + (len(em_array_rows), number_of_features), dtype=np.uint8 + ) + + for i, feats in em_array_rows.items(): + em_array.rows[i] = [feat_idx - 1 for feat_idx in feats] + em_array.data[i] = [1] * len(em_array.rows[i]) + + em_array = em_array.tocsr() + # save an array with features > 0, as in em_array order - tokeep = np.argwhere(np.any(em_array[..., :] > 0, axis=0))[:,0] + 1 + tokeep = np.flatnonzero(em_array.sum(axis=0)) # remove unmapped features from em_array - todel = np.argwhere(np.all(em_array[..., :] == 0, axis=0)) - em_array = np.delete(em_array, todel, axis=1) + em_array = em_array[:, tokeep] # run EM em_counts, em_stats = run_em( - em_array, - cycles=max_iters, - tolerance=tolerance + em_array, cycles=max_iters, tolerance=tolerance ) - em_counts = [x*em_array.shape[0] for x in em_counts] - for i, c in zip(tokeep, em_counts): + em_counts = em_counts * em_array.shape[0] + + for i, c in zip(tokeep + 1, em_counts): if c > 0: counts[i] += c return dict(counts), dump, em_stats + def split_barcodes(barcodes_file, n): """ barcodes_file : iterable @@ -250,27 +282,39 @@ def split_barcodes(barcodes_file, n): out : int, dict """ nBarcodes = getlen(barcodes_file) - with gzip.open(barcodes_file, 'rb') as f: + with gzip.open(barcodes_file, "rb") as f: for i, chunk in enumerate(get_ranges(nBarcodes, n)): - yield i, {next(f).strip(): x+1 for x in chunk} + yield i, {next(f).strip(): x + 1 for x in chunk} -def run_count(maps_file, features_index, tmpdir, no_umi, dumpEC, max_iters, - tolerance, verbose, barcodes_set): + +def run_count( + maps_file, + features_index, + tmpdir, + no_umi, + dumpEC, + max_iters, + tolerance, + verbose, + barcodes_set, +): # NB: keep args order consistent with main.countFun taskn, barcodes = barcodes_set - matrix_file = os.path.join(tmpdir, f'{taskn}_matrix.mtx.gz') - dump_file = os.path.join(tmpdir, f'{taskn}_EqCdump.tsv.gz') - with gzip.open(matrix_file, 'wb') as f, \ - gzip.open(dump_file, 'wb') if dumpEC \ - else gzip.open(os.devnull) as df: + matrix_file = os.path.join(tmpdir, f"{taskn}_matrix.mtx.gz") + dump_file = os.path.join(tmpdir, f"{taskn}_EqCdump.tsv.gz") + with ( + gzip.open(matrix_file, "wb") as f, + gzip.open(dump_file, "wb") if dumpEC else gzip.open(os.devnull) as df, + ): for cellbarcode, cellmaps in parse_maps(maps_file, features_index): if cellbarcode not in barcodes: continue cellidx = barcodes[cellbarcode] writerr( - f'[{taskn}] Run count for cell ' - f'{cellidx} ({cellbarcode.decode()})', - level=2, send=verbose + f"[{taskn}] Run count for cell " + f"{cellidx} ({cellbarcode.decode()})", + level=2, + send=verbose, ) cellcounts, dump, em_stats = compute_cell_counts( equivalence_classes=cellmaps, @@ -278,103 +322,132 @@ def run_count(maps_file, features_index, tmpdir, no_umi, dumpEC, max_iters, max_iters=max_iters, tolerance=tolerance, dumpEC=dumpEC, - no_umi=no_umi + no_umi=no_umi, ) writerr( f"[{taskn}] Write cell {cellidx} ({cellbarcode.decode()}). " f"EM cycles: {em_stats[0]}. Converged: {em_stats[1]}. " f"Log likelihood: {em_stats[2]}. Increment: {em_stats[3]}.", - level=1, send=verbose + level=1, + send=verbose, ) # round counts to 3rd decimal point and write to matrix file # only if count is at least 0.001 - lines = [f'{feature} {cellidx} {round(count, 3)}\n'.encode() - for feature, count in cellcounts.items() - if count >= 0.001] + lines = [ + f"{feature} {cellidx} {round(count, 3)}\n".encode() + for feature, count in cellcounts.items() + if count >= 0.001 + ] f.writelines(lines) if dumpEC: writerr( - f'[{taskn}] Write ECdump for cell ' - f'{cellidx} ({cellbarcode.decode()})', - level=1, send=verbose + f"[{taskn}] Write ECdump for cell " + f"{cellidx} ({cellbarcode.decode()})", + level=1, + send=verbose, ) # reverse features index to get names back - findex = dict(zip(features_index.values(), - features_index.keys())) - + findex = dict( + zip(features_index.values(), features_index.keys()) + ) + if not no_umi: dumplines = [ - b'\t'.join( - [str(cellidx).encode(), - cellbarcode, - str(i).encode(), + b"\t".join( + [ + str(cellidx).encode(), + cellbarcode, + str(i).encode(), + umi, + b",".join([findex[f] for f in feats]), + str(count).encode(), + pumi, + b",".join([findex[f] for f in pfeats]), + ] + ) + + b"\n" + for i, ( umi, - b','.join([findex[f] for f in feats]), - str(count).encode(), + feats, + count, pumi, - b','.join([findex[f] for f in pfeats])] - ) + b'\n' - for i, (umi, feats, count, pumi, pfeats) in dump.items() + pfeats, + ) in dump.items() ] else: dumplines = [ - b'\t'.join( - [str(cellidx).encode(), - cellbarcode, - str(i).encode(), - readname, - b','.join([findex[f] for f in feats]), - str(count).encode()] - ) + b'\n' + b"\t".join( + [ + str(cellidx).encode(), + cellbarcode, + str(i).encode(), + readname, + b",".join([findex[f] for f in feats]), + str(count).encode(), + ] + ) + + b"\n" for i, (readname, feats, count) in dump.items() ] df.writelines(dumplines) return matrix_file, dump_file + def formatMM(matrix_files, feature_index, barcodes_chunks, outdir): if type(matrix_files) is str: matrix_files = [matrix_files] - matrix_out = os.path.join(outdir, 'matrix.mtx.gz') + matrix_out = os.path.join(outdir, "matrix.mtx.gz") features_count = len(feature_index) barcodes_count = sum(len(x) for _, x in barcodes_chunks) mmsize = sum(getlen(f) for f in matrix_files) - mmheader = b'%%MatrixMarket matrix coordinate real general\n' - mmtotal = f'{features_count} {barcodes_count} {mmsize}\n'.encode() - with gzip.GzipFile(matrix_out, 'wb', mtime=0) as mmout: + mmheader = b"%%MatrixMarket matrix coordinate real general\n" + mmtotal = f"{features_count} {barcodes_count} {mmsize}\n".encode() + with gzip.GzipFile(matrix_out, "wb", mtime=0) as mmout: mmout.write(mmheader) mmout.write(mmtotal) - mtxstr = ' '.join(matrix_files) - cmd = f'zcat {mtxstr} | LC_ALL=C sort -k2,2n -k1,1n | gzip >> {matrix_out}' + mtxstr = " ".join(matrix_files) + cmd = f"zcat {mtxstr} | LC_ALL=C sort -k2,2n -k1,1n | gzip >> {matrix_out}" run_shell_cmd(cmd) - return(matrix_out) + return matrix_out + def writeEC(ecdump_files, no_umi, outdir): if type(ecdump_files) is str: ecdump_files = [ecdump_files] - ecdump_out = os.path.join(outdir, 'ec_dump.tsv.gz') - ecdumpstr = ' '.join(ecdump_files) + ecdump_out = os.path.join(outdir, "ec_dump.tsv.gz") + ecdumpstr = " ".join(ecdump_files) if not no_umi: - header = '\t'.join([ - 'Barcode_id', - 'Barcode', - 'EqClass', - 'UMI', - 'Features', - 'Read_count', - 'Dedup_UMI', - 'Dedup_feature' - ]) + '\n' + header = ( + "\t".join( + [ + "Barcode_id", + "Barcode", + "EqClass", + "UMI", + "Features", + "Read_count", + "Dedup_UMI", + "Dedup_feature", + ] + ) + + "\n" + ) else: - header = '\t'.join([ - 'Barcode_id', - 'Barcode', - 'EqClass', - 'Read_name', - 'Features', - 'Read_count' - ]) + '\n' - with gzip.GzipFile(ecdump_out, 'wb', mtime=0) as f: + header = ( + "\t".join( + [ + "Barcode_id", + "Barcode", + "EqClass", + "Read_name", + "Features", + "Read_count", + ] + ) + + "\n" + ) + with gzip.GzipFile(ecdump_out, "wb", mtime=0) as f: f.write(header.encode()) - cmd = f'zcat {ecdumpstr} | LC_ALL=C sort -k1,1n -k3,3n | gzip >> {ecdump_out}' + cmd = f"zcat {ecdumpstr} | LC_ALL=C sort -k1,1n -k3,3n | gzip >> {ecdump_out}" run_shell_cmd(cmd) - return ecdump_out \ No newline at end of file + return ecdump_out diff --git a/irescue/em.py b/irescue/em.py index 5625b44..32c94e7 100644 --- a/irescue/em.py +++ b/irescue/em.py @@ -1,31 +1,35 @@ import numpy as np + def e_step(matrix, counts): """ Performs E-step of EM algorithm: proportionally assigns reads to features based on relative feature abundances. """ - colsums = (matrix * counts).sum(axis=1)[:, np.newaxis] - out = matrix / colsums * counts - return(out) + colsums = matrix.dot(counts)[:, np.newaxis] + out = (matrix / colsums).multiply(counts) + return out + def m_step(matrix): """ Performs M-step of EM algorithm: calculates feature abundances from read counts proportionally distributed to features. """ - counts = matrix.sum(axis=0) / matrix.sum() - return(counts) + counts = np.ravel(matrix.sum(axis=0) / matrix.sum()) + return counts + def log_likelihood(matrix, counts): """ Compute log-likelihood of data. """ - likelihoods = (matrix * counts).sum(axis=1) + likelihoods = matrix.dot(counts) log_likelihood = np.sum(np.log(likelihoods + np.finfo(float).eps)) return log_likelihood -def run_em(matrix, cycles=100, tolerance=1e-5): + +def run_em(matrix, cycles=100, tolerance=1e-4): """ Run Expectation-Maximization (EM) algorithm to redistribute read counts across a set of features. @@ -53,7 +57,7 @@ def run_em(matrix, cycles=100, tolerance=1e-5): # (let the sum of counts of features be 1, # will be multiplied by the real UMI count later) nFeatures = matrix.shape[1] - counts = np.array([1 / nFeatures] * nFeatures) + counts = np.full(shape=nFeatures, fill_value=1 / nFeatures) # Initial log-likelihood prev_loglik = log_likelihood(matrix, counts) diff --git a/irescue/main.py b/irescue/main.py index d398555..9df89be 100644 --- a/irescue/main.py +++ b/irescue/main.py @@ -3,97 +3,234 @@ import argparse import os import sys -from multiprocessing import Pool from functools import partial +from multiprocessing import Pool from shutil import rmtree -from irescue._version import __version__ + from irescue._genomes import __genomes__ -from irescue.misc import writerr, versiontuple, run_shell_cmd -from irescue.misc import check_requirement, check_tags -from irescue.map import makeRmsk, getRefs, prepare_whitelist, isec, chrcat -from irescue.map import checkIndex -from irescue.count import split_barcodes, index_features, run_count, formatMM, writeEC +from irescue._version import __version__ +from irescue.count import ( + formatMM, + index_features, + run_count, + split_barcodes, + writeEC, +) +from irescue.map import ( + checkIndex, + chrcat, + getRefs, + isec, + makeRmsk, + prepare_whitelist, +) +from irescue.misc import ( + check_requirement, + check_tags, + run_shell_cmd, + versiontuple, + writerr, +) + def parseArguments(): parser = argparse.ArgumentParser( - prog='IRescue', + prog="IRescue", usage="irescue -b " - " [-g GENOME_ASSEMBLY | -r BED_FILE] [OPTIONS]", - description="IRescue (Interspersed Repeats single-cell quantifier):" - " a tool for quantifying transposable elements expression" - " in scRNA-seq.", - epilog="Home page: https://github.com/bodegalab/irescue" - ) - parser.add_argument('-b', '--bam', required=True, metavar='FILE', - help="scRNA-seq reads aligned to a reference genome (required).") - parser.add_argument('-r', '--regions', metavar='FILE', - help="Genomic TE coordinates in bed format (at least 4 columns with TE feature name " - "(e.g. subfamily) as the 4th column). Takes priority over --genome (default: %(default)s).") - parser.add_argument('-g', '--genome', metavar='STR', choices=__genomes__.keys(), - help="Genome assembly symbol. One of: {} (default: " - "%(default)s).".format(', '.join(__genomes__))) - parser.add_argument('-w', '--whitelist', metavar='FILE', - help="Text file of filtered cell barcodes by e.g. Cell Ranger, STARSolo " - "or your gene expression quantifier of choice (Recommended. default: %(default)s).") - parser.add_argument('-c', '--cb-tag', default='CB', metavar='STR', - help="BAM tag containing the cell barcode sequence (default: %(default)s).") - parser.add_argument('-u', '--umi-tag', default='UR', metavar='STR', - help="BAM tag containing the UMI sequence (default: %(default)s).") - parser.add_argument('--no-umi', action='store_true', - help="Ignore UMI sequence (for UMI-less technologies, such as SMART-seq).") - parser.add_argument('-p', '--threads', type=int, default=1, metavar='CPUS', - help="Number of cpus to use (default: %(default)s).") - parser.add_argument('-o', '--outdir', default='irescue_out', metavar='DIR', - help="Output directory name (default: %(default)s).") - parser.add_argument('--min-bp-overlap', type=int, metavar='INT', - help="Minimum overlap between read and TE as number of nucleotides (Default: disabled).") - parser.add_argument('--min-fraction-overlap', type=float, metavar='FLOAT', choices=[x/100 for x in range(101)], - help="Minimum overlap between read and TE as a fraction of read's alignment" - " (i.e. 0.00 <= NUM <= 1.00) (Default: disabled).") - parser.add_argument('--max-iters', type=int, metavar='INT', default=100, - help="Maximum number of EM iterations (Default: %(default)s).") - parser.add_argument('--tolerance', type=float, metavar='FLOAT', default=1e-4, - help="Log-likelihood change below which convergence is assumed (Default: %(default)s).") - parser.add_argument('--dump-ec', action='store_true', - help="Write a description log file of Equivalence Classes.") - parser.add_argument('--integers', action='store_true', - help="Use if integers count are needed for downstream analysis.") - parser.add_argument('--samtools', default='samtools', metavar='PATH', - help="Path to samtools binary, in case it's not in PATH (Default: %(default)s).") - parser.add_argument('--bedtools', default='bedtools', metavar='PATH', - help="Path to bedtools binary, in case it's not in PATH (Default: %(default)s).") - parser.add_argument('--no-tags-check', action='store_true', - help="Suppress checking for CBtag and UMItag presence in BAM file.") - parser.add_argument('--keeptmp', action='store_true', - help="Keep temporary files under /tmp.") - parser.add_argument('-v', '--verbose', action='count', default=0, - help="Writes additional logging to stderr. Use once for normal verbosity (-v), " - "twice for debugging (-vv).") - parser.add_argument('-V', '--version', action='version', version='%(prog)s {}'.format(__version__), - help="Print software's version and exit.") + " [-g GENOME_ASSEMBLY | -r BED_FILE] [OPTIONS]", + description=( + "IRescue (Interspersed Repeats single-cell quantifier): a tool" + "for quantifying transposable elements expression in scRNA-seq." + ), + epilog="Home page: https://github.com/bodegalab/irescue", + ) + parser.add_argument( + "-b", + "--bam", + required=True, + metavar="FILE", + help="scRNA-seq reads aligned to a reference genome (required).", + ) + parser.add_argument( + "-r", + "--regions", + metavar="FILE", + help=( + "Genomic TE coordinates in bed format (at least 4 columns with TE" + " feature name (e.g. subfamily) as the 4th column). " + "Takes priority over --genome (default: %(default)s)." + ), + ) + parser.add_argument( + "-g", + "--genome", + metavar="STR", + choices=__genomes__.keys(), + help=( + "Genome assembly symbol. One of: {} (default: " + "%(default)s).".format(", ".join(__genomes__)) + ), + ) + parser.add_argument( + "-w", + "--whitelist", + metavar="FILE", + help=( + "Text file of filtered cell barcodes by e.g. Cell Ranger, " + "STARSolo or your gene expression quantifier of choice " + "(Recommended. default: %(default)s)." + " Note: If not provided, all barcodes found in BAM will be used." + ), + ) + parser.add_argument( + "-c", + "--cb-tag", + default="CB", + metavar="STR", + help=( + "BAM tag containing the cell barcode sequence" + " (default: %(default)s)." + ), + ) + parser.add_argument( + "-u", + "--umi-tag", + default="UR", + metavar="STR", + help="BAM tag containing the UMI sequence (default: %(default)s).", + ) + parser.add_argument( + "--no-umi", + action="store_true", + help="Ignore UMI sequence (for UMI-less datasets, such as Smart-seq).", + ) + parser.add_argument( + "-p", + "--threads", + type=int, + default=1, + metavar="CPUS", + help="Number of cpus to use (default: %(default)s).", + ) + parser.add_argument( + "-o", + "--outdir", + default="irescue_out", + metavar="DIR", + help="Output directory name (default: %(default)s).", + ) + parser.add_argument( + "--min-bp-overlap", + type=int, + metavar="INT", + help=( + "Minimum overlap between read and TE as number of nucleotides " + "(Default: disabled)." + ), + ) + parser.add_argument( + "--min-fraction-overlap", + type=float, + metavar="FLOAT", + choices=[x / 100 for x in range(101)], + help=( + "Minimum overlap between read and TE as a fraction of " + "read's alignment (i.e. 0.00 <= NUM <= 1.00) (Default: disabled)." + ), + ) + parser.add_argument( + "--max-iters", + type=int, + metavar="INT", + default=100, + help="Maximum number of EM iterations (Default: %(default)s).", + ) + parser.add_argument( + "--tolerance", + type=float, + metavar="FLOAT", + default=1e-4, + help=( + "Log-likelihood change below which convergence is assumed " + "(Default: %(default)s)." + ), + ) + parser.add_argument( + "--dump-ec", + action="store_true", + help="Write a description log file of Equivalence Classes.", + ) + parser.add_argument( + "--integers", + action="store_true", + help="Use if integers count are needed for downstream analysis.", + ) + parser.add_argument( + "--samtools", + default="samtools", + metavar="PATH", + help=( + "Path to samtools binary, in case it's not in PATH " + "(Default: %(default)s)." + ), + ) + parser.add_argument( + "--bedtools", + default="bedtools", + metavar="PATH", + help=( + "Path to bedtools binary, in case it's not in PATH " + "(Default: %(default)s)." + ), + ) + parser.add_argument( + "--no-tags-check", + action="store_true", + help="Suppress checking for CBtag and UMItag presence in BAM file.", + ) + parser.add_argument( + "--keeptmp", + action="store_true", + help="Keep temporary files under /tmp.", + ) + parser.add_argument( + "-v", + "--verbose", + action="count", + default=0, + help=( + "Writes additional logging to stderr. " + "Use once for normal verbosity (-v), twice for debugging (-vv)." + ), + ) + parser.add_argument( + "-V", + "--version", + action="version", + version="%(prog)s {}".format(__version__), + help="Print software's version and exit.", + ) return parser def main(): - # Parse and print arguments parser = parseArguments() - args = parser.parse_args(args=None if sys.argv[1:] else ['--help']) - + args = parser.parse_args(args=None if sys.argv[1:] else ["--help"]) + if args.no_umi: args.umi_tag = "" - argstr = '\n'.join(f' {k}: {v}' for k, v in args.__dict__.items()) + argstr = "\n".join(f" {k}: {v}" for k, v in args.__dict__.items()) sys.stderr.write(f" IRescue version {__version__}\n{argstr}\n") - #__tmpdir__ = os.path.join(args.outdir, 'tmp') + # __tmpdir__ = os.path.join(args.outdir, 'tmp') dirs = { - 'out': args.outdir, - 'tmp': os.path.join(args.outdir, 'tmp'), - 'mex': os.path.join(args.outdir, 'counts') + "out": args.outdir, + "tmp": os.path.join(args.outdir, "tmp"), + "mex": os.path.join(args.outdir, "counts"), } - #################### # Preliminar steps # #################### @@ -102,22 +239,30 @@ def main(): # Check requirements check_requirement( - args.bedtools, '2.30.0', + args.bedtools, + "2.30.0", lambda: versiontuple( - run_shell_cmd('bedtools --version').split()[1][1:] + run_shell_cmd("bedtools --version").split()[1][1:] ), - args.verbose + args.verbose, ) check_requirement( - args.samtools, '1.11', - lambda: versiontuple(run_shell_cmd('samtools --version').split()[1]), - args.verbose + args.samtools, + "1.11", + lambda: versiontuple(run_shell_cmd("samtools --version").split()[1]), + args.verbose, ) # Check if the selected cell barcode and UMI tags are present in bam file. if not args.no_tags_check: - check_tags(bamFile=args.bam, CBtag=args.cb_tag, UMItag=args.umi_tag, - nLines=999999, exit_with_error=True, verbose=args.verbose) + check_tags( + bamFile=args.bam, + CBtag=args.cb_tag, + UMItag=args.umi_tag, + nLines=999999, + exit_with_error=True, + verbose=args.verbose, + ) # Check for bam index file. If not present, will build an index. checkIndex(args.bam, verbose=args.verbose) @@ -126,7 +271,6 @@ def main(): for v in dirs.values(): os.makedirs(v, exist_ok=True) - ########### # Mapping # ########### @@ -134,15 +278,19 @@ def main(): writerr("Running mapping step.") # set regions object (provided or downloaded bed file) - regions = makeRmsk(regions=args.regions, genome=args.genome, - genomes=__genomes__, tmpdir=dirs['tmp'], - outname='rmsk.bed') + regions = makeRmsk( + regions=args.regions, + genome=args.genome, + genomes=__genomes__, + tmpdir=dirs["tmp"], + outname="rmsk.bed", + ) # get list of reference names from bam chrNames = getRefs(args.bam, regions) # decompress whitelist if compressed - whitelist = prepare_whitelist(args.whitelist, dirs['tmp']) + whitelist = prepare_whitelist(args.whitelist, dirs["tmp"]) # Allocate threads if args.threads > 1: @@ -151,13 +299,23 @@ def main(): # Execute intersection between reads and TE coordinates writerr( "Computing overlap between reads and TEs coordinates in the " - "following references: {}".format(', '.join(chrNames)), - level=1, send=args.verbose + "following references: {}".format(", ".join(chrNames)), + level=1, + send=args.verbose, ) isecFun = partial( - isec, args.bam, regions, whitelist, args.cb_tag, args.umi_tag, - args.min_bp_overlap, args.min_fraction_overlap, dirs['tmp'], - args.samtools, args.bedtools, args.verbose + isec, + args.bam, + regions, + whitelist, + args.cb_tag, + args.umi_tag, + args.min_bp_overlap, + args.min_fraction_overlap, + dirs["tmp"], + args.samtools, + args.bedtools, + args.verbose, ) if args.threads > 1: isecFiles = pool.map(isecFun, chrNames) @@ -166,11 +324,14 @@ def main(): # concatenate intersection results mappings_file, barcodes_file, features_file = chrcat( - isecFiles, threads=args.threads, outdir=dirs['mex'], - tmpdir=dirs['tmp'], bedtools=args.bedtools, verbose=args.verbose + isecFiles, + threads=args.threads, + outdir=dirs["mex"], + tmpdir=dirs["tmp"], + bedtools=args.bedtools, + verbose=args.verbose, ) - ######### # Count # ######### @@ -184,8 +345,15 @@ def main(): # calculate TE counts countFun = partial( - run_count, mappings_file, feature_index, dirs['tmp'], args.no_umi, - args.dump_ec, args.max_iters, args.tolerance, args.verbose + run_count, + mappings_file, + feature_index, + dirs["tmp"], + args.no_umi, + args.dump_ec, + args.max_iters, + args.tolerance, + args.verbose, ) if args.threads > 1: mtxFiles = pool.map(countFun, bc_per_thread) @@ -198,18 +366,18 @@ def main(): pool.join() # concatenate matrix files chunks - matrix_files = [ i for i, j in mtxFiles] - ecdump_files = [ j for i, j in mtxFiles] + matrix_files = [i for i, j in mtxFiles] + ecdump_files = [j for i, j in mtxFiles] matrix_file = formatMM( - matrix_files, feature_index, bc_per_thread, dirs['mex'] + matrix_files, feature_index, bc_per_thread, dirs["mex"] ) - writerr(f'Writing sparse matrix to {matrix_file}') + writerr(f"Writing sparse matrix to {matrix_file}") if args.dump_ec: - ecdump_file = writeEC(ecdump_files, args.no_umi, outdir=dirs['out']) - writerr(f'Writing Equivalence Classes to {ecdump_file}') + ecdump_file = writeEC(ecdump_files, args.no_umi, outdir=dirs["out"]) + writerr(f"Writing Equivalence Classes to {ecdump_file}") if not args.keeptmp: - writerr('Cleaning up temporary files.', level=1, send=args.verbose) - rmtree(dirs['tmp']) + writerr("Cleaning up temporary files.", level=1, send=args.verbose) + rmtree(dirs["tmp"]) - writerr('Done.') + writerr("Done.") diff --git a/irescue/map.py b/irescue/map.py index e3509b4..bebc388 100644 --- a/irescue/map.py +++ b/irescue/map.py @@ -1,35 +1,38 @@ #!/usr/bin/env python -import requests import io import os -from pysam import idxstats, AlignmentFile, index from gzip import open as gzopen -from irescue.misc import testGz -from irescue.misc import writerr -from irescue.misc import unGzip -from irescue.misc import run_shell_cmd -from irescue.misc import getlen + +import requests +from pysam import AlignmentFile, idxstats, index + +from irescue.misc import getlen, run_shell_cmd, testGz, unGzip, writerr + # Check if bam file is indexed def checkIndex(bamFile, verbose): with AlignmentFile(bamFile) as bam: if not bam.has_index(): - writerr('BAM index not found. Attempting to index the BAM...') + writerr("BAM index not found. Attempting to index the BAM...") try: index(bamFile) except Exception as e: writerr( "ERROR: Couldn't index the BAM file. Is your BAM file " f"sorted? If not, please sort it by coordinate.\n\n{e}", - error=True + error=True, ) else: - writerr('BAM indexing done.') + writerr("BAM indexing done.") else: if verbose: - writerr(f'Found index for BAM file {bamFile}.', - level=1, send=verbose) + writerr( + f"Found index for BAM file {bamFile}.", + level=1, + send=verbose, + ) + # Check repeatmasker regions bed file format. Download if not provided. # Returns the path of the repeatmasker bed file. @@ -37,24 +40,27 @@ def makeRmsk(regions, genome, genomes, tmpdir, outname): # if a repeatmasker bed file is provided, use that if regions: if testGz(regions): - f = gzopen(regions, 'rb') + f = gzopen(regions, "rb") + def rl(x): return x.readline().decode() else: - f = open(regions, 'r') + f = open(regions, "r") + def rl(x): return x.readline() + # skip header line = rl(f) - while line[0] == '#': + while line[0] == "#": line = rl(f) # check for minimum column number - if len(line.strip().split('\t')) < 4: + if len(line.strip().split("\t")) < 4: writerr( "Error: please provide a tab-separated BED file with at " "least 4 columns and TE feature name (e.g. subfamily) " "in 4th column.", - error=True + error=True, ) f.close() out = regions @@ -64,89 +70,93 @@ def rl(x): url, header_lines = genomes[genome] writerr( "Downloading and parsing RepeatMasker annotation for " - f"assembly {genome} from {url} ...") + f"assembly {genome} from {url} ..." + ) try: response = requests.get(url, stream=True, timeout=60) except Exception as e: writerr( "ERROR: Download of RepeatMasker annotation failed. " f"Couldn't connect to host.\n\n{e}", - error=True + error=True, ) - rmsk = gzopen(io.BytesIO(response.content), 'rb') + rmsk = gzopen(io.BytesIO(response.content), "rb") out = os.path.join(tmpdir, outname) - with open(out, 'w') as f: + with open(out, "w") as f: # print header - h = ['#chr','start','end','name','score','strand'] - h = '\t'.join(h) - h += '\n' + h = ["#chr", "start", "end", "name", "score", "strand"] + h = "\t".join(h) + h += "\n" f.write(h) # skip rmsk header for _ in range(header_lines): next(rmsk) # parse rmsk fams_to_skip = [ - 'Low_complexity', - 'Simple_repeat', - 'rRNA', - 'scRNA', - 'srpRNA', - 'tRNA' + "Low_complexity", + "Simple_repeat", + "rRNA", + "scRNA", + "srpRNA", + "tRNA", ] for line in rmsk: - lst = line.decode('utf-8').strip().split() + lst = line.decode("utf-8").strip().split() strand, subfamily, famclass = lst[8:11] - if famclass.split('/')[0] in fams_to_skip: + if famclass.split("/")[0] in fams_to_skip: continue # concatenate family and class with subfamily - subfamily += '#' + famclass + subfamily += "#" + famclass score = lst[0] chr, start, end = lst[4:7] # make coordinates 0-based - start = str(int(start)-1) - if strand != '+': - strand = '-' - outl = '\t'.join([chr, start, end, subfamily, score, strand]) - outl += '\n' + start = str(int(start) - 1) + if strand != "+": + strand = "-" + outl = "\t".join([chr, start, end, subfamily, score, strand]) + outl += "\n" f.write(outl) else: writerr( "Error: it is mandatory to define either --regions OR " "--genome parameter.", - error=True + error=True, ) - return(out) + return out + # Uncompress the whitelist file if compressed. # Return the whitelist path, or False if not using a whitelist. def prepare_whitelist(whitelist, tmpdir): if whitelist and testGz(whitelist): - wlout = os.path.join(tmpdir, 'whitelist.tsv') + wlout = os.path.join(tmpdir, "whitelist.tsv") whitelist = unGzip(whitelist, wlout) return whitelist + # Get list of reference names from BAM file, skipping those without reads. def getRefs(bamFile, bedFile): chrNames = list() - for line in idxstats(bamFile).strip().split('\n'): - fields = line.strip().split('\t') - if int(fields[2])>0: + for line in idxstats(bamFile).strip().split("\n"): + fields = line.strip().split("\t") + if int(fields[2]) > 0: chrNames.append(fields[0]) bedChrNames = set() if testGz(bedFile): - with gzopen(bedFile, 'rb') as f: + with gzopen(bedFile, "rb") as f: for line in f: - bedChrNames.add(line.decode().split('\t')[0]) + bedChrNames.add(line.decode().split("\t")[0]) else: - with open(bedFile, 'r') as f: + with open(bedFile, "r") as f: for line in f: - bedChrNames.add(line.split('\t')[0]) + bedChrNames.add(line.split("\t")[0]) skipChr = [x for x in chrNames if x not in bedChrNames] if skipChr: writerr( "WARNING: The following references contain read alignments but " "are not found in the TE annotation and will be skipped: " - f"{', '.join(skipChr)}") + f"{', '.join(skipChr)}" + ) chrNames = [x for x in chrNames if x in bedChrNames] if chrNames: return chrNames @@ -158,39 +168,56 @@ def getRefs(bamFile, bedFile): you can either change it to UCSC (chr1, chr2, etc...), or use a custom TE annotation with ENSEMBL chromosome names. """, - error=True + error=True, ) + # Intersect reads with repeatmasker regions. Return the intersection file path. -def isec(bamFile, bedFile, whitelist, CBtag, UMItag, bpOverlap, fracOverlap, - tmpdir, samtools, bedtools, verbose, chrom): - refdir = os.path.join(tmpdir, 'refs') - isecdir = os.path.join(tmpdir, 'isec') +def isec( + bamFile, + bedFile, + whitelist, + CBtag, + UMItag, + bpOverlap, + fracOverlap, + tmpdir, + samtools, + bedtools, + verbose, + chrom, +): + refdir = os.path.join(tmpdir, "refs") + isecdir = os.path.join(tmpdir, "isec") os.makedirs(refdir, exist_ok=True) os.makedirs(isecdir, exist_ok=True) - refFile = os.path.join(refdir, chrom + '.bed.gz') - isecFile = os.path.join(isecdir, chrom + '.isec.txt.gz') + refFile = os.path.join(refdir, chrom + ".bed.gz") + isecFile = os.path.join(isecdir, chrom + ".isec.txt.gz") # split bed file by chromosome - sort = 'LC_ALL=C sort -k1,1 -k2,2n --buffer-size=1G' - if bedFile[-3:] == '.gz': - cmd0 = f'zcat {bedFile} | gawk \'$1=="{chrom}"\' ' - cmd0 += f' | {sort} | gzip > {refFile}' + sort = "LC_ALL=C sort -k1,1 -k2,2n --buffer-size=1G" + if bedFile[-3:] == ".gz": + cmd0 = f"zcat {bedFile} | gawk '$1==\"{chrom}\"' " + cmd0 += f" | {sort} | gzip > {refFile}" else: - cmd0 = f'gawk \'$1=="{chrom}"\' {bedFile} | {sort} | gzip > {refFile}' + cmd0 = f"gawk '$1==\"{chrom}\"' {bedFile} | {sort} | gzip > {refFile}" # command streaming alignments for intersection if whitelist: - stream = f' <({samtools} view -h {bamFile} -D {CBtag}:{whitelist} {chrom} | ' + stream = f" <({samtools} view -h {bamFile} -D {CBtag}:{whitelist} {chrom} | " else: - stream = f' <({samtools} view -h {bamFile} {chrom} | ' + stream = f" <({samtools} view -h {bamFile} {chrom} | " stream += ' gawk \'!($1~/^@/) { split("", tags); ' - stream += ' for (i=12;i<=NF;i++) {split($i,tag,":"); tags[tag[1]]=tag[3]}; ' + stream += ( + ' for (i=12;i<=NF;i++) {split($i,tag,":"); tags[tag[1]]=tag[3]}; ' + ) # Discard records without CB tag, unvalid STARSolo CBs, missing UMI tag, # UMIs with Ns and homopolymer UMIs if UMItag: - stream += f' if(tags["{CBtag}"]~/^(|-)$/ || tags["{UMItag}"]~/.*N.*/ || ' + stream += ( + f' if(tags["{CBtag}"]~/^(|-)$/ || tags["{UMItag}"]~/.*N.*/ || ' + ) stream += f' tags["{UMItag}"]~/^$|^(A+|G+|T+|C+)$/) {{next}}; ' # Append CB and UMI to read name stream += f' $1=$1"/"tags["{CBtag}"]"/"tags["{UMItag}"]; ' @@ -198,89 +225,90 @@ def isec(bamFile, bedFile, whitelist, CBtag, UMItag, bpOverlap, fracOverlap, stream += f' if(tags["{CBtag}"]~/^(|-)$/) {{next}}; ' # Append CB to read name, and read name again (replacing UMI) stream += f' $1=$1"/"tags["{CBtag}"]"/"$1; ' - stream += ' } ' + stream += " } " stream += ' { OFS="\\t"; print }\' | ' - stream += f' {samtools} view -u - | ' - stream += f' {bedtools} bamtobed -i stdin -bed12 -split -splitD) ' + stream += f" {samtools} view -u - | " + stream += f" {bedtools} bamtobed -i stdin -bed12 -split -splitD) " # filter by minimum overlap between read and feature, if set - ovfrac = f' -f {fracOverlap} ' if fracOverlap else '' - ovbp = f' $NF>={bpOverlap} ' if bpOverlap else '' + ovfrac = f" -f {fracOverlap} " if fracOverlap else "" + ovbp = f" $NF>={bpOverlap} " if bpOverlap else "" # intersection command - cmd = f'{bedtools} intersect -a {stream} -b {refFile} ' + cmd = f"{bedtools} intersect -a {stream} -b {refFile} " cmd += f' -split -bed -wo -sorted {ovfrac} | gawk -vOFS="\\t" \'{ovbp} ' # remove mate information from read name cmd += ' { sub(/\\/[12]$/,"",$4); ' # concatenate CB and UMI with feature name - cmd += ' n=split($4,qname,/\\//); ' + cmd += " n=split($4,qname,/\\//); " cmd += ' print qname[n-1]"\\t"qname[n]"\\t"qname[1]"\\t"$16 }\' ' - cmd += f' | gzip > {isecFile}' + cmd += f" | gzip > {isecFile}" - writerr(f'Extracting {chrom} reference', level=2, send=verbose) + writerr(f"Extracting {chrom} reference", level=2, send=verbose) run_shell_cmd(cmd0) - writerr(f'Mapping alignments to {chrom}', level=1, send=verbose) + writerr(f"Mapping alignments to {chrom}", level=1, send=verbose) run_shell_cmd(cmd) - writerr(f'Mapped {chrom}', level=1, send=verbose) + writerr(f"Mapped {chrom}", level=1, send=verbose) return isecFile + # Concatenate and sort data obtained from isec() def chrcat(filesList, threads, outdir, tmpdir, bedtools, verbose): os.makedirs(outdir, exist_ok=True) - mappings_file = os.path.join(tmpdir, 'mappings.tsv.gz') - barcodes_file = os.path.join(outdir, 'barcodes.tsv.gz') - features_file = os.path.join(outdir, 'features.tsv.gz') - bedFiles = ' '.join(filesList) + mappings_file = os.path.join(tmpdir, "mappings.tsv.gz") + barcodes_file = os.path.join(outdir, "barcodes.tsv.gz") + features_file = os.path.join(outdir, "features.tsv.gz") + bedFiles = " ".join(filesList) sort_threads = int(threads / 2 - 1) - sort_threads = sort_threads if sort_threads>0 else 1 + sort_threads = sort_threads if sort_threads > 0 else 1 # sort and summarize UMI-READ-TE mappings # (if --no-umi, UMI is replaced by READ) - sort_res = f'--parallel {sort_threads} --buffer-size 2G' - cmd0 = f'zcat {bedFiles}' - # input: "CB UMI READ FEAT" - cmd0 += f' | LC_ALL=C sort -u {sort_res}' - cmd0 += f' | {bedtools} groupby -g 1,2,3 -c 4 -o distinct' - # result: "CB UMI READ FEATs" - cmd0 += f' | LC_ALL=C sort -k1,2 -k4,4 {sort_res}' - cmd0 += f' | {bedtools} groupby -g 1,2,4 -c 3 -o count_distinct' - # result: "CB UMI FEATs count" - cmd0 += f' | gzip > {mappings_file}' + sort_res = f"--parallel {sort_threads} --buffer-size 2G" + cmd0 = f"zcat {bedFiles}" + # input: "CB UMI READ FEAT" + cmd0 += f" | LC_ALL=C sort -u {sort_res}" + cmd0 += f" | {bedtools} groupby -g 1,2,3 -c 4 -o distinct" + # result: "CB UMI READ FEATs" + cmd0 += f" | LC_ALL=C sort -k1,2 -k4,4 {sort_res}" + cmd0 += f" | {bedtools} groupby -g 1,2,4 -c 3 -o count_distinct" + # result: "CB UMI FEATs count" + cmd0 += f" | gzip > {mappings_file}" # write barcodes.tsv.gz file - cmd1 = f'zcat {mappings_file} | cut -f1 | uniq | gzip > {barcodes_file} ' + cmd1 = f"zcat {mappings_file} | cut -f1 | uniq | gzip > {barcodes_file} " # write features.tsv.gz file - cmd2 = f'zcat {mappings_file} ' - cmd2 += ' | cut -f3 | sed \'s/,/\\n/g\' | gawk \'!x[$1]++ { ' + cmd2 = f"zcat {mappings_file} " + cmd2 += " | cut -f3 | sed 's/,/\\n/g' | gawk '!x[$1]++ { " cmd2 += ' print $1"\\t"gensub(/#.+/,"",1,$1)"\\tGene Expression" }\' ' - cmd2 += f' | LC_ALL=C sort -u | gzip > {features_file} ' + cmd2 += f" | LC_ALL=C sort -u | gzip > {features_file} " - writerr('Concatenating mappings', level=1, send=verbose) + writerr("Concatenating mappings", level=1, send=verbose) run_shell_cmd(cmd0) if getlen(mappings_file) == 0: writerr( - f'No read-TE mappings found in {mappings_file}.' - ' Check annotation and temporary files to troubleshoot.', - error=True + f"No read-TE mappings found in {mappings_file}." + " Check annotation and temporary files to troubleshoot.", + error=True, ) - writerr(f'Writing mapped barcodes to {barcodes_file}') + writerr(f"Writing mapped barcodes to {barcodes_file}") run_shell_cmd(cmd1) if getlen(barcodes_file) == 0: writerr( - f'No features written in {features_file}.' - ' Check BAM format and reference annotation (e.g. chr names)' - ' to troubleshoot.', - error=True + f"No features written in {features_file}." + " Check BAM format and reference annotation (e.g. chr names)" + " to troubleshoot.", + error=True, ) - writerr(f'Writing mapped features to {features_file}') + writerr(f"Writing mapped features to {features_file}") run_shell_cmd(cmd2) if getlen(features_file) == 0: writerr( - f'No features written in {features_file}.' - ' Check annotation and temporary files to troubleshoot.', - error=True + f"No features written in {features_file}." + " Check annotation and temporary files to troubleshoot.", + error=True, ) return mappings_file, barcodes_file, features_file diff --git a/irescue/misc.py b/irescue/misc.py index 44556b7..5c94e62 100644 --- a/irescue/misc.py +++ b/irescue/misc.py @@ -1,29 +1,32 @@ #!/usr/bin/env python -import subprocess +import gzip import os +import subprocess import sys -import gzip from datetime import datetime from shutil import which + import pysam + def run_shell_cmd(cmd): """ Execute a command on bash shell with subprocess. """ p = subprocess.Popen( - ['/bin/bash', '-o', 'pipefail'], + ["/bin/bash", "-o", "pipefail"], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE, preexec_fn=os.setsid, - text=True + text=True, ) pid = p.pid os.getpgid(pid) stdout, stdin = p.communicate(cmd) - return stdout.strip('\n') + return stdout.strip("\n") + def check_path(cmdname): """ @@ -31,11 +34,13 @@ def check_path(cmdname): """ return which(cmdname) is not None + def versiontuple(version): """ Convert a semver string "X.Y.Z" to a tuple of integers (X, Y, Z). """ - return tuple(map(int, version.split('.'))) + return tuple(map(int, version.split("."))) + def check_requirement(cmd, required_version, parser, verbose): """ @@ -45,7 +50,7 @@ def check_requirement(cmd, required_version, parser, verbose): writerr( f"ERROR: Couldn't find {cmd} in PATH. Please install " f"{cmd} >={required_version} and try again.", - error=True + error=True, ) else: try: @@ -58,7 +63,8 @@ def check_requirement(cmd, required_version, parser, verbose): else: writerr( f"Found {cmd} version {version}. Proceeding.", - level=1, send=verbose + level=1, + send=verbose, ) except Exception as e: writerr( @@ -67,6 +73,7 @@ def check_requirement(cmd, required_version, parser, verbose): f"not supported.\n\n{e}" ) + # Small function to write a message to stderr with timestamp def writerr(msg, error=False, level=0, send=0): """ @@ -83,28 +90,30 @@ def writerr(msg, error=False, level=0, send=0): send: int If >=verbosity, the message will be sent. """ - if send>=level or error: + if send >= level or error: timelog = datetime.now().strftime("%Y/%m/%d - %H:%M:%S") - message = f'[{timelog}] ' - if not msg[-1]=='\n': - msg += '\n' + message = f"[{timelog}] " + if not msg[-1] == "\n": + msg += "\n" message += msg if error: sys.exit(message) else: sys.stderr.write(message) + def testGz(input_file): """ Check if a file is gzip compressed. """ - with gzip.open(input_file, 'rb') as f: + with gzip.open(input_file, "rb") as f: try: f.read(1) return True except gzip.BadGzipFile: return False + # Uncompress gzipped file def unGzip(input_file, output_file): """ @@ -115,27 +124,27 @@ def unGzip(input_file, output_file): input_file: gzip file to decompress. output_file: uncompressed file to write. """ - with gzip.open(input_file, 'rb') as fin,\ - open(output_file, 'w') as fout: + with gzip.open(input_file, "rb") as fin, open(output_file, "w") as fout: for line in fin: - fout.write(line.decode('utf-8')) + fout.write(line.decode("utf-8")) return output_file + def getlen(file): """ Count the number of lines in a file (plain or gzip-compressed). """ if testGz(file): - f = gzip.open(file, 'rb') + f = gzip.open(file, "rb") else: - f = open(file, 'r') + f = open(file, "r") out = sum(1 for line in f) f.close() return out + def check_tags( - bamFile, CBtag, UMItag, - nLines=None, exit_with_error=True, verbose=False + bamFile, CBtag, UMItag, nLines=None, exit_with_error=True, verbose=False ): """ Check if BAM file contains barcode and UMI tags. @@ -155,12 +164,13 @@ def check_tags( verbose: bool Write progress info to stderr. """ - with pysam.AlignmentFile(bamFile, 'rb') as f: + with pysam.AlignmentFile(bamFile, "rb") as f: c = 1 writerr( f"Testing bam file for {CBtag} {'and ' + UMItag if UMItag else ''}" "tags presence. Will stop at the first occurrence.", - level=1, send=verbose + level=1, + send=verbose, ) for read in f: if nLines and c >= nLines: @@ -181,10 +191,8 @@ def check_tags( ) else: read.get_tag(CBtag) - writerr( - f"Found {CBtag} tag occurrence in BAM's line {c}." - ) - return(True) + writerr(f"Found {CBtag} tag occurrence in BAM's line {c}.") + return True except Exception: c += 1 pass @@ -201,24 +209,25 @@ def check_tags( If you do not expect that all alignments contain the tags, you can suppress this check with --no-tags-check. """.format( - f"{CBtag} and/or {UMItag} tags" if UMItag else CBtag + ' tag', - f"the first {nLines} lines of " if nLines else '', - '--CBtag and --UMItag flags' if UMItag else '--CBtag flag' + f"{CBtag} and/or {UMItag} tags" if UMItag else CBtag + " tag", + f"the first {nLines} lines of " if nLines else "", + "--CBtag and --UMItag flags" if UMItag else "--CBtag flag", ), - error=True + error=True, ) else: - return(False) + return False + def get_ranges(num, div): """ Splits an integer X into N integers whose sum is equal to X. """ div = num if div > num else div - split = int(num/div) + split = int(num / div) for i in range(0, num, split): j = i + split - if j > num-split: + if j > num - split: j = num yield range(i, j) break diff --git a/irescue/network.py b/irescue/network.py index f274e3a..f5b4594 100644 --- a/irescue/network.py +++ b/irescue/network.py @@ -29,11 +29,12 @@ from collections import defaultdict + def get_substr_slices(umi_length, idx_size): - ''' + """ Create slices to split a UMI into approximately equal size substrings Returns a list of tuples that can be passed to slice function - ''' + """ cs, r = divmod(umi_length, idx_size) sub_sizes = [cs + 1] * r + [cs] * (idx_size - r) offset = 0 @@ -43,11 +44,12 @@ def get_substr_slices(umi_length, idx_size): offset += s return slices + def build_substr_idx(equivalence_classes, length, threshold): - ''' + """ Group equivalence classes into subgroups having a common substring - ''' - slices = get_substr_slices(length, threshold+1) + """ + slices = get_substr_slices(length, threshold + 1) substr_idx = {k: defaultdict(set) for k in slices} for idx in slices: for ec in equivalence_classes: @@ -55,10 +57,11 @@ def build_substr_idx(equivalence_classes, length, threshold): substr_idx[idx][sub].add(ec) return substr_idx + def gen_ec_pairs(equivalence_classes, substr_idx): - ''' + """ Yields equivalence classes pairs from build_substr_idx() - ''' + """ for i, ec in enumerate(equivalence_classes, start=1): neighbours = set() for idx, substr_map in substr_idx.items(): @@ -66,4 +69,4 @@ def gen_ec_pairs(equivalence_classes, substr_idx): neighbours = neighbours.union(substr_map[sub]) neighbours.difference_update(equivalence_classes[:i]) for nbr in neighbours: - yield ec, nbr \ No newline at end of file + yield ec, nbr diff --git a/pyproject.toml b/pyproject.toml index 61434eb..6527b1d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,6 +34,7 @@ dependencies = [ "pysam >= 0.16.0.1", "requests >= 2.27.1", "networkx >= 3.1", + "scipy >= 1.11.4" ] dynamic = ["version"]