Skip to content

Commit

Permalink
fix bug in graph prunning
Browse files Browse the repository at this point in the history
  • Loading branch information
@iprada
iprada committed Mar 13, 2020
1 parent fd191d6 commit 8b80422
Show file tree
Hide file tree
Showing 7 changed files with 59 additions and 41 deletions.
23 changes: 12 additions & 11 deletions circlemap/bam2bam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,17 +95,22 @@ def __init__(self, input_bam,output,qname_bam,genome_fasta,directory,mapq_cutoff

def listener_writer(self,bam):

f = open('test.sam',"w")

while True:
try:
# Read from the queue and do nothing
read = self.queue.get()


if read == "DONE":
f.close()
break
else:
pysam_read = ps.AlignedSegment.fromstring(read,bam.header)
f.write(read + "\n")
bam.write(pysam_read)

except BaseException as e:
print("Error on the writer process. Shutting down")
print(e)
Expand Down Expand Up @@ -169,7 +174,7 @@ def realign(self,peaks):



if len(candidate_mates) > 0:
if len(candidate_mates) > 0 or candidate_mates != None:


realignment_interval_extended = get_realignment_intervals(candidate_mates,extension,self.interval_p,
Expand Down Expand Up @@ -319,17 +324,13 @@ def realign(self,peaks):


except BaseException as e:
traceback.print_exc(file=sys.stdout)
warnings.warn(
"Failed on interval %s due to the error %s" % (
str(interval), str(e)))
return([1,1])


if self.verbose < 2:
print(interval)
traceback.print_exc(file=sys.stdout)
warnings.warn(
"Failed on interval %s due to the error %s" % (
str(interval), str(e)))

else:
pass


ecc_dna.close()
Expand All @@ -338,7 +339,7 @@ def realign(self,peaks):
except:
print("Failed on cluster:")
print(traceback.print_exc(file=sys.stdout))
return([0,0])
return([1,1])

genome_fa.close()
ecc_dna.close()
Expand Down
Empty file added circlemap/call.py
View file
Empty file.
18 changes: 18 additions & 0 deletions circlemap/circle_map.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@
import glob
from tqdm import *
from circlemap.__version__ import __version__ as cm_version
import datetime

class circle_map:

Expand Down Expand Up @@ -166,10 +167,19 @@ def __init__(self):

pool = mp.Pool(processes=self.args.threads)


#progress bar
with tqdm(total=len(splitted)) as pbar:
for i,exits in tqdm(enumerate(pool.imap_unordered(object.realign, splitted))):
pbar.update()
#kill if process returns 1,1
if exits == [1,1]:
pool.close()
pool.terminate()
pbar.close()
print(datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S:"),
"An error happenend during execution. Exiting")
sys.exit()

pbar.close()
pool.close()
Expand Down Expand Up @@ -247,6 +257,14 @@ def __init__(self):
with tqdm(total=len(splitted)) as pbar:
for i,exits in tqdm(enumerate(pool.imap_unordered(object.realign,splitted))):
pbar.update()
# kill if process returns 1,1
if exits == [1, 1]:
pool.close()
pool.terminate()
pbar.close()
print(datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S:"),
"An error happenend during execution. Exiting")
sys.exit()

pbar.close()
pool.close()
Expand Down
24 changes: 15 additions & 9 deletions circlemap/extract_circle_SV_reads.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,7 @@ def extract_sv_circleReads(self):

# gets its on mapq since mate is unmapped
if read1.mapq >= self.mapq_cutoff:

read1.tags += [('MQ', read1.mapq)]
circle_sv_reads.write(read1)

Expand All @@ -185,6 +186,7 @@ def extract_sv_circleReads(self):

# gets its on mapq since mate is unmapped
if read2.mapq >= self.mapq_cutoff:

read2.tags += [('MQ', read2.mapq)]
circle_sv_reads.write(read2)

Expand All @@ -197,7 +199,7 @@ def extract_sv_circleReads(self):
# gets its on mapq since mate is unmapped
if read2.mapq >= self.mapq_cutoff:
read2.tags += [('MQ', read2.mapq)]
circle_sv_reads.write(read1)
circle_sv_reads.write(read2)

else:

Expand Down Expand Up @@ -338,6 +340,7 @@ def extract_sv_circleReads(self):
if self.no_hard_clipped == False:

if read1.mapq >= self.mapq_cutoff:

read1.tags += [('MQ', read2.mapq)]
circle_sv_reads.write(read1)

Expand All @@ -346,25 +349,26 @@ def extract_sv_circleReads(self):
pass


if is_soft_clipped(read2) == True:
if is_soft_clipped(read2) == True:


if self.no_soft_clipped == False:
if self.no_soft_clipped == False:

if read2.mapq >= self.mapq_cutoff:
if read2.mapq >= self.mapq_cutoff:

read2.tags += [('MQ', read1.mapq)]
circle_sv_reads.write(read2)
read2.tags += [('MQ', read1.mapq)]
circle_sv_reads.write(read2)


else:
pass
else:
pass
else:
if read1.is_unmapped == False:
if is_soft_clipped(read1) == True:
if self.no_soft_clipped == False:

if read1.mapq >= self.mapq_cutoff:

read1.tags += [('MQ', read2.mapq)]
circle_sv_reads.write(read1)

Expand All @@ -380,6 +384,7 @@ def extract_sv_circleReads(self):

#gets its on mapq since mate is unmapped
if read1.mapq >= self.mapq_cutoff:

read1.tags += [('MQ', read1.mapq)]
circle_sv_reads.write(read1)

Expand All @@ -393,6 +398,7 @@ def extract_sv_circleReads(self):

#gets its on mapq since mate is unmapped
if read2.mapq >= self.mapq_cutoff:

read2.tags += [('MQ', read2.mapq)]
circle_sv_reads.write(read2)

Expand All @@ -405,7 +411,7 @@ def extract_sv_circleReads(self):
# gets its on mapq since mate is unmapped
if read2.mapq >= self.mapq_cutoff:
read2.tags += [('MQ', read2.mapq)]
circle_sv_reads.write(read1)
circle_sv_reads.write(read2)

else:

Expand Down
20 changes: 7 additions & 13 deletions circlemap/realigner.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ def realign(self,peaks):



if len(candidate_mates) > 0:
if len(candidate_mates) > 0 or candidate_mates != None:


realignment_interval_extended = get_realignment_intervals(candidate_mates,extension,self.interval_p,
Expand All @@ -202,7 +202,6 @@ def realign(self,peaks):
continue



iteration_results = []
iteration_discordants = []
disorcordants_per_it = 0
Expand Down Expand Up @@ -386,16 +385,11 @@ def realign(self,peaks):

except BaseException as e:


if self.verbose < 2:
print(interval)
traceback.print_exc(file=sys.stdout)
warnings.warn(
"Failed on interval %s due to the error %s" % (
str(interval), str(e)))

else:
pass
traceback.print_exc(file=sys.stdout)
warnings.warn(
"Failed on interval %s due to the error %s" % (
str(interval), str(e)))
return([1,1])


ecc_dna.close()
Expand All @@ -413,7 +407,7 @@ def realign(self,peaks):
except:
print("Failed on cluster:")
print(traceback.print_exc(file=sys.stdout))
return([0,0])
return([1,1])

sorted_bam.close()
genome_fa.close()
Expand Down
15 changes: 7 additions & 8 deletions circlemap/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -517,9 +517,7 @@ def get_mate_intervals(sorted_bam,interval,mapq_cutoff,verbose,only_discordants)

except BaseException as e:

if verbose < 2:

warnings.warn(
warnings.warn(
"WARNING: Could not get mate interval priors for the interval %s due to the following error %s \n Skipping interval" % (str(interval),str(e)))


Expand Down Expand Up @@ -639,7 +637,7 @@ def get_realignment_intervals(bed_prior,interval_extension,interval_p_cutoff,ver
extended = []




#if argmax is turn on interval_p is 0
if interval_p_cutoff == 0:
Expand All @@ -648,6 +646,7 @@ def get_realignment_intervals(bed_prior,interval_extension,interval_p_cutoff,ver
candidate_mates = candidate_mates.loc[candidate_mates['probability'] == candidate_mates['probability'].max()]

for item,row in candidate_mates.iterrows():

if ('LR' in orientation) or ('L' and 'R' in orientation):


Expand Down Expand Up @@ -685,7 +684,8 @@ def get_realignment_intervals(bed_prior,interval_extension,interval_p_cutoff,ver

#small pseudocount to denominator to avoid div by zero

if interval['probability']/(sum+0.00000001) >= interval_p_cutoff:

if interval['probability'] >= interval_p_cutoff:

if ('LR' in orientation) or ('L' and 'R' in orientation):

Expand Down Expand Up @@ -722,9 +722,8 @@ def get_realignment_intervals(bed_prior,interval_extension,interval_p_cutoff,ver

except BaseException as e:

if verbose < 2:

warnings.warn(
warnings.warn(
"WARNING: Could not compute the probability for the mate interval priors %s due to the following error %s \n Skipping intervals" % (
str(bed_prior), str(e)))

Expand Down Expand Up @@ -1467,7 +1466,7 @@ def assign_discordants(split_bed,discordant_bed,insert_mean,insert_std):
"""Function that takes as input the the discordant reads supporting an interval and assigns them to the
interval if they are close by (using the insert size estimate)"""

max_dist = (insert_mean / 2) + (3 * insert_std)
max_dist = (insert_mean) + (5 * insert_std)

splits = pd.DataFrame.from_records(split_bed, columns=['chrom', 'start', 'end', 'read', 'iteration',
'score']).sort_values(['chrom', 'start', 'end'])
Expand Down
Empty file added tests/run_call.py
View file
Empty file.

0 comments on commit 8b80422

Please sign in to comment.