#!/usr/bin/env python

from __future__ import print_function

import logging

import os

import random

import re

import sys

from argparse import SUPPRESS

import matplotlib

matplotlib.use("Agg") #must be before imports of submodules in matplotlib

import matplotlib.gridspec as gridspec

import matplotlib.patches as mpatches

import matplotlib.path as mpath

import matplotlib.pyplot as plt

import matplotlib.ticker as ticker

import numpy as np

import pysam

import warnings

warnings.filterwarnings('ignore', 'FixedFormatter should only be used together with FixedLocator')

from matplotlib.offsetbox import AnchoredText

logger = logging.getLogger(__name__)

INTERCHROM_YAXIS = 5000

COLORS = {

"Deletion/Normal": "black",

"Deletion": "black",

"Duplication": "red",

"Inversion": "blue",

"InterChrmInversion": "blue",

"InterChrm": "black",

}

READ_TYPES_USED = {

"Deletion/Normal": False,

"Duplication": False,

"Inversion": False,

"Aligned long read": False,

"Linked read": False,

"Split-read": False,

"Paired-end read": False,

}

# pysam.readthedocs.io/en/latest/api.html#pysam.AlignedSegment.cigartuples

CIGAR_MAP = {

"M": 0,

"I": 1,

"D": 2,

"N": 3,

"S": 4,

"H": 5,

"P": 6,

"=": 7,

"X": 8,

"B": 9,

}

def strip_chr(chrom):

"""

safer way to replace chr string, to support non-human genomes

"""

if chrom[:3] == "chr":

chrom = chrom[3:]

return chrom

# {{{class plan_step:

class plan_step:

step_events = ["Align", "ANNOTATION"]

def __init__(self, start_pos, end_pos, event, info=None):

self.start_pos = start_pos

self.end_pos = end_pos

self.event = event

self.info = info

def __str__(self):

if self.info:

return (

"Step("

+ str(self.start_pos)

+ ", "

+ str(self.end_pos)

+ ", "

+ self.event

+ ", "

+ str(self.info)

+ ")"

)

else:

return (

"Step("

+ str(self.start_pos)

+ ", "

+ str(self.end_pos)

+ ", "

+ self.event

+ ")"

)

def __repr__(self):

return str(self)

# }}}

# {{{class genome_interval:

class genome_interval:

def __init__(self, chrm, start, end):

self.chrm = chrm

self.start = start

self.end = end

def __str__(self):

return "(" + self.chrm + "," + str(self.start) + "," + str(self.end) + ")"

def __repr__(self):

return str(self)

def __eq__(self, gi2):

return self.chrm == gi2.chrm and self.start == gi2.start and self.end == gi2.end

""" return -1 if before, 0 if in, 1 if after """

def intersect(self, gi):

if strip_chr(gi.chrm) < strip_chr(self.chrm) or gi.end < self.start:

return -1

elif strip_chr(gi.chrm) > strip_chr(self.chrm) or gi.start > self.end:

return 1

else:

return 0

# }}}

# {{{def get_range_hit(ranges, chrm, point):

def get_range_hit(ranges, chrm, point):

for j in range(len(ranges)):

r = ranges[j]

if (

strip_chr(r.chrm) == strip_chr(chrm)

and r.start <= point

and r.end >= point

):

return j

return None

# }}}

# {{{def map_genome_point_to_range_points(ranges, chrm, point):

def map_genome_point_to_range_points(ranges, chrm, point):

range_hit = get_range_hit(ranges, chrm, point)

if range_hit == None:

return None

p = 1.0 / len(ranges) * range_hit + (1.0 / len(ranges)) * (

float(point - ranges[range_hit].start)

/ float(ranges[range_hit].end - ranges[range_hit].start)

)

return p

# }}}

# {{{def points_in_window(points):

def points_in_window(points):

"""Checks whether these points lie within the window of interest

Points is a list of one start, one end coordinate (ints)

"""

if (

None in points

or points[0] < -5

or points[1] < -5

or points[0] > 5

or points[1] > 5

):

return False

return True

# }}}

# {{{ def get_tabix_iter(chrm, start, end, datafile):

def get_tabix_iter(chrm, start, end, datafile):

"""Gets an iterator from a tabix BED/GFF3 file

Used to avoid chrX vs. X notation issues when extracting data from

annotation files

"""

try:

tbx = pysam.TabixFile(datafile)

except:

tbx = pysam.TabixFile(datafile, index=datafile+".csi")

itr = None

try:

itr = tbx.fetch(chrm, max(0, start - 1000), end + 1000)

except ValueError:

# try and account for chr/no chr prefix

if chrm[:3] == "chr":

chrm = chrm[3:]

else:

chrm = "chr" + chrm

try:

itr = tbx.fetch(chrm, max(0, start - 1000), end + 1000)

except ValueError as e:

logger.warning(

"Could not fetch {}:{}-{} from {}".format(

chrm,

start,

end,

datafile

)

)

print(e)

return itr

# }}}

##Coverage methods

# {{{def add_coverage(bam_file, read, coverage, separate_mqual):

def add_coverage(read, coverage_matrix, offset, column):

"""Adds a read to the known coverage

Coverage from Pysam read is added to coverage_matrix.

offset defines the start position of the current range

column specifies which column to add to.

"""

curr_pos = read.reference_start

if not read.cigartuples:

return

for op, length in read.cigartuples:

if op in [CIGAR_MAP["M"], CIGAR_MAP["="], CIGAR_MAP["X"]]:

coverage_matrix[curr_pos - offset: curr_pos + length - offset, column] += 1

curr_pos += length

elif op == CIGAR_MAP["I"]:

curr_pos = curr_pos

elif op == CIGAR_MAP["D"]:

curr_pos += length

elif op == CIGAR_MAP["N"]:

curr_pos = length

elif op == CIGAR_MAP["S"]:

curr_pos = curr_pos

elif op == CIGAR_MAP["H"]:

curr_pos = curr_pos

else:

curr_pos += length

# }}}

# {{{def plot_coverage(coverage,

def plot_coverage(

coverage,

ax,

ranges,

hp_count,

max_coverage,

tracktype,

yaxis_label_fontsize,

max_coverage_points,

):

"""Plots high and low quality coverage for the region

User may specify a preference between stacked and superimposed

superimposed may cause unexpected behavior if low-quality depth is

greater than high

"""

cover_x = []

cover_y_lowqual = []

cover_y_highqual = []

cover_y_all = []

for i in range(len(ranges)):

r = ranges[i]

region_len = r.end-r.start

downsample = 1

if region_len > max_coverage_points:

downsample = int(region_len / max_coverage_points)

for i,pos in enumerate(range(r.start, r.end + 1)):

if i%downsample != 0:

continue

cover_x.append(map_genome_point_to_range_points(ranges, r.chrm, pos))

if r.chrm in coverage and pos in coverage[r.chrm]:

cover_y_all.append(coverage[r.chrm][pos][0] + coverage[r.chrm][pos][1])

cover_y_highqual.append(coverage[r.chrm][pos][0])

cover_y_lowqual.append(coverage[r.chrm][pos][1])

else:

cover_y_lowqual.append(0)

cover_y_highqual.append(0)

cover_y_all.append(0)

cover_y_lowqual = np.array(cover_y_lowqual)

cover_y_highqual = np.array(cover_y_highqual)

cover_y_all = np.array(cover_y_all)

if max_coverage > 0:

max_plot_depth = max_coverage

elif cover_y_all.max() > 3 * cover_y_all.mean():

max_plot_depth = max(

np.percentile(cover_y_all, 99.5), np.percentile(cover_y_all, 99.5)

)

else:

max_plot_depth = np.percentile(cover_y_all.max(), 99.5)

ax2 = ax.twinx()

ax2.set_xlim([0, 1])

if 0 == max_plot_depth:

max_plot_depth = 0.01

ax2.set_ylim([0, max(1, max_plot_depth)])

bottom_fill = np.zeros(len(cover_y_all))

if tracktype == "stack":

ax2.fill_between(

cover_x,

cover_y_highqual,

bottom_fill,

color="darkgrey",

step="pre",

alpha=0.4,

)

ax2.fill_between(

cover_x, cover_y_all, cover_y_highqual, color="grey", step="pre", alpha=0.15

)

elif tracktype == "superimpose":

ax2.fill_between(

cover_x, cover_y_lowqual, bottom_fill, color="grey", step="pre", alpha=0.15

)

ax2.fill_between(

cover_x,

cover_y_highqual,

cover_y_lowqual,

color="darkgrey",

step="pre",

alpha=0.4,

)

ax2.fill_between(

cover_x, cover_y_lowqual, bottom_fill, color="grey", step="pre", alpha=0.15

)

## tracktype==None also allowed

# number of ticks should be 6 if there's one hp, 3 otherwise

tick_count = 5 if hp_count == 1 else 2

tick_count = max(int(max_plot_depth / tick_count), 1)

# set axis parameters

#ax2.yaxis.set_major_locator(ticker.FixedLocator(tick_count))

ax2.yaxis.set_major_locator(ticker.MultipleLocator(tick_count))

ax2.tick_params(axis="y", colors="grey", labelsize=yaxis_label_fontsize)

ax2.spines["top"].set_visible(False)

ax2.spines["bottom"].set_visible(False)

ax2.spines["left"].set_visible(False)

ax2.spines["right"].set_visible(False)

ax2.tick_params(axis="x", length=0)

ax2.tick_params(axis="y", length=0)

# break the variant plot when we have multiple ranges

for i in range(1, len(ranges)):

ax2.axvline(x=1.0 / len(ranges), color="white", linewidth=5)

return ax2

# }}}

##Pair End methods

# {{{class PairedEnd:

class PairedEnd:

"""container of paired-end read info

Contains start(int), end(int), strand(bool True=forward), MI (int

molecular identifier), HP (int haplotype)

"""

def __init__(self, chrm, start, end, is_reverse, MI_tag, HP_tag):

"""Create PairedEnd instance

Genomic interval is defined by start and end integers

Strand is opposite of is_reverse

Molecular identifier and Haplotype are integers if present, else

False

"""

self.pos = genome_interval(chrm, start, end)

self.strand = not (is_reverse)

# molecular identifier - linked reads only

self.MI = None

# haplotype - phased reads only

self.HP = 0

if MI_tag:

self.MI = MI_tag

if HP_tag:

self.HP = HP_tag

def __repr__(self):

return "PairedEnd(%s,%s,%s,%s,%s,%s)" % (

self.pos.chrm,

self.pos.start,

self.pos.end,

self.strand,

self.MI,

self.HP,

)

# }}}

# {{{ def add_pair_end(bam_file, read, pairs, linked_reads):

def add_pair_end(bam_file, read, pairs, linked_reads, ignore_hp):

"""adds a (mapped, primary, non-supplementary, and paired) read to the

pairs list

Pysam read is added as simpified PairedEnd instance to pairs

Also added to linked_reads list if there is an associated MI tag

"""

if read.is_unmapped:

return

if not (read.is_paired):

return

if read.is_secondary:

return

if read.is_supplementary:

return

MI_tag = False

HP_tag = False

if read.has_tag("MI"):

MI_tag = int(read.get_tag("MI"))

if not ignore_hp and read.has_tag("HP"):

HP_tag = int(read.get_tag("HP"))

pe = PairedEnd(

bam_file.get_reference_name(read.reference_id),

read.reference_start,

read.reference_end,

read.is_reverse,

MI_tag,

HP_tag,

)

if pe.HP not in pairs:

pairs[pe.HP] = {}

if read.query_name not in pairs[pe.HP]:

pairs[pe.HP][read.query_name] = []

if pe.MI:

if pe.HP not in linked_reads:

linked_reads[pe.HP] = {}

if pe.MI not in linked_reads[pe.HP]:

linked_reads[pe.HP][pe.MI] = [[], []]

linked_reads[pe.HP][pe.MI][0].append(read.query_name)

pairs[pe.HP][read.query_name].append(pe)

pairs[pe.HP][read.query_name].sort(key=lambda x: x.pos.start)

# }}}

# {{{def sample_normal(max_depth, pairs, z):

def sample_normal(max_depth, pairs, z):

"""Downsamples paired-end reads

Selects max_depth reads

Does not remove discordant pairs, those with insert distance greater

than z stdevs from mean

Returns downsampled pairs list

"""

sampled_pairs = {}

plus_minus_pairs = {}

if max_depth == 0:

return sampled_pairs

for read_name in pairs:

pair = pairs[read_name]

if len(pair) != 2:

continue

if pair[0].strand == True and pair[1].strand == False:

plus_minus_pairs[read_name] = pair

else:

sampled_pairs[read_name] = pair

if len(plus_minus_pairs) > max_depth:

lens = np.array(

[pair[1].pos.end - pair[0].pos.start for pair in plus_minus_pairs.values()]

)

mean = np.mean(lens)

stdev = np.std(lens)

inside_norm = {}

for read_name in pairs:

pair = pairs[read_name]

if len(pair) != 2:

continue

if pair[1].pos.end - pair[0].pos.start >= mean + z * stdev:

sampled_pairs[read_name] = pair

else:

inside_norm[read_name] = pair

if len(inside_norm) > max_depth:

for read_name in random.sample(list(inside_norm.keys()), max_depth):

sampled_pairs[read_name] = inside_norm[read_name]

else:

for read_name in inside_norm:

sampled_pairs[read_name] = inside_norm[read_name]

else:

for read_name in plus_minus_pairs:

sampled_pairs[read_name] = plus_minus_pairs[read_name]

return sampled_pairs

# }}}

# {{{def get_pairs_insert_sizes(pairs):

def get_pairs_insert_sizes(ranges, pairs):

"""Extracts the integer insert sizes for all pairs

Return list of integer insert sizes

"""

pair_insert_sizes = []

for hp in pairs:

for read_name in pairs[hp]:

if len(pairs[hp][read_name]) == 2:

size = get_pair_insert_size(ranges, pairs[hp][read_name])

if size:

pair_insert_sizes.append(size)

return pair_insert_sizes

# }}}

# {{{def get_pair_insert_size(ranges, pair):

def get_pair_insert_size(ranges, pair):

""" Gives the outer distance

"""

first = pair[0]

second = pair[1]

# make sure both sides are in range

if (

get_range_hit(ranges, first.pos.chrm, first.pos.start) != None

or get_range_hit(ranges, first.pos.chrm, first.pos.end) != None

) and (

get_range_hit(ranges, second.pos.chrm, second.pos.start) != None

or get_range_hit(ranges, second.pos.chrm, second.pos.end) != None

):

if first.pos.chrm == second.pos.chrm:

return abs(second.pos.end - first.pos.start)

else:

return INTERCHROM_YAXIS

else:

return None

# }}}

# {{{ def get_pairs_plan(ranges, pairs, linked_plan=False):

def get_pairs_plan(ranges, pairs, linked_plan=False):

steps = []

max_event = 0

insert_sizes = []

for read_name in pairs:

pair = pairs[read_name]

plan = get_pair_plan(ranges, pair)

if plan:

insert_size, step = plan

insert_sizes.append(insert_size)

steps.append(step)

if len(insert_sizes) > 0:

max_event = max(insert_sizes)

plan = [max_event, steps]

return plan

# }}}

# {{{def get_pair_plan(ranges, pair, linked_plan=False):

def get_pair_plan(ranges, pair, linked_plan=False):

if pair == None or len(pair) != 2:

return None

first = pair[0]

second = pair[1]

# see if they are part of a linked read

if not linked_plan and (first.MI or second.MI):

return None

# make sure both ends are in the plotted region

first_s_hit = get_range_hit(ranges, first.pos.chrm, first.pos.start)

first_e_hit = get_range_hit(ranges, first.pos.chrm, first.pos.end)

second_s_hit = get_range_hit(ranges, second.pos.chrm, second.pos.start)

second_e_hit = get_range_hit(ranges, second.pos.chrm, second.pos.end)

if (first_s_hit == None and first_e_hit == None) or (

second_s_hit == None and second_e_hit == None

):

return None

insert_size = get_pair_insert_size(ranges, pair)

first_hit = first_s_hit if first_s_hit != None else first_e_hit

second_hit = second_e_hit if second_e_hit != None else second_s_hit

start = genome_interval(

first.pos.chrm,

max(first.pos.start, ranges[first_hit].start),

max(first.pos.start, ranges[first_hit].start),

)

end = genome_interval(

second.pos.chrm,

min(second.pos.end, ranges[second_hit].end),

min(second.pos.end, ranges[second_hit].end),

)

step = plan_step(start, end, "PAIREND")

event_type = get_pair_event_type(pair)

step.info = {"TYPE": event_type, "INSERTSIZE": insert_size}

return insert_size, step

# }}}

# {{{def get_pair_event_type(pe_read):

def get_pair_event_type(pe_read):

"""Decide what type of event the read supports (del/normal, dup, inv)

"""

event_by_strand = {

(True, False): "Deletion/Normal",

(False, True): "Duplication",

(False, False): "Inversion",

(True, True): "Inversion",

}

event_type = event_by_strand[pe_read[0].strand, pe_read[1].strand]

return event_type

# }}}

def jitter(value, bounds: float = 0.1) -> float:

"""

Offset value by a random value within the defined bounds

"""

assert 0.0 <= bounds < 1.0

return value * (1 + bounds * random.uniform(-1, 1))

# {{{def plot_pair_plan(ranges, step, ax):

def plot_pair_plan(ranges, step, ax, marker_size, jitter_bounds):

p = [

map_genome_point_to_range_points(

ranges, step.start_pos.chrm, step.start_pos.start

),

map_genome_point_to_range_points(ranges, step.end_pos.chrm, step.end_pos.end),

]

if None in p:

return False

# some points are far outside of the printable area, so we ignore them

if not points_in_window(p):

return False

READ_TYPES_USED["Paired-end read"] = True

y = step.info["INSERTSIZE"]

# Offset y-values using jitter to avoid overlapping lines

y = jitter(y, bounds=jitter_bounds)

event_type = step.info["TYPE"]

READ_TYPES_USED[event_type] = True

color = COLORS[event_type]

# plot the individual pair

ax.plot(

p,

[y, y],

"-",

color=color,

alpha=0.25,

lw=0.5,

marker="s",

markersize=marker_size,

zorder=10,

)

return True

# }}}

# {{{def plot_pairs(pairs,

def plot_pairs(

pairs, ax, ranges, curr_min_insert_size, curr_max_insert_size, marker_size, jitter_bounds,

):

"""Plots all PairedEnd reads for the region

"""

plan = get_pairs_plan(ranges, pairs)

if not plan:

[curr_min_insert_size, curr_max_insert_size]

max_event, steps = plan

for step in steps:

plot_pair_plan(ranges, step, ax, marker_size, jitter_bounds)

if not curr_min_insert_size or curr_min_insert_size > max_event:

curr_min_insert_size = max_event

if not curr_max_insert_size or curr_max_insert_size < max_event:

curr_max_insert_size = max_event

return [curr_min_insert_size, curr_max_insert_size]

# }}}

##Split Read methods

# {{{class SplitRead:

class SplitRead:

"""container of split read info

Contains start(int), end(int), strand(bool True=forward), query

position (int), MI (int molecular identifier), HP (int haplotype)

"""

def __init__(self, chrm, start, end, strand, query_pos, MI_tag=None, HP_tag=None):

"""Create SplitRead instance

Genomic interval is defined by start, end, and query_pos integers

Strand is opposite of is_reverse

Molecular identifier and Haplotype are integers if present, else

False

"""

self.pos = genome_interval(chrm, start, end)

self.strand = strand

self.query_pos = query_pos

# molecular identifier - linked reads only

self.MI = None

# haplotype - phased reads only

self.HP = 0

if MI_tag:

self.MI = MI_tag

if HP_tag:

self.HP = HP_tag

def __repr__(self):

return "SplitRead(%s,%s,%s,%s,%s,%s,%s)" % (

self.pos.chrm,

self.pos.start,

self.pos.end,

self.strand,

self.query_pos,

self.MI,

self.HP,

)

# }}}

# {{{def calc_query_pos_from_cigar(cigar, strand):

def calc_query_pos_from_cigar(cigar, strand):

"""Uses the CIGAR string to determine the query position of a read

The cigar arg is a string like the following: 86M65S

The strand arg is a boolean, True for forward strand and False for

reverse

Returns pair of ints for query start, end positions

"""

cigar_ops = [[int(op[0]), op[1]] for op in re.findall("(\d+)([A-Za-z])", cigar)]

order_ops = cigar_ops

if not strand: # - strand

order_ops = order_ops[::-1]

qs_pos = 0

qe_pos = 0

q_len = 0

for op_position in range(len(cigar_ops)):

op_len = cigar_ops[op_position][0]

op_type = cigar_ops[op_position][1]

if op_position == 0 and (op_type == "H" or op_type == "S"):

qs_pos += op_len

qe_pos += op_len

q_len += op_len

elif op_type == "H" or op_type == "S":

q_len += op_len

elif op_type == "M" or op_type == "I" or op_type == "X":

qe_pos += op_len

q_len += op_len

return qs_pos, qe_pos

# }}}

# {{{def add_split(read, splits, bam_file, linked_reads):

def add_split(read, splits, bam_file, linked_reads, ignore_hp):

"""adds a (primary, non-supplementary) read to the splits list

Pysam read is added as simpified SplitRead instance to splits

Also added to linked_reads list if there is an associated MI tag

"""

if read.is_secondary:

return

if read.is_supplementary:

return

if not read.has_tag("SA"):

return

qs_pos, qe_pos = calc_query_pos_from_cigar(read.cigarstring, (not read.is_reverse))

HP_tag = False

MI_tag = False

if read.has_tag("MI"):

MI_tag = int(read.get_tag("MI"))

if not ignore_hp and read.has_tag("HP"):

HP_tag = int(read.get_tag("HP"))

sr = SplitRead(

bam_file.get_reference_name(read.reference_id),

read.reference_start,

read.reference_end,

not (read.is_reverse),

qs_pos,

MI_tag,

HP_tag,

)

if sr.MI:

if sr.HP not in linked_reads:

linked_reads[sr.HP] = {}

if sr.MI not in linked_reads[sr.HP]:

linked_reads[sr.HP][sr.MI] = [[], []]

linked_reads[sr.HP][sr.MI][1].append(read.query_name)

if sr.HP not in splits:

splits[sr.HP] = {}

splits[sr.HP][read.query_name] = [sr]

for sa in read.get_tag("SA").split(";"):

if len(sa) == 0:

continue

A = sa.split(",")

chrm = A[0]

pos = int(A[1])

strand = A[2] == "+"

cigar = A[3]

#mapq and nm are never used, annotating this for code readability

mapq = int(A[4])

nm = int(A[5])

qs_pos, qe_pos = calc_query_pos_from_cigar(cigar, strand)

splits[sr.HP][read.query_name].append(

SplitRead(chrm, pos, pos + qe_pos, strand, qs_pos)

)

if len(splits[sr.HP][read.query_name]) == 1:

del splits[sr.HP][read.query_name]

else:

splits[sr.HP][read.query_name].sort(key=lambda x: x.pos.start)

# }}}

# {{{def get_split_plan(ranges, split):

def get_split_plan(ranges, split, linked_plan=False):

"""

There can be 2 or more alignments in a split. Plot only those that are in a

range, and set the insert size to be the largest gap

A split read acts like a long read, so we will covert the split read

to a long read, then convert the long read plan back to a split read plan

"""

alignments = []

for s in split:

# see if they are part of a linked read

if not linked_plan and (s.MI):

return None

alignment = Alignment(s.pos.chrm, s.pos.start, s.pos.end, s.strand, s.query_pos)

alignments.append(alignment)

long_read = LongRead(alignments)

long_reads = {}

long_reads["convert"] = [long_read]

plan = get_long_read_plan("convert", long_reads, ranges)

if not plan:

return None

max_gap, lr_steps = plan

if len(lr_steps) < 3:

return None

sr_steps = []

# a split read will include 3 long read steps, align, event, align

for i in range(0, len(lr_steps), 2):

if i + 2 > len(lr_steps):

break

if (

lr_steps[i].info["TYPE"] == "Align"

and lr_steps[i + 1].info["TYPE"] != "Align"

and lr_steps[i + 2].info["TYPE"] == "Align"

):

start = genome_interval(

lr_steps[i].end_pos.chrm,

lr_steps[i].end_pos.end,

lr_steps[i].end_pos.end,

)

end = genome_interval(

lr_steps[i + 2].start_pos.chrm,

lr_steps[i + 2].start_pos.start,

lr_steps[i + 2].start_pos.start,

)

sr_steps.append(

plan_step(

start,

end,

"SPLITREAD",

info={"TYPE": lr_steps[i + 1].info["TYPE"], "INSERTSIZE": max_gap},

)

)

return max_gap, sr_steps

# }}}

# {{{def get_splits_plan(ranges, splits, linked_plan=False):

def get_splits_plan(ranges, splits, linked_plan=False):

steps = []

max_event = 0

insert_sizes = []

for read_name in splits:

split = splits[read_name]

plan = get_split_plan(ranges, split)

if plan:

insert_size, step = plan

insert_sizes.append(insert_size)

steps += step

if len(insert_sizes) > 0:

max_event = max(insert_sizes)

plan = [max_event, steps]

return plan