vcfs_comp.py

import argparse
from pathlib import Path

import vcf


def save_combination(v1, v2):
    print(f"Similar SV found at {v1.CHROM}:{v1.POS} ({v1.INFO['SVTYPE']}) vs. \
            {v2.CHROM}:{v2.POS} ({v2.INFO['SVTYPE']})")

def make_reader(f: Path) -> vcf.Reader:
    if f.suffix == ".gz":
        return vcf.Reader(open(f, 'rb'), compressed=True)
    else:
        return vcf.Reader(open(f, 'r'))

def compare(vcf1, vcf2, t_lower=50, t_upper=100):
    r1, r2 = make_reader(vcf1), make_reader(vcf2)

    var1, var2 = next(r1), next(r2)

    while True:
        if var1.CHROM < var2.CHROM:
            var1 = next(r1)
            continue
        elif var1.CHROM > var2.CHROM:
            var2 = next(r2)
            continue

        dist = var1.POS - var2.POS
        if t_lower < abs(dist) < t_upper:
            save_combination(var1, var2)

        try:
            if var1 < var2:
                var1 = next(r1)
            else:
                var2 = next(r2)
        except StopIteration:
            print("Iteration stopped, as one of the files ended!")
            break


if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument("vcf1", help="First VCF file", type=str)
    parser.add_argument("vcf2", help="Second VCF file", type=str)
    parser.add_argument("-l", "--lower_threshold",
                        help="The minimum distance for reported common SVs",
                        type=int)
    parser.add_argument("-u", "--upper_threshold",
                        help="The maximum distance for reported common SVs",
                        type=int)

    args = parser.parse_args()
    print(f"You have chosen to compare the files {args.vcf1} and {args.vcf2}!")

    compare(args.vcf1, args.vcf2, args.lower_threshold, args.upper_threshold)