'Deletion/Normal': 'black',

'Deletion': 'black',

'Duplication': 'red',

'Inversion': 'blue',

'InterChrmInversion': 'blue',

'InterChrm': 'black'

"Deletion/Normal":False,

"Duplication":False,

"Inversion":False,

"Aligned long read": False,

"Linked read": False,

"Split-read": False,

"Paired-end read": False

'M' : 0,

'I' : 1,

'D' : 2,

'N' : 3,

'S' : 4,

'H' : 5,

'P' : 6,

'=' : 7,

'X' : 8,

'B' : 9

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):

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 gi.chrm < self.chrm or gi.end < self.start:

return -1

elif gi.chrm > self.chrm or gi.start > self.end:

return 1

else:

for j in range(len(ranges)):

r = ranges[j]

if r.chrm.strip("chr") == chrm.strip("chr") and r.start <= point and r.end >= point:

return j

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))

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

if None in points or points[0] < -5 or points[1] < -5 or points[0] > 5 or points[1] > 5:

return False

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

tbx = pysam.TabixFile(datafile)

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:

chrom = 'chr' + chrm

try:

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

except ValueError:

sys.exit('Warning: Could not fetch ' + \

chrm + ':' + start + '-' + end + \

' from ' + datafile)

"""Adds a read to the known coverage

chrm = bam_file.get_reference_name(read.reference_id)

hp = 0

if read.has_tag('HP'):

hp = int(read.get_tag('HP'))

if hp not in coverage:

coverage[hp] = {}

if chrm not in coverage[hp]:

coverage[hp][chrm] = {}

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']]:

for pos in range(curr_pos, curr_pos+length):

if pos not in coverage[hp][chrm]:

coverage[hp][chrm][pos] = [0,0]

#the two coverage tracks are [0] high-quality and [1]

#low-quality

if minq and (read.mapping_quality > minq):

coverage[hp][chrm][pos][0] += 1

else:

coverage[hp][chrm][pos][1] += 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:

ax,

ranges,

hp_count,

max_coverage,

tracktype,

yaxis_label_fontsize,

same_yaxis_labels=False):

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

cover_x = []

cover_y_lowqual = []

cover_y_highqual = []

cover_y_all = []

for i in range(len(ranges)):

r = ranges[i]

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

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 and same_yaxis_labels:

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 = max(cover_y_all.max(), cover_y_all.max())

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=.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=.15)

ax2.fill_between(cover_x,

cover_y_highqual,

cover_y_lowqual,

color='darkgrey',

step="pre",

alpha=.4)

ax2.fill_between(cover_x,

cover_y_lowqual,

bottom_fill,

color='grey',

step="pre",

alpha=0.15)

#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.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)

"""container of paired-end read info

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

"""Create PairedEnd instance

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,

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

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 read.has_tag('HP'):

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

READ_TYPES_USED["Paired-end read"] = True

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:

READ_TYPES_USED["Linked read"] = True

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 )

"""Downsamples paired-end reads

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 = [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(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]

"""Extracts the integer insert sizes for all pairs

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)

''' 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:

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]

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}

"""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]

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

y = step.info['INSERTSIZE']

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)

ax,

ranges,

curr_min_insert_size,

curr_max_insert_size):

"""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)

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

"""container of split read info

def __init__(self,

chrm,

start,

end,

strand,

query_pos,

MI_tag=None,

HP_tag=None):

"""Create SplitRead instance

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,

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

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

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

if read.is_secondary: return

if read.is_supplementary: return

if not read.has_tag('SA'): return

READ_TYPES_USED["Split-read"] = True

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 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_read[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 = 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:

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

first = split[0]

second = split[1]

if first.pos.start > second.pos.end or first.pos.chrm > second.pos.chrm:

second = split[0]

first = split[1]

# first.strand, second.strand,

# first.query<second.query,first.start<second.start

event_type_by_strand_and_order = {

(True, False) : 'Inversion', #mixed strands

(False, True) : 'Inversion', #mixed strands

(True, True, True) : 'Deletion/Normal', #forward strand

(True, True, False) : 'Duplication', #forward strand

(False, False, False, False): 'Deletion/Normal', #reverse strand

(False, False, False, True) : 'Duplication', #reverse strand

(False, False, True, True) : 'Deletion/Normal', #reverse strand

(False, False, True, False) : 'Duplication' #reverse strand

}

orientations = [first.strand, second.strand]

#if same strand, need query position info

if orientations[0] == orientations[1]:

#first query position smaller than second query position,

#normal for forward strand

orientations.append(first.query_pos < second.query_pos)

#reverse strand requires start position info

if False in orientations[:2]:

#first start smaller than second start, normal for forward strand

orientations.append(first.pos.start < second.pos.start)

event_type = event_type_by_strand_and_order[tuple(orientations)]

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

if event_type == 'Inversion':

event_type = 'InterChrmInversion'

else:

event_type = 'InterChrm'

"""Extracts the integer gap sizes for all split reads

split_insert_sizes = []

for hp in splits:

for read_name in splits[hp]:

sizes = get_split_insert_sizes(ranges, splits[hp][read_name])

if sizes:

split_insert_sizes += sizes

split_insert_sizes = []

srs = sorted(srs, key = lambda x: x.query_pos)

#only keep the alignments that are in a range

srs = [sr for sr in srs if \

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

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

if len(srs) < 2:

return None

last = srs[0]

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

curr = srs[i]

#INTERCHROM

if last.pos.chrm != curr.pos.chrm:

split_insert_sizes.append(INTERCHROM_YAXIS)

#Inversion

elif last.strand != curr.strand:

split_insert_sizes.append(abs(curr.pos.end - last.pos.end))

#Duplication

elif curr.pos.start < last.pos.end:

split_insert_sizes.append(\

abs(last.pos.end - curr.pos.start))

#Deletion

elif curr.pos.start > last.pos.end:

split_insert_sizes.append(\

abs(curr.pos.start - last.pos.end))

else:

sys.stderr.write('WARNING: Could not classify event:' \

+ str(last) + str(curr))

last = curr

'''

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}))

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]