From cf73458d6b93945e5ae22bbb4eb7377538e8ada3 Mon Sep 17 00:00:00 2001 From: Anastasia Gorbarenko <131146976+CaptnClementine@users.noreply.github.com> Date: Wed, 8 May 2024 12:35:39 -0700 Subject: [PATCH 1/2] Add new version of flexible split-merge detection --- src/tads_find_overlaps.py | 148 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 148 insertions(+) create mode 100644 src/tads_find_overlaps.py diff --git a/src/tads_find_overlaps.py b/src/tads_find_overlaps.py new file mode 100644 index 0000000..109b438 --- /dev/null +++ b/src/tads_find_overlaps.py @@ -0,0 +1,148 @@ +import os + +import numpy as np +import pandas as pd +from concurrent.futures import ThreadPoolExecutor + + +def find_overlaps(table1: pd.DataFrame, table2: pd.DataFrame, bins: int = 100_000) -> pd.DataFrame: + """ + Finds overlaps between two tables of TADs with a tolerance of 1.5 bins. + + Args: + table1 (pd.DataFrame): DataFrame representing the first table of TADs. Must include specified columns: + - 'chrom' with chromosome number (int ot str) + - 'start' - coordinate of TAD start + - 'end' - coordinate of TAD end + table2 (pd.DataFrame): DataFrame representing the second table of TADs. Must include same specified columns + bins (int, optional): Size of bins. Defaults to 100'000. + + Returns: + pd.DataFrame: DataFrame containing overlaps between the two tables. Columns include: + - Region1 - chromosome number from first table + - Region2 - chromosome number from second table + - Table1_length - length for each TAD from first table + - Table2_length - length for each TAD from second table + - Table1_coordinates - coordinates for each TAD from first table + - Table2_coordinates - coordinates for each TAD from second table + Example: + results = find_overlaps(table1, table2) + """ + bins_tolerance = bins*1.5 + overlaps = [] + for idx1, row1 in table1.iterrows(): + for idx2, row2 in table2.iterrows(): + if row1['chrom'] != row2['chrom']: #проверка по хромосоме + continue + if ( + (abs(row1['start'] - row2['start']) <= bins_tolerance) and (row2['end'] < row1['end']) #первый микротад + ) or ( + (row1['start'] < row2['start']) and (row2['end'] < row1['end']) #средний микротад + ) or ( + (row1['start'] < row2['start']) and (abs(row1['end'] - row2['end']) <= bins_tolerance) #последний микротад + ): + overlaps.append((row1['chrom'], row2['chrom'], + abs(row1['end'] - row1['start']), + abs(row2['end'] - row2['start']), + (row1['start'], row1['end']), #tuple или pandas.interval (?) + (row2['start'], row2['end']))) + return pd.DataFrame(overlaps, columns=['Region1', 'Region2', 'Table1_length', 'Table2_length', 'Table1_coordinates', 'Table2_coordinates']) + + +def filter_results_diff_split(results: pd.DataFrame, epsilon: int = 200_000) -> pd.DataFrame: + """ + Filters regions with different TAD sizes between the same Region1 values. + + Args: + results (pd.DataFrame): Output DataFrame from 'find_overlaps' fuction. Must include specified columns: + - Region1 - chromosome number from first table + - Region2 - chromosome number from second table + - Table1_length - length for each TAD from first table + - Table2_length - length for each TAD from second table + - Table1_coordinates - coordinates for each TAD from first table + - Table2_coordinates - coordinates for each TAD from second table + epsilon (int, optional): Threshold for considering TAD size differences significant. Defaults to 200000. + + Returns: + pd.DataFrame: DataFrame containing filtered regions with significant TAD size differences. + Example: + filtered_results = filter_results_diff_split(results) + """ + filtered_regions = [] + for region1 in results['Region1'].unique(): + df_region1 = results[results['Region1'] == region1] + unique_lengths = df_region1['Table2_length'].unique() + + if len(unique_lengths) == 1: #если одно значение - выкидываем его из результатов + continue + + for length in unique_lengths: + min_length = np.min(df_region1['Table2_length']) + max_length = np.max(df_region1['Table2_length']) + if max_length - min_length > epsilon: # если разница между двумя тадами больше, чем epsilon - то это наши разноразмерные тады + filtered_regions.extend(df_region1.to_dict('records')) + break + filtered_regions = pd.DataFrame(filtered_regions) + #следующая строка - для проверки странных ТАДов (может они и есть шифты, а может для обратной таблицы это сплит) + index_not_unique_TAD = filtered_regions.groupby(by=['Region1', 'Table1_coordinates'])['Region2'].transform('count') != 1 + filtered_regions = filtered_regions.loc[index_not_unique_TAD] + return filtered_regions + + +def make_correct_names(result_table: pd.DataFrame) -> pd.DataFrame: + """ + Rename columns and extract start and end coordinates from 'result_table'. + + Parameters: + result_table (pd.DataFrame): DataFrame obtained from 'filter_results_diff_split' function. + + Returns: + pd.DataFrame: DataFrame with corrected column names and extracted coordinates. + """ + result_table = result_table.rename(columns={'Region1': 'chrom'}) + result_table['start_1'], result_table['end_1'] = zip(*result_table['Table1_coordinates']) + result_table['start_2'], result_table['end_2'] = zip(*result_table['Table2_coordinates']) + result_table = result_table.drop(columns=['Region2', 'Table1_length', 'Table2_length', + 'Table1_coordinates', 'Table2_coordinates']) + return result_table + + +def main_flex_split_merge(table1: pd.DataFrame, table2: pd.DataFrame, save_directory: os.path = './') -> pd.DataFrame: + """ + Merge and split overlapping regions with flexible size from two input tables and create a final DataFrame. + + Parameters: + table1 (pd.DataFrame): First input DataFrame containing regions. + table2 (pd.DataFrame): Second input DataFrame containing regions. + + Returns: + pd.DataFrame: Final DataFrame with merged and split regions. + """ + with ThreadPoolExecutor() as executor: + # Submit tasks for parallel execution + future_overlaps_table1 = executor.submit(find_overlaps, table1, table2) + future_overlaps_table2 = executor.submit(find_overlaps, table2, table1) + + # Wait for tasks to complete and get results + overlaps_table1 = future_overlaps_table1.result() + overlaps_table2 = future_overlaps_table2.result() + + with ThreadPoolExecutor() as executor: + # Submit tasks for parallel execution + future_result_split = executor.submit(filter_results_diff_split, overlaps_table1) + future_result_merge = executor.submit(filter_results_diff_split, overlaps_table2) + + # Wait for tasks to complete and get results + result_table_split = future_result_split.result() + result_table_merge = future_result_merge.result() + + result_table_split = make_correct_names(result_table_split) + result_table_merge = make_correct_names(result_table_merge) + result_table_merge = result_table_merge.rename(columns = {'start_1':'start_2','end_1':'end_2','start_2':'start_1','end_2':'end_1'}) + + result_table_split.to_csv(f'{save_directory}/flex_split.csv') + result_table_merge.to_csv(f'{save_directory}/flex_merge.csv') + + + + From e4b40274835ffd59a079778b6231d9e30280bb1f Mon Sep 17 00:00:00 2001 From: Anastasia Gorbarenko <131146976+CaptnClementine@users.noreply.github.com> Date: Wed, 8 May 2024 22:40:14 +0300 Subject: [PATCH 2/2] Add fles split merge function --- CTADO.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/CTADO.py b/CTADO.py index e965b41..197696f 100644 --- a/CTADO.py +++ b/CTADO.py @@ -5,10 +5,13 @@ from src.calculate_intensity_change import count_tads_change_intensity from src.tads_plot import visualisation from src.split_merge_detect import main_split_merge_detection +from src.tads_find_overlaps import main_flex_split_merge INTENSITY = 'intensity_change_result.csv' SPLIT = 'split_coords.csv' MERGE = 'merge_coords.csv' +FSPLIT = 'flex_split.csv' +FMERGE = 'flex_merge.csv' if __name__ == "__main__": @@ -49,8 +52,11 @@ f'{args.save_directory}/{args.clr1_filename}_{args.window}_result_df.csv', f'{args.save_directory}/{args.clr2_filename}_{args.window}_result_df.csv', args.save_directory) + main_flex_split_merge(f'{args.save_directory}/{args.clr1_filename}_{args.window}_result_df.csv', + f'{args.save_directory}/{args.clr2_filename}_{args.window}_result_df.csv', + args.save_directory) sys.stdout.write(f'Visualising...\n'); sys.stdout.flush() - for file in [INTENSITY, SPLIT, MERGE]: + for file in [INTENSITY, SPLIT, MERGE, FSPLIT, FMERGE]: type_of_change = file[:file.find('_')] visualisation(args.clr1_filename, args.clr2_filename, args.clr1_boundaries_name, args.clr2_boundaries_name, args.resolution, args.binsize, args.window, f'{args.save_directory}/{file}', type_of_change, args.save_directory) - sys.stdout.write(f'CTADO completed successfully! Output location:\n{os.path.abspath(args.save_directory)}\n'); sys.stdout.flush() \ No newline at end of file + sys.stdout.write(f'CTADO completed successfully! Output location:\n{os.path.abspath(args.save_directory)}\n'); sys.stdout.flush()