diff --git a/irescue/main.py b/irescue/main.py index dc7d5e5..226197e 100644 --- a/irescue/main.py +++ b/irescue/main.py @@ -100,12 +100,15 @@ def parseArguments(): help="BAM tag containing the UMI sequence (default: %(default)s).", ) parser.add_argument( - "-l", - "--locus", - action="store_true", + "--locus-level", + type=str, + metavar="STR", + choices=["disabled", "fragment", "instance"], + default="disabled", help=( - "Perform locus-level quantification, instead of subfamily-level" - " (default: %(default)s)." + "Perform locus-level quantification. " + "One of: disabled, fragment, instance (default: %(default)s). " + "Disabled means subfamily-level quantification." ), ) parser.add_argument( @@ -306,7 +309,7 @@ def main(): genome=args.genome, genomes=__genomes__, outdir=dirs["out"], - locus=args.locus, + locus=args.locus_level, outname="rmsk.bed.gz", ) @@ -353,7 +356,7 @@ def main(): threads=args.threads, outdir=dirs["mex"], tmpdir=dirs["tmp"], - locus=args.locus, + locus=args.locus_level, bedtools=args.bedtools, verbose=args.verbose, ) diff --git a/irescue/map.py b/irescue/map.py index 06ddecf..bda5a93 100644 --- a/irescue/map.py +++ b/irescue/map.py @@ -35,9 +35,7 @@ def checkIndex(bamFile, verbose): ) -def makeRmsk( - regions, genome, genomes, outdir, locus=False, outname="rmsk.bed.gz" -): +def makeRmsk(regions, genome, genomes, outdir, locus="disabled", outname="rmsk.bed.gz"): """Format and/or download RepeatMasker annotation. Check repeatmasker regions bed file format. Download if not provided. @@ -45,11 +43,11 @@ def makeRmsk( Args: regions (str): Path to repeatmasker bed file. - Takes priority over genome. + Ignores "genome", "genomes" and "locus" parameters. genome (str): Genome assembly name. genomes (dict): Dictionary of genome assembly names and URLs. outdir (str): Path to output directory. - locus (bool): If True, prepare for locus-level quantification. + locus (str): If not disabled, prepare for locus-level quantification. outname (str): Name of the output repeatmasker bed file. Returns: @@ -59,8 +57,9 @@ def makeRmsk( SystemExit: If neither regions nor genome is provided, or if the regions file is not properly formatted. """ - # if a repeatmasker bed file is provided, use that + if regions: + # if a repeatmasker bed file is provided, use that is_gz = testGz(regions) f = gzip.open(regions, "rb") if is_gz else open(regions, "r") @@ -75,15 +74,18 @@ def rl(x, decode=False): 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. locus or subfamily)" - " in 4th column.", + " 4 columns with the 4th column being the TE feature name, " + "such as the TE subfamily or locus " + "(as in fragment or instance ID).", error=True, ) f.close() - out = regions - # if no repeatmasker file is provided, and a genome assembly name is - # provided, download and prepare a rmsk.bed file + writerr(f"Using provided RepeatMasker annotation from {regions}.") + return regions + elif genome: + # if no repeatmasker file is provided, and a genome assembly name is + # provided, download and prepare a rmsk.bed file url, header_lines = genomes[genome] writerr( "Downloading and parsing RepeatMasker annotation for " @@ -97,18 +99,20 @@ def rl(x, decode=False): f"Couldn't connect to host.\n\n{e}", error=True, ) - rmsk = gzip.open(io.BytesIO(response.content), "rb") out = os.path.join(outdir, outname) - with gzip.GzipFile(out, "wb", mtime=0) as f: + with ( + gzip.open(io.BytesIO(response.content), "rb") as rmsk_in, + gzip.GzipFile(out, "wb", mtime=0) as rmsk_out, + ): # print header - h = ["#chr", "start", "end", "name", "locus_index", "strand"] + h = ["#chr", "start", "end", "name", "repInstance", "strand"] h = "\t".join(h) + "\n" - f.write(h.encode()) + rmsk_out.write(h.encode()) # skip rmsk header for _ in range(header_lines): - next(rmsk) + next(rmsk_in) # parse rmsk - fams_to_skip = [ + repeats_to_skip = [ "Low_complexity", "Simple_repeat", "rRNA", @@ -117,36 +121,82 @@ def rl(x, decode=False): "tRNA", ] subfamilies = defaultdict(int) - for line in rmsk: - lst = line.decode("utf-8").strip().split() - strand, repname, famclass = lst[8:11] - if famclass.split("/")[0] in fams_to_skip: + rmsk_dict = dict() + for line in rmsk_in: + fields = line.decode("utf-8").strip().split() + chr, start, end = fields[4:7] + strand, subfamily, clasfam = fields[8:11] + instance = fields[14] + clas, family = clasfam.split("/") if "/" in clasfam else (clasfam, "") + # skip problematic repeats + if clas in repeats_to_skip: continue # concatenate family and class with subfamily - repname += "#" + famclass - subfamilies[repname] += 1 - locus_index = subfamilies[repname] - if locus: - # make unique locus names - repname += f"~{locus_index}" - chr, start, end = lst[4:7] - # make coordinates 0-based - start = str(int(start) - 1) + name = f"{subfamily}#{clasfam}" + + # repeat names for locus-level quantifications + if locus != "disabled": + subfamilies[name] += 1 + locus_index = subfamilies[name] + name += f"~{locus_index}" + if locus == "instance": + if instance in rmsk_dict: + rmsk_dict[instance]["repSubfam"].add(subfamily) + rmsk_dict[instance]["repFam"].add(family) + rmsk_dict[instance]["repClass"].add(clas) + else: + rmsk_dict[instance] = { + "repSubfam": {subfamily}, + "repFam": {family}, + "repClass": {clas}, + } + + start = str(int(start) - 1) # make coordinates 0-based if strand != "+": strand = "-" - outl = "\t".join( - [chr, start, end, repname, str(locus_index), strand] - ) + outl = "\t".join([chr, start, end, name, instance, strand]) outl += "\n" - f.write(outl.encode()) - writerr(f"Wrote RepeatMasker annotation to {out}.") + rmsk_out.write(outl.encode()) + if locus == "instance": + # if locus-level quantification by instance is requested, + # prepare a separate file with instance-level annotation + outname_instance = "rmsk_instance.bed.gz" + out_instance = os.path.join(outdir, outname_instance) + with ( + gzip.open(out, "rb") as rmsk_in, + gzip.GzipFile(out_instance, "wb", mtime=0) as rmsk_out, + ): + next(rmsk_in) # skip header row + # print header + h = ["#chr", "start", "end", "name", "repInstance", "strand"] + h = "\t".join(h) + "\n" + rmsk_out.write(h.encode()) + for line in rmsk_in: + chr, start, end, name, instance, strand = ( + line.decode("utf-8").strip().split("\t") + ) + repSubfam = "|".join(sorted(rmsk_dict[instance]["repSubfam"])) + repFam = "|".join( + sorted(f for f in rmsk_dict[instance]["repFam"] if f) + ) + repFam = ( + "" if not repFam else "/" + repFam + ) # to take into account cases in which repFamily is not annotated + repClass = "|".join(sorted(rmsk_dict[instance]["repClass"])) + repname = f"{repSubfam}#{repClass}{repFam}~{instance}" + outl = "\t".join([chr, start, end, repname, instance, strand]) + outl += "\n" + rmsk_out.write(outl.encode()) + level = "subfamily" if locus == "disabled" else locus + outfile = out if locus != "instance" else out_instance + writerr(f"Wrote RepeatMasker {level}-level annotation to {outfile}.") + return outfile + else: writerr( - "Error: it is mandatory to define either --regions OR " - "--genome parameter.", + "Error: it is mandatory to define either --regions OR --genome parameter.", error=True, ) - return out def prepare_whitelist(whitelist, tmpdir): @@ -242,15 +292,11 @@ def isec( else: 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}"]; ' @@ -300,7 +346,7 @@ def chrcat( threads, outdir, tmpdir, - locus=False, + locus="disabled", bedtools="bedtools", verbose=0, ): @@ -337,7 +383,7 @@ def chrcat( cmd2 = f"zcat {mappings_file} " cmd2 += " | cut -f3 | sed 's/,/\\n/g' | gawk '!x[$1]++ { " cmd2 += ' print $1"\\t"gensub(/#' - cmd2 += "[^~]" if locus else "." + cmd2 += "[^~]" if locus != "disabled" else "." cmd2 += '+/,"",1,$1)"\\tGene Expression" }\' ' cmd2 += f" | LC_ALL=C sort -u | gzip > {features_file} " diff --git a/tests/test.yml b/tests/test.yml index a743008..88cb7d3 100644 --- a/tests/test.yml +++ b/tests/test.yml @@ -28,7 +28,7 @@ - path: "irescue_out/tmp/mappings.tsv.gz" md5sum: d404e6c3123f8cfde7689b5fb7763a89 - path: "irescue_out/rmsk.bed.gz" - md5sum: 169571f538624496a00f189771be2f5e + md5sum: 663134e9bb5519d136201d2bd5d35d63 - name: multi tags: @@ -126,10 +126,10 @@ - path: "irescue_out/tmp/mappings.tsv.gz" md5sum: f5aae354f59bef7a4bdb9cf3c5c8dafe -- name: locus +- name: fragment tags: - locus - command: irescue --keeptmp --dump-ec -vv -b ./tests/data/Aligned.sortedByCoord.out.bam -g test --locus + command: irescue --keeptmp --dump-ec -vv -b ./tests/data/Aligned.sortedByCoord.out.bam -g test --locus-level fragment files: - path: "irescue_out/counts/barcodes.tsv.gz" md5sum: 1a74fa12e65ac1703bbe61282854f151 @@ -142,4 +142,24 @@ - path: "irescue_out/tmp/mappings.tsv.gz" md5sum: d7db121c92c04b36e3cc26ae0223426d - path: "irescue_out/rmsk.bed.gz" - md5sum: 75ed0333b049a02672da8b1379cfa6bc + md5sum: 050b0a484a04e5d466afd54c794e35fd + +- name: instance + tags: + - locus + command: irescue --keeptmp --dump-ec -vv -b ./tests/data/Aligned.sortedByCoord.out.bam -g test --locus-level instance + files: + - path: "irescue_out/counts/barcodes.tsv.gz" + md5sum: 1a74fa12e65ac1703bbe61282854f151 + - path: "irescue_out/counts/features.tsv.gz" + md5sum: 045273e35f85522ec70a7be127fd5d7f + - path: "irescue_out/counts/matrix.mtx.gz" + md5sum: a63c59d45b79bd4bfae1d7ffb26b6c6e + - path: "irescue_out/ec_dump.tsv.gz" + md5sum: dc7a39ce64959f0b83dd53c158041746 + - path: "irescue_out/tmp/mappings.tsv.gz" + md5sum: e0986db4bea3929a8f14eb61f2596a41 + - path: "irescue_out/rmsk.bed.gz" + md5sum: 050b0a484a04e5d466afd54c794e35fd + - path: "irescue_out/rmsk_instance.bed.gz" + md5sum: 7ba1ba91b19292a6d49a34e818470630