bedpull 0.3.0

bedpull - Pull the query sequence from bam or fasta references using a bed file
Documentation 
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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

#!/usr/bin/env python3
"""Turn a bedpull --paf mode FASTA + --unmapped output into a results.tsv.

PAF-mode headers look like:
  contig_name|chr:region_start-region_end|region_name|contig_name:query_start-query_end|strand
  contig_name|chr:region_start-region_end|region_name|contig_name:query_start-query_end|strand|h1

(--partial, and its missing_left/missing_right header suffix, is BAM/CRAM-only
and never appears in PAF-mode output, so this parser doesn't need to handle it.)

A window (region_name) can produce zero, one, or several FASTA records — zero
if --unmapped logged it, several if more than one alignment overlapped the
window. This script keeps every hit; the scorer picks the best one per name
across haplotype runs.
"""

import argparse
import sys


def parse_fasta_lengths(fasta_path):
    """Yield (header, sequence_length) for each record in a FASTA file."""
    header = None
    seq_len = 0
    with open(fasta_path) as f:
        for line in f:
            line = line.rstrip("\n")
            if line.startswith(">"):
                if header is not None:
                    yield header, seq_len
                header = line[1:]
                seq_len = 0
            else:
                seq_len += len(line)
    if header is not None:
        yield header, seq_len


def parse_unmapped_bed(unmapped_path):
    """Yield (reason, chrom, start, end, name) from a bedpull --unmapped file."""
    reason = None
    try:
        f = open(unmapped_path)
    except FileNotFoundError:
        return
    with f:
        for line in f:
            line = line.rstrip("\n")
            if not line:
                continue
            if line.startswith("#"):
                reason = line[1:]
                continue
            fields = line.split("\t")
            if len(fields) < 4:
                continue
            yield reason, fields[0], int(fields[1]), int(fields[2]), fields[3]


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument("--fasta", required=True, help="bedpull FASTA output")
    p.add_argument("--unmapped", required=True, help="bedpull --unmapped output (may not exist)")
    p.add_argument("--hap", required=True, help="haplotype tag to stamp on every row, e.g. PATERNAL")
    p.add_argument("--out", required=True, help="output results.tsv path")
    args = p.parse_args()

    cols = [
        "name",
        "hap",
        "status",
        "unmapped_reason",
        "extracted_len",
        "contig_name",
        "query_start",
        "query_end",
        "strand",
    ]

    rows = []
    for header, seq_len in parse_fasta_lengths(args.fasta):
        fields = header.split("|")
        if len(fields) < 5:
            print(f"warning: unrecognised header shape, skipping: {header}", file=sys.stderr)
            continue
        contig_name = fields[0]
        region_name = fields[2]
        query_coords = fields[3]
        strand = fields[4][:1]
        query_start, query_end = "", ""
        if ":" in query_coords:
            _, coord_part = query_coords.split(":", 1)
            if "-" in coord_part:
                query_start, query_end = coord_part.split("-", 1)
        rows.append(
            {
                "name": region_name,
                "hap": args.hap,
                "status": "ok",
                "unmapped_reason": "",
                "extracted_len": seq_len,
                "contig_name": contig_name,
                "query_start": query_start,
                "query_end": query_end,
                "strand": strand,
            }
        )

    for reason, _chrom, _start, _end, name in parse_unmapped_bed(args.unmapped):
        rows.append(
            {
                "name": name,
                "hap": args.hap,
                "status": "unmapped",
                "unmapped_reason": reason or "",
                "extracted_len": 0,
                "contig_name": "",
                "query_start": "",
                "query_end": "",
                "strand": "",
            }
        )

    with open(args.out, "w") as out:
        out.write("\t".join(cols) + "\n")
        for r in rows:
            out.write("\t".join(str(r[c]) for c in cols) + "\n")

    n_ok = sum(1 for r in rows if r["status"] == "ok")
    n_unmapped = sum(1 for r in rows if r["status"] == "unmapped")
    print(f"{args.hap}: {n_ok} hit(s), {n_unmapped} unmapped window(s)", file=sys.stderr)


if __name__ == "__main__":
    main()