convert_stripmask_to_bed.py

#!/usr/bin/env python3

# This script takes a mask file generated by Genome STRiP ComputeGenomeMask:
# http://gatkforums.broadinstitute.org/gatk/discussion/1499/computegenomemask
# and outputs a BED3 file containing the masked regions to stdout
 
import sys

CURRENT_CHROMOSOME = None
CURRENT_POSITION = 0
IN_MASK_INTERVAL = False
IN_MASK_CHAR = "1"
NOT_IN_MASK_CHAR = "0"
CURRENT_INTERVAL = None

def main():
    with open(sys.argv[1], "r") as input_file:
        for line in input_file:
            line = line.rstrip()
            if line.startswith(">"):
                start_new_chromosome(line)
            else:
                parse_mask_line(line)
        if IN_MASK_INTERVAL:
            exit_mask_interval_and_print()

def start_new_chromosome(line):
    global CURRENT_CHROMOSOME
    global CURRENT_POSITION
    if IN_MASK_INTERVAL:
        exit_mask_interval_and_print()
    CURRENT_CHROMOSOME = get_chromosome_name_from_line(line)
    CURRENT_POSITION = 0

def get_chromosome_name_from_line(line):
    return line.rstrip()[1:]

def parse_mask_line(line):
    for char in line:
        parse_mask_char(char)

def parse_mask_char(char):
    global CURRENT_POSITION
    if char == NOT_IN_MASK_CHAR:
        if IN_MASK_INTERVAL:
            exit_mask_interval_and_print()
    else:
        assert char == IN_MASK_CHAR
        if not IN_MASK_INTERVAL:
            enter_mask_interval()
    CURRENT_POSITION += 1

def enter_mask_interval():
    global IN_MASK_INTERVAL
    global CURRENT_INTERVAL
    assert not IN_MASK_INTERVAL
    IN_MASK_INTERVAL = True
    CURRENT_INTERVAL = CURRENT_CHROMOSOME + "\t" + str(CURRENT_POSITION)

def exit_mask_interval_and_print():
    global IN_MASK_INTERVAL
    global CURRENT_INTERVAL
    assert IN_MASK_INTERVAL
    IN_MASK_INTERVAL = False
    CURRENT_INTERVAL += ("\t" + str(CURRENT_POSITION))
    print(CURRENT_INTERVAL)

if __name__ == "__main__":
    main()