<commit_before>#! /usr/bin/env python

from __future__ import print_function
import sys, os
sys.path.append(os.path.join(os.path.dirname(__file__), '..'))
from varpkg.overlap import overlapValid
import argparse
import vcf
import re
import banyan
import collections
import numpy


# Parse command line
parser = argparse.ArgumentParser(description='Filter for somatic SVs.')
parser.add_argument('-v', '--vcf', metavar='variants.vcf', required=True, dest='vcfFile', help='input vcf file (required)')
parser.add_argument('-o', '--out', metavar='out.vcf', required=True, dest='outFile', help='output vcf file (required)')
parser.add_argument('-t', '--type', metavar='DEL', required=True, dest='svType', help='SV type [DEL, DUP, INV, INS, TRA] (required)')
parser.add_argument('-a', '--altaf', metavar='0.1', required=False, dest='altAF', help='min. alt. AF (optional)')
parser.add_argument('-m', '--minsize', metavar='500', required=False, dest='minSize', help='min. size (optional)')
parser.add_argument('-n', '--maxsize', metavar='500000000', required=False, dest='maxSize', help='max. size (optional)')
parser.add_argument('-r', '--ratioGeno', metavar='0.75', required=False, dest='ratioGeno', help='min. fraction of genotyped samples (optional)')
parser.add_argument('-f', '--filter', dest='siteFilter', action='store_true', help='Filter sites for PASS')
args = parser.parse_args()

# Command-line args
minSize = 500
if args.minSize:
    minSize = int(args.minSize)
maxSize = 500000000
if args.maxSize:
    maxSize = int(args.maxSize)
altAF = 0.1
if args.altAF:
    altAF = float(args.altAF)
ratioGeno = 0.75
if args.ratioGeno:
    ratioGeno = float(args.ratioGeno)
traWindow = 2500  # 2.5kb translocation window


# Collect high-quality SVs
sv = dict()
rdRat = dict()
svDups = collections.defaultdict(list)
validRecordID = set()
if args.vcfFile:
    vcf_reader = vcf.Reader(open(args.vcfFile), 'r', compressed=True) if args.vcfFile.endswith('.gz') else vcf.Reader(open(args.vcfFile), 'r', compressed=False)
    for record in vcf_reader:
        svlen = record.INFO['END'] - record.POS
        if ((args.svType == 'TRA') or ((svlen >= minSize) and (svlen <= maxSize))) and ((not args.siteFilter) or (len(record.FILTER) == 0)):
            precise = False
            if 'PRECISE' in record.INFO.keys():
                precise = record.INFO['PRECISE']
            rcRef = []
            rcAlt = []
            nCount = 0
            tCount = 0
            for call in record.samples:
                if call.called:
                    if re.search(r"[Nn]ormal", call.sample) != None:
                        nCount += 1
                        if call.gt_type == 0:
                            if ((not precise) and (call['DV'] == 0)) or ((precise) and (call['RV'] == 0)):
                                rcRef.append(call['RC'])
                    if re.search(r"[Tt]umo[ur]", call.sample) != None:
                        tCount += 1
                        if call.gt_type != 0:  # Misses subclonal SVs with call['DV']>0 and GT=0/0
                            if ((not precise) and (call['DV'] >= 2) and (float(call['DV'])/float(call['DV']+call['DR']) >= altAF)) or ((precise) and (call['RV'] >= 2) and (float(call['RV'])/float(call['RR'] + call['RV']) >= altAF)):
                                rcAlt.append(call['RC'])
                else:
                    if (re.search(r"[Nn]ormal", call.sample) != None) and (precise):
                        if (call['DR'] > 2) and (call['DV'] == 0):
                            nCount += 1
                            rcRef.append(call['RC'])
            genotypeRatio = float(nCount + tCount) /  float(len(record.samples))
            if (nCount > 0) and (tCount > 0) and (len(rcRef) == nCount) and (len(rcAlt) > 0) and (genotypeRatio >= ratioGeno):
                rdRatio = 1
                if numpy.median(rcRef):
                    rdRatio = round(numpy.median(rcAlt)/numpy.median(rcRef), 4)
                rdRat[record.ID] = rdRatio
                validRecordID.add(record.ID)
                if not sv.has_key(record.CHROM):
                    sv[record.CHROM] = banyan.SortedDict(key_type=(int, int), alg=banyan.RED_BLACK_TREE, updator=banyan.OverlappingIntervalsUpdator)
                if not sv.has_key(record.INFO['CHR2']):
                    sv[record.INFO['CHR2']] = banyan.SortedDict(key_type=(int, int), alg=banyan.RED_BLACK_TREE, updator=banyan.OverlappingIntervalsUpdator)
                if args.svType == 'TRA':
                    if (record.POS - traWindow, record.POS + traWindow) in sv[record.CHROM]:
                        svDups[(record.CHROM, record.POS - traWindow, record.POS + traWindow)].append((record.ID, record.INFO['PE']))
                    else:<commit_msg>Add arguments to specify sample names for tumor/normal to somaticFilter.py<commit_after>#! /usr/bin/env python

from __future__ import print_function
import sys, os
sys.path.append(os.path.join(os.path.dirname(__file__), '..'))
from varpkg.overlap import overlapValid
import argparse
import vcf
import re
import banyan
import collections
import numpy


# Parse command line
parser = argparse.ArgumentParser(description='Filter for somatic SVs.')
parser.add_argument('-v', '--vcf', metavar='variants.vcf', required=True, dest='vcfFile', help='input vcf file (required)')
parser.add_argument('-o', '--out', metavar='out.vcf', required=True, dest='outFile', help='output vcf file (required)')
parser.add_argument('-t', '--type', metavar='DEL', required=True, dest='svType', help='SV type [DEL, DUP, INV, INS, TRA] (required)')
parser.add_argument('-a', '--altaf', metavar='0.1', required=False, dest='altAF', help='min. alt. AF (optional)')
parser.add_argument('-m', '--minsize', metavar='500', required=False, dest='minSize', help='min. size (optional)')
parser.add_argument('-n', '--maxsize', metavar='500000000', required=False, dest='maxSize', help='max. size (optional)')
parser.add_argument('-r', '--ratioGeno', metavar='0.75', required=False, dest='ratioGeno', help='min. fraction of genotyped samples (optional)')
parser.add_argument('-f', '--filter', dest='siteFilter', action='store_true', help='Filter sites for PASS')

# edits (meiers):
parser.add_argument('-N', '--normal', metavar='name', required=False, dest='nameNormal', help='Name of normal sample in VCF (optional)')
parser.add_argument('-T', '--tumor',  metavar='name', required=False, dest='nameTumor', help='Name of tumor sample in VCF (optional)')

args = parser.parse_args()

# Command-line args
minSize = 500
if args.minSize:
    minSize = int(args.minSize)
maxSize = 500000000
if args.maxSize:
    maxSize = int(args.maxSize)
altAF = 0.1
if args.altAF:
    altAF = float(args.altAF)
ratioGeno = 0.75
if args.ratioGeno:
    ratioGeno = float(args.ratioGeno)
traWindow = 2500  # 2.5kb translocation window


# Collect high-quality SVs
sv = dict()
rdRat = dict()
svDups = collections.defaultdict(list)
validRecordID = set()
if args.vcfFile:
    vcf_reader = vcf.Reader(open(args.vcfFile), 'r', compressed=True) if args.vcfFile.endswith('.gz') else vcf.Reader(open(args.vcfFile), 'r', compressed=False)
    for record in vcf_reader:
        svlen = record.INFO['END'] - record.POS
        if ((args.svType == 'TRA') or ((svlen >= minSize) and (svlen <= maxSize))) and ((not args.siteFilter) or (len(record.FILTER) == 0)):
            precise = False
            if 'PRECISE' in record.INFO.keys():
                precise = record.INFO['PRECISE']
            rcRef = []
            rcAlt = []
            nCount = 0
            tCount = 0
            for call in record.samples:
                if call.called:
                    if args.nameNormal and call.sample == args.nameNormal or re.search(r"[Nn]ormal", call.sample) != None:
                        nCount += 1
                        if call.gt_type == 0:
                            if ((not precise) and (call['DV'] == 0)) or ((precise) and (call['RV'] == 0)):
                                rcRef.append(call['RC'])
                    if args.nameTumor and call.sample == args.nameTumor or re.search(r"[Tt]umo[ur]", call.sample) != None:
                        tCount += 1
                        if call.gt_type != 0:  # Misses subclonal SVs with call['DV']>0 and GT=0/0
                            if ((not precise) and (call['DV'] >= 2) and (float(call['DV'])/float(call['DV']+call['DR']) >= altAF)) or ((precise) and (call['RV'] >= 2) and (float(call['RV'])/float(call['RR'] + call['RV']) >= altAF)):
                                rcAlt.append(call['RC'])
                else:
                    if (args.nameNormal and call.sample == args.nameNormal or re.search(r"[Nn]ormal", call.sample) != None) and (precise):
                        if (call['DR'] > 2) and (call['DV'] == 0):
                            nCount += 1
                            rcRef.append(call['RC'])
            genotypeRatio = float(nCount + tCount) /  float(len(record.samples))
            if (nCount > 0) and (tCount > 0) and (len(rcRef) == nCount) and (len(rcAlt) > 0) and (genotypeRatio >= ratioGeno):
                rdRatio = 1
                if numpy.median(rcRef):
                    rdRatio = round(numpy.median(rcAlt)/numpy.median(rcRef), 4)
                rdRat[record.ID] = rdRatio
                validRecordID.add(record.ID)
                if not sv.has_key(record.CHROM):
                    sv[record.CHROM] = banyan.SortedDict(key_type=(int, int), alg=banyan.RED_BLACK_TREE, updator=banyan.OverlappingIntervalsUpdator)
                if not sv.has_key(record.INFO['CHR2']):
                    sv[record.INFO['CHR2']] = banyan.SortedDict(key_type=(int, int), alg=banyan.RED_BLACK_TREE, updator=banyan.OverlappingIntervalsUpdator)
                if args.svType == 'TRA':
                    if (record.POS - traWindow, record.POS + traWindow) in sv[record.CHROM]:
                        svDups[(record.CHROM, record.POS - traWindow, record.POS + traWindow)].append((record.ID, record.INFO['PE']))
                    else:
{'fraction_non_alphanumeric': 0.15942897074972548, 'fraction_numerical': 0.012778276929220326, 'mean_word_length': 2.620527647271413, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1364199', 'n_tokens_mistral': 3391, 'n_tokens_neox': 3170, 'n_words': 733}