Skip to content
Jump to bottom

Peak calling to dev after resolving overlaps #41

Merged
merged 4 commits into from
Jan 8, 2019
Merged

Peak calling to dev after resolving overlaps #41

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
aa1e9f0
fixing the bug with footprint length
anastasiia Jan 7, 2019
c452934
fixing the bug with overlaps + adding a merging function if the footp…
anastasiia Jan 7, 2019
7a5d6df
allowing user to set the parameter for max bp allowed in between the …
anastasiia Jan 7, 2019
0308a0c
changing the name of footprints from footprint_ to f_ to make sure th…
anastasiia Jan 8, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
67 changes: 62 additions & 5 deletions bin/1.1_footprint_extraction/footprints_extraction.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ def parse_args():
parser.add_argument('--window_length', default=200, type=int, help='Please enter the length for a window, by defauld 200 bp.')
parser.add_argument('--step', default=100, type=int, help='Please enter a step to move the window, by default 100 bp.')
parser.add_argument('--percentage', default=0, type=int, help='Please enter a percentage to be added to background while searching for footprints, by default 0%%.')
parser.add_argument('--max_bp_between', default=6, type=int, help='Please enter the number of bp allowed to be in between two footprints, by default 6 bp.')
parser.add_argument('--silent', action='store_true', help='While working with data write the information only into ./footprints_extraction.log.')
args = parser.parse_args()

Expand Down Expand Up @@ -154,6 +155,7 @@ def save_footprint(footprint_count, footprint_scores, peak_footprints, chromosom
old_start = peak_footprints[existing_footprint_name]['start']
old_end = peak_footprints[existing_footprint_name]['end']
old_score = peak_footprints[existing_footprint_name]['score']
old_length = peak_footprints[existing_footprint_name]['len']

if footprint_start >= old_start and footprint_start <= old_end: #the start of the new footprint is between the start and end of an old footprint
if footprint_end > old_end: #the new footprint is not completely inside the old one
Expand All @@ -163,19 +165,21 @@ def save_footprint(footprint_count, footprint_scores, peak_footprints, chromosom

peak_footprints[existing_footprint_name]['end'] = footprint_end
peak_footprints[existing_footprint_name]['score'] = footprint_score
peak_footprints[existing_footprint_name]['len'] = footprint_end - old_start
#we can not update the max_pos as we do not have the information about scores array of the existing footprint
#else: the new footprint is completely inside the old one, do nothing
save_current_footprint = False
break

elif footprint_end >= old_start and footprint_end <= old_end: #the end of the new footprint is between the start and end of an old footprint
if footprint_start < old_start: #the new footprint is not completely inside the old one
if footprint_start < old_start: #the new footprint is not completely inside the old one
#update the information about the existing footprint
#find the average of both scores
footprint_score = (peak_footprints[existing_footprint_name]['score'] + footprint_score) / 2

peak_footprints[existing_footprint_name]['start'] = footprint_start
peak_footprints[existing_footprint_name]['score'] = footprint_score
peak_footprints[existing_footprint_name]['len'] = old_end - footprint_start
#else do nothing
save_current_footprint = False
break
Expand All @@ -191,7 +195,7 @@ def save_footprint(footprint_count, footprint_scores, peak_footprints, chromosom

if save_current_footprint == True:
#make sure to save this footprint
footprint_name = "footprint_" + str(footprint_count)
footprint_name = "f_" + str(footprint_count)

peak_footprints[footprint_name] = peak_footprints.get(footprint_name, {})
peak_footprints[footprint_name] = {'chromosom': chromosom, 'start': footprint_start, 'end': footprint_end, 'score': footprint_score, 'len': len(footprint_scores), 'bonus': bonus_info_from_bed, 'max_pos': max_pos}
Expand Down Expand Up @@ -270,11 +274,62 @@ def search_in_window(peak_footprints, footprint_count, chromosom, peak_start, pe

return peak_footprints, footprint_count

#this function checks if there are footprints with a small amount of bp in between, and if so, merges them together
#if at this point there are still overlaps of footprints, this function will delete them
#the input parameter are: dictionary with footprints within one peak, and the max number of bp allowed to be in between the footprints
#the output is the renewed dictionary containing only the best footprints for the output file
def check_and_merge(peak_footprints, max_bp_between):
peak_footprints_new = {}

for footprint_to_check in peak_footprints.keys():
start_to_check = peak_footprints[footprint_to_check]['start']
end_to_check = peak_footprints[footprint_to_check]['end']

merge_footprints_left = None

for compared_footprint in peak_footprints.keys():

if start_to_check > peak_footprints[compared_footprint]['start'] and start_to_check - peak_footprints[compared_footprint]['end'] < max_bp_between:
#make compared_footprint longer
merge_footprints_left = False
break
elif end_to_check < peak_footprints[compared_footprint]['end'] and peak_footprints[compared_footprint]['start'] - end_to_check < max_bp_between:
#make footprint_to_check longer
merge_footprints_left = True
break

if merge_footprints_left: #if the merging left is enabled
#check if this footprint can be merged with the compared_footprint
#if compared footprint is not in peak_footprint_new.keys(), the next loop will check for this footprint. There is no need for doulbe check now
if start_to_check < peak_footprints[compared_footprint]['start'] and compared_footprint in peak_footprints_new.keys():
#update the start position
peak_footprints_new[compared_footprint]['start'] = start_to_check
#update the length
peak_footprints_new[compared_footprint]['len'] = peak_footprints[compared_footprint]['end'] - start_to_check
#update the score
peak_footprints_new[compared_footprint]['score'] = (peak_footprints[footprint_to_check]['score'] + peak_footprints[compared_footprint]['score']) / 2

elif merge_footprints_left == False: #otherwise merge right
#check if the merging is possible
if end_to_check > peak_footprints[compared_footprint]['end'] and compared_footprint in peak_footprints_new.keys():
#update the end position
peak_footprints_new[compared_footprint]['end'] = end_to_check
#update the length
peak_footprints_new[compared_footprint]['len'] = end_to_check - peak_footprints[compared_footprint]['end']
#update the score
peak_footprints_new[compared_footprint]['score'] = (peak_footprints[footprint_to_check]['score'] + peak_footprints[compared_footprint]['score']) / 2

else: #save the current footprint, as it should not be merged
peak_footprints_new[footprint_to_check] = peak_footprints_new.get(footprint_to_check, [])
peak_footprints_new[footprint_to_check] = peak_footprints[footprint_to_check]

return peak_footprints_new

#this function uses the information provided from the .bed file to look for footprints within the peaks of interest
#as input the information from the original bed file, as well as bigwig file is needed
#the optional parameters window_length, step and percentage are needed as well to use the sliding window algorithm and work with the "background" score
#the output of this function is a dictionary contains all the found footprints ready to write to the output file
def find_peaks_from_bw(bed_dictionary, bw_file, window_length, step, percentage):
def find_peaks_from_bw(bed_dictionary, bw_file, window_length, step, percentage, max_bp_between):

logger.info('Looking for footprints within peaks...')

Expand All @@ -295,9 +350,11 @@ def find_peaks_from_bw(bed_dictionary, bw_file, window_length, step, percentage)
peak_end = int(positions[1])

scores_in_peak = np.nan_to_num(np.array(list(bw_open.values(chromosom, peak_start, peak_end)))) #save the scores to an array

peak_footprints, footprint_count = search_in_window(peak_footprints, footprint_count, chromosom, peak_start, peak_end, scores_in_peak, window_length, bed_dictionary[header], step, percentage)

#double check for overlaps and possibly merging of footprints having up to 5 bp in between
peak_footprints = check_and_merge(peak_footprints, max_bp_between)

for footprint_name in peak_footprints.keys():
all_footprints[footprint_name] = all_footprints.get(footprint_name, {})
all_footprints[footprint_name] = peak_footprints[footprint_name]
Expand Down Expand Up @@ -381,7 +438,7 @@ def main():
logger.info("The script footprints_extraction.py was called using these parameters: " + str(vars(args)))

bed_dictionary = make_bed_dictionary(args.bed)
all_footprints = find_peaks_from_bw(bed_dictionary, args.bigwig, args.window_length, args.step, args.percentage)
all_footprints = find_peaks_from_bw(bed_dictionary, args.bigwig, args.window_length, args.step, args.percentage, args.max_bp_between)
write_to_bed_file(all_footprints, args.output_file)

logger.info("the number of peaks: " + str(len(bed_dictionary)))
Expand Down