[Python] Efficient Mapping of Average Scores Across Interval Windows

I have a df of genomic intervals with millions of rows for example:

chromosome start end
1 300 500
1 400 600
... ... ...

1. find the center of each interval (+/- 250 from the center the start/end)

chromosome start end center
1 300 500 400
1 400 600 500
... ... ...

1. Create new windows (-/+ 250 around the center)

chromosome start end center window_start window_end
1 300 500 400 50 750
1 400 600 500 150 850
... ... ...

1. get the scores for every single position within these windows (using my own function). You must first recreate the df with only Chr, window_start, and window_end.

chromosome start end
1 50 750
1 150 850
... ... ...


My score function will output the following df (much bigger, there will be 500 rows for each interval because 1 row=1 position. For this example, it would be 1000 rows because we have 2 intervals). :

chromosome start end score
1 50 50 .8
1 51 51 .2
1 52 52 .12
...
1 750 750 .43
... ... ...

1. FINAL: I want a plot such that I will be able to have the x axis represent the position relative to the center (-250th position to +250th position) and the y axis will be the average score. So I would take the -250th position across all of my intervals (here I only have 2 intervals) and I would average across those scores.

Example of what I want is here:

I mostly need help with getting from step 4 to 5 in the most efficient way (step 4 isn't the time-consuming part, step 5 is-- according to the way I'm doing my code right now).

This code here works perfectly, but it is way too slow.

def get_postional_avgs(promoter_df):
promoter_df['center'] = (promoter_df['Start'] + promoter_df['End']) // 2

# Calculate window start and end
promoter_df['window_start'] = promoter_df['center'] - 250
promoter_df['window_end'] = promoter_df['center'] + 250

# Prepare a new DataFrame for conservation score requests
conservation_intervals = promoter_df[['Chromosome', 'window_start', 'window_end']]
conservation_intervals.columns = ['Chromosome', 'Start', 'End'] #create a new start/end

# Retrieve base-level conservation scores
base_level_scores = get_base_scores(conservation_intervals)

# Initialize a DataFrame to hold average scores
average_scores = []

for index, row in promoter_df.iterrows():
center = row['center']

# Extract scores for the current cCRE
scores_for_cCRE = base_level_scores[
(base_level_scores['Chromosome'] == row['Chromosome']) &
(base_level_scores['Start'] >= (center - 250)) &
(base_level_scores['End'] <= (center + 250))
]

# Initialize a list to store average scores for each position
avg_scores = {}

for pos in range(-250, 251):
current_position = center + pos
current_scores = scores_for_cCRE[scores_for_cCRE['Start'] == current_position]

if not current_scores.empty:
avg_scores[pos] = current_scores['Score'].mean()
else:
avg_scores[pos] = np.nan # Handle missing scores

# Append to the results list
average_scores.append(avg_scores)

# Convert average scores to a DataFrame
avg_score_df = pd.DataFrame(average_scores)



If you want to see the code for the get_score():

import pyBigWig

def get_base_level_conservation_scores(df):
"""
Retrieve all conservation scores that occur within a given genomic interval.
The output DataFrame will be expanded to include scores for each base.
Output DataFrame has 1-indexed coordinates.
"""
df = df.reset_index(drop=True)

bw = pyBigWig.open("hg38.phyloP100way.bw")

# Convert to 0-indexing (assuming input is 1-indexed)
df['Start'] = df['Start'] - 1

# Initialize score data list
score_data = {
'Chromosome': [],
'Start': [],
'End': [],
'Score': []
}

# Process each interval
for i in range(len(df)):
chrom = df['Chromosome']
start = df['Start']
end = df['End']

# Check for invalid intervals
if start < 0 or end <= 0 or start >= end:
print(f"Invalid interval for {chrom}: start={start}, end={end}")
continue # Skip this iteration

# Get scores for the interval
try:
scores_per_base = bw.intervals(chrom, start, end)
except RuntimeError as e:
print(f"Error for {chrom}: start={start}, end={end} - {e}")
continue # Skip this iteration if there's an error

if scores_per_base is None or len(scores_per_base) == 0:
# If there are no scores in the specified interval, add NaN entries
for position in range(start, end):
score_data['Chromosome'].append(chrom)
score_data['Start'].append(position)
score_data['End'].append(position)
score_data['Score'].append(np.nan)
else:
for start_pos, end_pos, score in scores_per_base:
for position in range(start_pos, end_pos):
score_data['Chromosome'].append(chrom)
score_data['Start'].append(position)
score_data['End'].append(position)
score_data['Score'].append(score)

# Convert to DataFrame, drop duplicates, and drop NaN scores
score_data_df = pd.DataFrame(score_data)
score_data_df = score_data_df.drop_duplicates(subset=['Chromosome', 'Start', 'End'])
score_data_df = score_data_df.dropna(subset=['Score'])

# Close the bigWig file
bw.close()

return score_data_df

Continue reading...