#! /usr/bin/env python
import argparse, sys
import os.path
from math import fabs
import numpy
import pysam
import vcf
from svreader.vcf_utils import get_template
from svreader import SVInter, SVRecord, SVReader
# A small wrapper to print to stderr
def eprint(*args, **kwargs):
print(*args, file=sys.stderr, **kwargs)
class CNVR(SVInter):
'''
CNVR object :
Copy Number Variation Region
'''
def __init__(self,svs):
svinter = next(iter(svs))
start = int(numpy.median(numpy.array([x.start for x in svs])))
end = int(numpy.median(numpy.array([x.end for x in svs])) )
super(CNVR, self).__init__(svinter.chrom,start,end,".")
self.__svs = {x.id:x for x in svs}
self._repr_cnv = self._repr_cnv()
self._bp_precision()
self._CIntervals()
if self._repr_cnv:
eprint(self._repr_cnv.id,"Precise")
self.start = self._repr_cnv.record.pos
self.end = self._repr_cnv.record.info['END']
@property
def svs(self):
return self.__svs
def __str__(self):
return self.chrom+":"+str(self.start)+"-"+str(self.end)+"\t"+self.id+"\t"+str(self.length())
def overlaps(self):
return [ self.svs[s].length()/float(self.length()) for s in sorted(self.svs)]
def intervals(self):
return [ self.svs[s].id for s in sorted(self.svs) ]
def CopyNumbers(self):
svs = self.__svs
cn = []
for sv in svs:
cn.append(svs[sv].CN())
return ",".join(cn)
def NumCNV(self):
return len(self.__svs)
def Callers(self):
callers=list(set([c.split('_')[0] for c in self.__svs]))
return callers
def CallersVsamples(self):
vsamples = []
for caller in self.Callers():
cnv=self.gettoolCNV(caller)
if 'VSAMPLES' in cnv.record.info:
vsamples.append(caller+":"+'-'.join(cnv.record.info['VSAMPLES']))
return vsamples
def NumCallers(self):
return len(set([c.split('_')[0] for c in self.__svs]))
def CNV(self):
return list(self.__svs.keys())
def cipos(self):
return ",".join([str(x) for x in self._cipos])
def ciend(self):
return ",".join([str(x) for x in self._ciend])
def bedformat(self):
return "\t".join(map(str,[self.chrom,self.start,self.end,self.name,len(self.__svs),".",",".join(self.intervals()),",".join(["%0.2f" % x for x in self.overlaps()]),]))
def precision(self):
return self._precision
def IsPrecise(self):
if self._repr_cnv:
return True
else:
return False
def repr_cnv(self):
return self._repr_cnv.record.id
def _bp_precision(self):
# Trying to infer the precision breakpoints according to the the different CNVs merged
svs = self.__svs
#print(self)
#print(svs.keys())
starts = [ svs[x].start for x in svs ]
ends = [ svs[x].end for x in svs ]
start_range = max(starts) - min(starts)
end_range = max(ends) - min(ends)
#print(",".join([ str(svs[x].start) for x in svs]),start_range)
#print(",".join([ str(svs[x].end) for x in svs]),end_range)
self._precision = [start_range,end_range]
def _repr_cnv(self):
'''
Trying to identify a cnv among the cnv of this CNVR that is a representative
either Pindel (more than 10 supp reads), or PRECISE, delly, lumpy or genomeSTRIP
'''
selected = None
if "pindel" in self.Callers():
cnv = self.gettoolCNV("pindel")
if cnv.record.info["SU"][0]>10:
selected = cnv
elif "delly" in self.Callers():
cnv = self.gettoolCNV("delly")
if "PRECISE" in cnv.record.info and cnv.record.info["SR"]>10:
selected = cnv
elif "lumpy" in self.Callers():
cnv = self.gettoolCNV("lumpy")
if not "IMPRECISE" in cnv.record.info and cnv.record.info["SR"][0]>10:
selected = cnv
elif "genomestrip" in self.Callers():
cnv = self.gettoolCNV("genomestrip")
if not "IMPRECISE" in cnv.record.info:
selected = cnv
return selected
def _CIntervals(self):
# Trying to infer CIPOS and CIEND from the called CNVs
# we search for the CNV CIPOS and CIEND from the most confident
# to the least confident SV detection tool
if self._repr_cnv:
self._cipos = self._repr_cnv.record.info['CIPOS']
self._ciend = self._repr_cnv.record.info['CIEND']
return
# default is large confidence interval
self._cipos = (-50,50)
self._ciend = (-50,50)
tools=["genomestrip","lumpy","delly","pindel"]
for tool in tools:
cnv = self.gettoolCNV(tool)
if cnv:
self._cipos = cnv.record.info['CIPOS']
self._ciend = cnv.record.info['CIEND']
break # if a cnv from that tool was found we keep the walues from this call
def gettoolCNV(self,toolname):
cnvs = self.CNV()
# Find CNV made by a specific tool
indices = [i for i,s in enumerate(cnvs) if toolname in s]
if len(indices):
return self.__svs[cnvs[indices[0]]]
else:
return None
def to_vcf_record(cnv):
# Copy numbers
info = {"SVLEN" : cnv.length(), "SVTYPE" : "DEL", "END" : cnv.end, "NB_CNV" : cnv.NumCNV(), "NAMES":cnv.CNV(),"PRECISION":",".join(self._precision)}
alt = [vcf.model._SV("DEL")]
vcf_record = vcf.model._Record(cnv.chrom,
cnv.start,
cnv.name,
"N",
alt,
".",
".",
info,
"",
[0],
[])
return vcf_record
def trim_column_header(vcf_reader):
vcf_reader._column_headers.pop()
def passed_variant(record):
"""Did this variant pass?"""
return record.filter is None or len(record.filter) == 0
def intersect(a,b):
if a.chrom != b.chrom:
return 0
(low,high) = (a,b) if a.start< b.start else (b,a)
intersect = 0
if low.end < high.start: # low << high
intersect = 0
elif low.end > high.end: # low includes high
intersect = high.end - high.start
else: # high.start is within low.start and low.end
intersect = low.end - high.start
return intersect
def roverlap(a,b,cutoff):
'''
Returns true if the two intervals overlap (> cutoff and reciprocally)
'''
inter = intersect(a,b)
if inter:
pa = float(inter)/a.length()
pb = float(inter)/b.length()
else:
pa=0
pb=0
return True if pa >= cutoff and pb >= cutoff else False
def breakpoint_left_precision(a,b):
'''
Returns the precision of the left breakpoint
'''
left_precision = fabs(a.start-b.start)
return left_precision
def breakpoint_right_precision(a,b):
'''
Returns the precision of the right breakpoint
'''
right_precision = fabs(a.end-b.end)
return right_precision
def construct_overlap_graph(data,cutoff,left_precision,right_precision):
'''
Graph construction : add a link whenever two SV overlap
'''
nodes = set()
for i in range(len(data)):
for j in range(i+1,len(data)):
if ( roverlap(data[i],data[j],cutoff) and
breakpoint_left_precision(data[i],data[j]) <= left_precision and
breakpoint_right_precision(data[i],data[j]) <= right_precision ):
#print("%s %s %d %d %d %d %d" % (data[i],data[j],data[i].length(),data[j].length(),intersect(data[i],data[j]),breakpoint_left_precision(data[i],data[j]),breakpoint_right_precision(data[i],data[j])))
data[i].add_link(data[j])
def connected_components(nodes):
'''
Constructing connected component of the graph
input : a set of nodes with linkage information
output : a set of groups, one group for each connected component
'''
# List of connected components found. The order is random.
result = []
# Make a copy of the set, so we can modify it.
nodes = set(nodes)
# Iterate while we still have nodes to process.
while nodes:
n = nodes.pop() # Get a random node and remove it from the global set.
group = {n} # This set will contain the next group of nodes connected to each other.
queue = [n] # Build a queue with this node in it.
# Iterate the queue.
# When it's empty, we finished visiting a group of connected nodes.
while queue:
n = queue.pop(0) # Consume the next item from the queue.
neighbors = n.links # Fetch the neighbors.
neighbors.difference_update(group) # Remove the neighbors we already visited.
nodes.difference_update(neighbors) # Remove the remaining nodes from the global set.
group.update(neighbors) # Add them to the group of connected nodes.
queue.extend(neighbors) # Add them to the queue, so we visit them in the next iterations.
result.append(group) # Add the group to the list of groups.
# Return the list of gimport vcfroups.
return result
# --------------------------------------
# main function
def main():
# parse the command line args
args = get_args()
infiles = args.vcf
prefix = args.prefix
no_index = args.no_index
overlap_cutoff = args.overlap_cutoff
left_precision = args.left_precision
right_precision = args.right_precision
filenames=infiles.split(",")
# Checking the existence of the files
# Reading all the vcf files
SVSet=[]
for infile in filenames:
eprint(" Reading file %s" % (infile))
for record in SVReader(infile):
if not passed_variant(record):
continue
SVSet.append(record)
# Computing connected components according to reciprocal overlaps, left and right precision
# reciprocal overlaps : only SV intervals with reciprocal overlap > = overlap_cutoff are linked
# left, right precision : only SV intervals with left (right) within left_precision are linked
eprint("Constructing conected component")
construct_overlap_graph(SVSet,overlap_cutoff,left_precision,right_precision)
number = 1
cnvr = []
for components in connected_components(SVSet):
names = sorted(node.name for node in components)
names = ", ".join(names)
eprint("Group #%i: %s" % (number, names))
cnv = CNVR(components)
cnvr.append(cnv)
number += 1
# Writing the merge file in a single vcf file
vcf_template_reader = get_template( "merge" )
vcf_writer = vcf.Writer( sys.stdout, vcf_template_reader)
number = 1
prefix_name = prefix.split(".")[0]
for cnv in sorted(cnvr, key=lambda k: k.start):
cnv.name = prefix_name+"_"+str(number)
record = to_vcf_record(cnv)
vcf_writer.write_record(record)
number += 1
vcf_writer.close()
# initialize the script
if __name__ == '__main__':
try:
sys.exit(main())
except:
raise