diff --git a/mumdia.py b/mumdia.py index 5d8cb4a..8f7dc6e 100644 --- a/mumdia.py +++ b/mumdia.py @@ -63,6 +63,8 @@ def __getattr__(self, name): from prediction_wrappers.wrapper_ms2pip import ( get_predictions_fragment_intensity_main_loop, ) +from quantification.lfq import quantify_fragments +from utilities.plotting import plot_XIC_with_margins, plot_rt_margin_histogram from utilities.logger import log_info # Re-export for backward compatibility @@ -283,7 +285,7 @@ def run_mokapot(output_dir="results/") -> None: f"mokapot is not installed or failed to import ({e}). Skipping mokapot run." ) return None - psms = mokapot.read_pin(f"{output_dir}outfile.pin") + psms = mokapot.read_pin(f"{output_dir}/outfile.pin") model = KerasClassifier( build_fn=create_model, epochs=100, batch_size=1000, verbose=10 @@ -911,6 +913,233 @@ def extract_intensities(scannr, charge, calcmass): return df_psms +def calculate_rt_margins_intensity_based(df_fragments: pl.DataFrame, intensity_threshold: float, output_dir='xics') -> pl.DataFrame: + """ + Calculate retention time margins based on a relative intensity threshold of the apex intensity fragment. + The margins are determined by finding the retention times where the fragment intensity + drops below the specified fraction of the apex intensity on both sides of the apex. + If the intensity never drops below the threshold on one side, the margin is set to the + first/last retention time where the most intense fragment was detected. + The function also generates and saves a plot of the XIC with the calculated margins. + + Parameters + ---------- + df_fragments : pl.DataFrame + DataFrame containing fragment ion information for a single peptidoform. + intensity_threshold : float + Intensity threshold (as a fraction of apex intensity) to define retention time margins. + output_dir : str + Directory to save the XIC plots with margins. + Returns + ------- + left_bound : float + Left retention time margin. + right_bound : float + Right retention time margin. + apex_rt : float + Retention time at apex intensity. + """ + + # Sort by rt + df_sorted = df_fragments.sort("rt") + # Find apex + apex_idx = df_sorted["fragment_intensity"].arg_max() + apex_rt = df_sorted["rt"][apex_idx] + apex_intensity = df_sorted["fragment_intensity"][apex_idx] + # Threshold value + cutoff = intensity_threshold * apex_intensity + apex_fragment_name = df_sorted["fragment_name"][apex_idx] + + # Left of apex + left_df = df_sorted.filter(pl.col("fragment_name") == apex_fragment_name) # only consider the apex fragment + apex_idx_left = left_df["fragment_intensity"].arg_max() + left_df = left_df[:apex_idx_left][::-1] # reverse to go from apex down + left_bound = apex_rt + + for rt, intensity in zip(left_df["rt"], left_df["fragment_intensity"]): + if intensity < cutoff: + left_bound = rt + break + + # if the left bound is still the apex rt, set it to the first rt where fragment was detected + if left_bound == apex_rt and len(left_df) > 0: + left_bound = left_df["rt"][-1] + + # Right of apex + right_df = df_sorted.filter(pl.col("fragment_name") == apex_fragment_name) # only consider the apex fragment + apex_idx_right = right_df["fragment_intensity"].arg_max() + right_df = right_df[apex_idx_right+1:] + right_bound = apex_rt + for rt, intensity in zip(right_df["rt"], right_df["fragment_intensity"]): + if intensity < cutoff: + right_bound = rt + break + + # if the right bound is still the apex rt, set it to the last rt where fragment was detected + if right_bound == apex_rt and len(right_df) > 0: + right_bound = right_df["rt"][-1] + + # plot XIC with the margins + # plot_XIC_with_margins(df_sorted, output_dir=output_dir, adapted_interval=(left_bound, right_bound), apex_rt=apex_rt, cutoff=cutoff) + + return left_bound, right_bound, apex_rt + + +def calculate_min_max_margins(df_psms: pl.DataFrame, df_fragments: pl.DataFrame, top_n: int = 100, intensity_threshold: float = 0.01) -> dict: + """ + Calculate the retention time distribution of the top N peptidoforms (with at least 6 PSMs, and then ranked by spectrum peptide q value) + Min and max margins are defined as the 5th and 95th percentiles of the distribution of retention time margins + across the top N peptidoforms. + Returns a tuple with (min_diff, max_diff). + + Parameters + ---------- + df_psms : pl.DataFrame + DataFrame containing PSM information + df_fragments : pl.DataFrame + DataFrame containing fragment ion information + top_n : int, optional + Number of top peptidoforms to consider based on the lowest 'peptide_q' value (default is 100). + intensity_threshold : float, optional + Intensity threshold (as a fraction of apex intensity) to define retention time margins (default is 0.01). + """ + + # Step 1: Identify the 100 best scoring peptidoforms based on sage qvalue + # group by peptide and charge to get unique peptidoforms, aggregate number of PSMs, keep min peptide_q + + df_top_peptidoforms = ( + df_psms.group_by(["peptide", "charge"]) + .agg( + [pl.count().alias("num_psms"), pl.min("peptide_q").alias("min_peptide_q")] + ) + .sort("min_peptide_q") + ) + + # filter for peptidoforms with at least 6 PSMs + df_top_peptidoforms = df_top_peptidoforms.filter(pl.col("num_psms") >= 6) + + # get the top N peptidoforms + df_top_peptidoforms = df_top_peptidoforms.head(top_n) + + # Step 2: Extract the retention times of the entire XICs from df_fragments of these peptidoforms + df_fragments_top100 = df_fragments.filter(pl.col("peptide").is_in(df_top_peptidoforms["peptide"]) & pl.col("charge").is_in(df_top_peptidoforms["charge"])) + diffs = [] + + for (peptidoform, charge), df_fragments_top100_sub in tqdm( + df_fragments_top100.group_by(["peptide", "charge"]) + ): + left_bound, right_bound, apex_rt = calculate_rt_margins_intensity_based(df_fragments_top100_sub, intensity_threshold, output_dir='debug/calibration_xics') + left_diff = apex_rt - left_bound + right_diff = right_bound - apex_rt + diffs.append(left_diff) + diffs.append(right_diff) + + # remove 0 diffs (if the apex is at the start or end of the XIC) + diffs = [d for d in diffs if d > 0] + + # Step 3: Calculate the min and max retention times across all these XICs + if len(diffs) == 0: + log_info("Could not calibrate retention time margins, using default values.") + min_diff = 0.02 + max_diff = 0.2 + else: + # get 5th and 95th percentiles + min_diff = np.percentile(diffs, 5) + max_diff = np.percentile(diffs, 95) + log_info(f"Using min and max retention time margins: {min_diff}, {max_diff}") + + # plot histogram of diffs + plot_rt_margin_histogram(diffs, output_dir='debug/calibration_xics', min_diff=min_diff, max_diff=max_diff) + + return min_diff, max_diff + + +def add_retention_time_margins(df_psms: pl.DataFrame, df_fragment: pl.DataFrame, min_diff: float, max_diff: float, intensity_threshold: float) -> pl.DataFrame: + """ + Add retention time margin features to the PSM DataFrame. + """ + + pept2lowermargins = {} + pept2highermargins = {} + + log_info("Calculating adapted retention time margins based on intensity for all peptides") + + for (peptidoform, charge), df_fragments_sub in tqdm( + df_fragment.group_by(["peptide", "charge"]) + ): + + # speed up: skip peptidoforms with only 1 PSM + if df_fragments_sub['psm_id'].n_unique() < 2: + pept2lowermargins[(peptidoform, charge)] = np.nan + pept2highermargins[(peptidoform, charge)] = np.nan + continue + + intensity_based_margins = calculate_rt_margins_intensity_based(df_fragments_sub, intensity_threshold, output_dir='xics') + left_bound, right_bound, apex_rt = intensity_based_margins + + # check if the intensity based margins are higher than max or lower than min + left_diff = apex_rt - left_bound + right_diff = right_bound - apex_rt + + if left_diff < min_diff: + left_bound = apex_rt - min_diff + if right_diff < min_diff: + right_bound = apex_rt + min_diff + if left_diff > max_diff: + left_bound = apex_rt - max_diff + if right_diff > max_diff: + right_bound = apex_rt + max_diff + + pept2lowermargins[(peptidoform, charge)] = left_bound + pept2highermargins[(peptidoform, charge)] = right_bound + + log_info("Adding retention time margin features to PSM DataFrame...") + + # add rt_lower_margin and rt_higher_margin to df_psms + df_psms = df_psms.with_columns( + [ + pl.struct(["peptide", "charge"]) + .map_elements(lambda row: pept2lowermargins.get((row["peptide"], row["charge"]), np.nan)) + .alias("rt_lower_margin"), + pl.struct(["peptide", "charge"]) + .map_elements(lambda row: pept2highermargins.get((row["peptide"], row["charge"]), np.nan)) + .alias("rt_higher_margin") + ] + ) + + log_info("Adding retention time margin features to Fragment DataFrame...") + + # add rt_lower_margin and rt_higher_margin to df_fragment + df_fragment = df_fragment.with_columns( + [ + pl.struct(["peptide", "charge"]) + .map_elements(lambda row: pept2lowermargins.get((row["peptide"], row["charge"]), np.nan)) + .alias("rt_lower_margin"), + pl.struct(["peptide", "charge"]) + .map_elements(lambda row: pept2highermargins.get((row["peptide"], row["charge"]), np.nan)) + .alias("rt_higher_margin") + ] + ) + + return df_psms, df_fragment + + +def add_retention_time_margins_loop(df_psms: pl.DataFrame, df_fragment: pl.DataFrame, top_n: int = 10, intensity_threshold: float = 0.05) -> pl.DataFrame: + """ + Add retention time margin features to the PSM DataFrame. + """ + log_info("Calculating min max retention time margins based on intensity...") + # Step 1: Calculate min and max retention time window based on top 100 peptidoforms + min_diff, max_diff = calculate_min_max_margins(df_psms, df_fragment, top_n, intensity_threshold) + + # Step 2: Calculate adapted margins for each PSM based on the intensity of the most intense fragment + # and use the retention time distribution as min and max + log_info("Adding retention time margin features to PSM DataFrame...") + df_psms, df_fragment = add_retention_time_margins(df_psms, df_fragment, min_diff, max_diff, intensity_threshold) + + return df_psms, df_fragment + + def calculate_features( df_psms: pl.DataFrame, df_fragment: pl.DataFrame, @@ -993,6 +1222,23 @@ def calculate_features( .unique(subset=["peptide", "charge"], keep="first") ) + log_info("Regenerated df_fragment_max_peptide:") + log_info(" Shape: {}".format(df_fragment_max_peptide.shape)) + log_info(" Sample entries:") + #for row in df_fragment_max_peptide.head(3).iter_rows(named=True): + # log_info( + # " Peptide: {}, Charge: {}, PSM ID: {}, RT: {}, Fragment Intensity: {}".format( + # row["peptide"], + # row["charge"], + # row["psm_id"], + # row["rt"], + # row["fragment_intensity"], + # ) + # ) + + log_info( + "Counting individual peptides per MS2 and filtering by minimum occurrences" + ) df_psms = add_count_and_filter_peptides(df_psms, min_occurrences) # Filter df_fragment to only include PSMs that passed all filtering @@ -1147,6 +1393,7 @@ def calculate_features( f"{config['mumdia']['result_dir']}/outfile.pin", separator="\t" ) + return df_fragment, df_psms def main( df_fragment: Optional[pl.DataFrame] = None, @@ -1181,7 +1428,7 @@ def main( df_psms = df_psms.filter(~df_psms["peptide"].str.contains("U")) df_psms = df_psms.sort("rt") - calculate_features( + df_fragment, df_psms = calculate_features( df_psms, df_fragment, df_fragment_max, @@ -1193,8 +1440,21 @@ def main( ) log_info("Done running MuMDIA...") - # run_mokapot(output_dir=config["mumdia"]["result_dir"]) - + mokapot_results = run_mokapot(output_dir=config["mumdia"]["result_dir"]) + + df_fragment.write_csv("debug/df_fragment_before_quant.tsv", separator="\t") + df_psms.write_csv("debug/df_psms_before_quant.tsv", separator="\t") + + # this file will later be used for quantification of proteins with directLFQ (combined with all runs) + if mokapot_results is not None and isinstance(mokapot_results, (list, tuple)) and len(mokapot_results) > 1: + df_quant_fragment = quantify_fragments( + df_fragment, + mokapot_results[1], + config=config, + output_dir=config["mumdia"]["result_dir"] + ) + else: + logging.warning("mokapot_results is None or does not have enough elements; skipping quantification step.") if __name__ == "__main__": # In practice, load your input DataFrames (e.g., from parquet files) and then call main(). diff --git a/prediction_wrappers/wrapper_ms2pip.py b/prediction_wrappers/wrapper_ms2pip.py index f0c38c0..9456891 100644 --- a/prediction_wrappers/wrapper_ms2pip.py +++ b/prediction_wrappers/wrapper_ms2pip.py @@ -177,14 +177,4 @@ def get_predictions_fragment_intensity_main_loop( log_info("Df_fragment shape after filtering: {}".format(df_fragment.shape)) - df_fragment = df_fragment.with_columns( - pl.Series( - "fragment_name", - df_fragment["fragment_type"] - + df_fragment["fragment_ordinals"] - + "/" - + df_fragment["fragment_charge"], - ) - ) - return df_fragment, ms2pip_predictions diff --git a/quantification/lfq.py b/quantification/lfq.py new file mode 100644 index 0000000..d691317 --- /dev/null +++ b/quantification/lfq.py @@ -0,0 +1,162 @@ +import os +import logging + +import polars as pl +from tqdm import tqdm +from Bio import SeqIO +from pyteomics import proforma +import directlfq.config as lfqconfig +import directlfq.protein_intensity_estimation as lfqprot_estimation +import directlfq.utils as lfqutils +from utilities.logger import log_info + + +from scipy import integrate + + +def quantify_fragments(df_fragment, mokapot_psm_path, config, output_dir: str = None): + """Quantify fragment ions based on their intensities over retention time. + This function processes a DataFrame containing fragment ion data, filters it based on mokapot results, + and quantifies the fragment ions using both peak intensity and integrated intensity methods. + The results are saved to a CSV file in the specified output directory. + + TODO: These files should be combined for multiple runs and then used for directLFQ on fragment level. + TODO: Compare which quantification method works best, remove the others. + + Args: + df_fragment (pl.DataFrame): DataFrame containing fragment ion data with columns: + - 'peptide': Peptide sequence with modifications. + - 'stripped_peptide': Peptide sequence without modifications. + - 'charge': Charge state of the peptide. + - 'fragment_name': Name of the fragment ion (e.g., 'b3', 'y7'). + - 'proteins': Protein(s) the peptide maps to. + - 'rt': Retention time of the fragment ion measurement. + - 'fragment_intensity': Intensity of the fragment ion. + - 'psm_id': Unique identifier for the PSM. + - 'rt_lower_margin': Lower margin of the retention time window. + - 'rt_higher_margin': Upper margin of the retention time window. + mokapot_psm_path (str): Path to the mokapot PSM results file. + config (dict): Configuration dictionary containing at least the key 'sage' with subkey ' + mzml_paths' (list of str): List of mzML file paths. + output_dir (str, optional): Directory to save the output CSV file. Defaults to None. + Returns: + pl.DataFrame: DataFrame with quantified fragment ion intensities. + """ + # TODO: adapt this so it works for multiple runs + mzml_filename = os.path.basename(config["sage"]["mzml_paths"][0]).split('.')[0] + + #filter df_fragment to peptides that survived mokapot + mokapot_peptides = pl.read_csv(mokapot_psm_path, separator="\t") + mokapot_peptides = mokapot_peptides.filter(pl.col("mokapot q-value") < 0.01) + df_fragment_mokapot_filtered = df_fragment.filter(pl.col("peptide").is_in(mokapot_peptides["Peptide"])) + + results = [] + + logging.info(f"Quantifying fragments") + + for (peptidoform, stripped_sequence, charge, fragment_name, proteins), df_fragment_mokapot_filtered_sub in tqdm( + df_fragment_mokapot_filtered.group_by(["peptide", "stripped_peptide", "charge", "fragment_name", "proteins"]) + ): + proteins = ';'.join(proteinstring.split('|')[2] for proteinstring in proteins.split(';')) + + results.append({ + "protein": proteins, + "ion": "SEQ_" + stripped_sequence + "_MOD" + peptidoform + "_CHARGE_" + str(charge) + "_" + fragment_name, + #mzml_filename + "_Intensity_peak": quantify_fragment_peak_intensity(df_fragment_mokapot_filtered_sub, margin=False), + #mzml_filename + "_Intensity_peak_margin": quantify_fragment_peak_intensity(df_fragment_mokapot_filtered_sub, margin=True), + mzml_filename + "_Intensity_integrated": quantify_fragment_integrated_intensity(df_fragment_mokapot_filtered_sub, margin=False), + #mzml_filename + "_Intensity_integrated_margin": quantify_fragment_integrated_intensity(df_fragment_mokapot_filtered_sub, margin=True), + }) + + df_quant_fragment = pl.DataFrame(results) + df_quant_fragment.write_csv(os.path.join(output_dir, mzml_filename + "_fragment_level_intensities.csv")) + + return df_quant_fragment + +def quantify_fragment_peak_intensity(df_fragment_ion_psms, margin: bool = False): + # return highest intensity of fragment over RT + if margin: + df_fragment_ion_psms = _filter_margin_rt(df_fragment_ion_psms) + + return df_fragment_ion_psms["fragment_intensity"].max() + +def quantify_fragment_integrated_intensity(df_fragment_ion_psms, margin: bool = False): + # integrate fragment intensities over RT + if margin: + df_fragment_ion_psms = _filter_margin_rt(df_fragment_ion_psms) + + # for fragments only measured at one time point, return that intensity + if df_fragment_ion_psms.shape[0] == 1: + return df_fragment_ion_psms["fragment_intensity"].item() + + # sort by RT to ensure correct integration + df_fragment_ion_psms = df_fragment_ion_psms.sort("rt", descending=False) + + # approximate integration using trapezoidal rule + aoc = integrate.trapezoid( + y=df_fragment_ion_psms["fragment_intensity"].to_numpy(), + x=df_fragment_ion_psms["rt"].to_numpy() + ) + return aoc + +def _filter_margin_rt(df): + if df['psm_id'].n_unique() > 1: + df = df.filter( + (pl.col("rt") >= pl.col("rt_lower_margin")) & + (pl.col("rt") <= pl.col("rt_higher_margin")) + ) + return df + + +def quantify_proteins(df_fragment_quant_folder, output_dir: str = None): + """Estimate protein intensities from fragment ion quantifications using directLFQ. + This function is still TODO and not yet implemented. + + Args: + df_fragment_quant_folder (str): Path to the folder containing fragment ion quantification CSV files. + output_dir (str, optional): Directory to save the output CSV file. Defaults to None. + Returns: + directLFQ protein intensity DataFrame. + + """ + + # copy pasterino from alphadia: + #log_info.info("Performing label-free protein quantification using directLFQ") + + # combine all fragment quantification files in the folder + #df_results_fragment = pd.concat( + # [ + # pd.read_csv(os.path.join(df_fragment_quant_folder, f)) + # for f in os.listdir(df_fragment_quant_folder) + # if f.endswith("_fragment_level_intensities.csv") + # ], + # ignore_index=True, + #) + + # extract intensity columns + # intensity_cols = [col for col in df_results_fragment.columns if col.endswith("_Intensity_integrated_margin")] + #_intensity_df = df_results_fragment.select("Proteins", "ion", *intensity_cols) + + #lfqconfig.set_global_protein_and_ion_id(protein_id='Proteins', quant_id=mzml_filename + "_Intensity_integrated_margin") + #lfqconfig.set_compile_normalized_ion_table( + # compile_normalized_ion_table=False + #) # save compute time by avoiding the creation of a normalized ion table + #lfqconfig.check_wether_to_copy_numpy_arrays_derived_from_pandas() # avoid read-only pandas bug on linux if applicable + #lfqconfig.set_log_processed_proteins( + # log_processed_proteins=True + #) # here you can chose wether to log the processed proteins or not + + #_intensity_df.sort_values(by='_Intensity_integrated_margin', inplace=True, ignore_index=True) + + #lfq_df = lfqutils.index_and_log_transform_input_df(_intensity_df) + #lfq_df = lfqutils.remove_allnan_rows_input_df(lfq_df) + + #protein_df, _ = lfqprot_estimation.estimate_protein_intensities( + # lfq_df, + # min_nonan=1, + # num_samples_quadratic=50, + #) + + #protein_df.to_csv(os.path.join(output_dir, "protein_level_intensities_directlfq_out.csv")) + + #return protein_df \ No newline at end of file diff --git a/quantification/test_ipynbs/test_lfq.ipynb b/quantification/test_ipynbs/test_lfq.ipynb new file mode 100644 index 0000000..633f2ae --- /dev/null +++ b/quantification/test_ipynbs/test_lfq.ipynb @@ -0,0 +1,339 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "3752a435", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/caro/MuMDIA/.venv/lib/python3.11/site-packages/psims/mzmlb/writer.py:33: UserWarning: hdf5plugin is missing! Only the slower GZIP compression scheme will be available! Please install hdf5plugin to be able to use Blosc.\n", + " warnings.warn(\n" + ] + } + ], + "source": [ + "import sys\n", + "import os\n", + "import json\n", + "\n", + "from directlfq.lfq_manager import run_lfq\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import polars as pl\n", + "\n", + "sys.path.append('../..') \n", + "from quantification.lfq import quantify_fragments" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "8dc888ff", + "metadata": {}, + "outputs": [], + "source": [ + "def run_quant(condition: str):\n", + " df_fragment = pl.read_csv(\"../testfiles/condition{}_fullmzml/df_fragment_after_ms2pip.tsv\".format(condition), separator=\"\\t\")\n", + " mokapot_results = \"../testfiles/condition{}_fullmzml/mokapot.peptides.txt\".format(condition)\n", + " config = json.load(open(\"../../configs/config_condition{}_fullmzml.json\".format(condition)))\n", + "\n", + " quantify_fragments(df_fragment, mokapot_results, config, output_dir=\"quant_out/\")\n", + "\n", + "def assign_species(protein_str):\n", + " if \"HUMAN\" in protein_str:\n", + " return \"HUMAN\"\n", + " elif \"ECOLI\" in protein_str:\n", + " return \"ECOLI\"\n", + " elif \"YEAST\" in protein_str:\n", + " return \"YEAST\"\n", + " else:\n", + " print(\"Warning: Protein not assigned to any species: \", protein_str)\n", + " return \"OTHER\"\n", + "\n", + "def plot_logfc_histograms(directLFQ_protein_intensities):\n", + " # plot log2FC distribution by species\n", + " # create a new column \"Species\" based on the \"protein\" column\n", + " directLFQ_protein_intensities = directLFQ_protein_intensities.with_columns(\n", + " pl.col(\"protein\").map_elements(assign_species).alias(\"Species\")\n", + " )\n", + " species_colors = {\n", + " \"HUMAN\": \"green\",\n", + " \"ECOLI\": \"blue\",\n", + " \"YEAST\": \"red\",\n", + " \"OTHER\": \"gray\"\n", + " } \n", + "\n", + " quant_pd = directLFQ_protein_intensities.to_pandas()\n", + " plt.figure(figsize=(10, 6))\n", + " for species, color in species_colors.items():\n", + " subset = quant_pd[quant_pd[\"Species\"] == species]\n", + " plt.hist(subset[\"log2FC_A_vs_B\"], bins=150, alpha=0.5, label=species, color=color) \n", + " plt.xlabel(\"log2 Fold Change (A vs B)\")\n", + " plt.ylabel(\"Frequency\")\n", + " plt.title(\"Distribution of log2 Fold Changes by Species\")\n", + " plt.legend()\n", + " plt.axvline(x=1, color='red', linestyle='--')\n", + " plt.axvline(x=-2, color='blue', linestyle='--')\n", + " plt.axvline(x=0, color='green', linestyle='--')\n", + " plt.tight_layout()\n", + " plt.savefig(\"quant_out/log2FC_histogram_by_species.png\")\n", + " plt.show()\n", + "\n", + " return" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "d47630bd", + "metadata": {}, + "outputs": [], + "source": [ + "def quantify_proteins(df_fragment_quant_folder, output_dir: str = None):\n", + " \"\"\"Estimate protein intensities from fragment ion quantifications using directLFQ.\n", + "\n", + " Args:\n", + " df_fragment_quant_folder (str): Path to the folder containing fragment ion quantification CSV files.\n", + " output_dir (str, optional): Directory to save the output CSV files. Defaults to None.\n", + " \"\"\"\n", + " print(\"Performing label-free protein quantification using directLFQ\")\n", + "\n", + " # combine all fragment quantification files in the folder\n", + " fragment_files = [f for f in os.listdir(df_fragment_quant_folder) if f.endswith(\"_fragment_level_intensities.csv\")]\n", + "\n", + " if len(fragment_files) == 0:\n", + " raise ValueError(\"No fragment quantification files found in the specified folder.\")\n", + "\n", + " df_results_fragment = pd.read_csv(os.path.join(df_fragment_quant_folder, fragment_files[0]))\n", + "\n", + " if len(fragment_files) > 1:\n", + " # combine into one file\n", + " for f in fragment_files[1:]:\n", + " df_results_fragment = pd.merge(\n", + " df_results_fragment,\n", + " pd.read_csv(os.path.join(df_fragment_quant_folder, f)),\n", + " on=[\"protein\", \"ion\"],\n", + " how=\"inner\"\n", + " )\n", + "\n", + " # extract intensity columns\n", + " intensity_cols = [col for col in df_results_fragment.columns if col.endswith(\"_Intensity_integrated\")]\n", + " _intensity_df = df_results_fragment[[\"protein\", \"ion\", *intensity_cols]]\n", + "\n", + " # drop ions that match to multiple proteins\n", + " _intensity_df = _intensity_df[~_intensity_df[\"protein\"].str.contains(\";\")]\n", + "\n", + " # save combined fragment intensities\n", + " _intensity_df.to_csv(os.path.join(output_dir, \"combined_fragment_intensities.aq_reformat.tsv\"), index=False, sep=\"\\t\")\n", + "\n", + " # run directLFQ\n", + " run_lfq(\n", + " input_file=os.path.join(output_dir, \"combined_fragment_intensities.aq_reformat.tsv\"),\n", + " maximum_number_of_quadratic_ions_to_use_per_protein=10,\n", + " number_of_quadratic_samples=50\n", + " )\n", + "\n", + " return" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "a1a4ee12", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/tmp/ipykernel_1475011/136330505.py:4: ResourceWarning: unclosed file <_io.TextIOWrapper name='../../configs/config_conditionA_fullmzml.json' mode='r' encoding='UTF-8'>\n", + " config = json.load(open(\"../../configs/config_condition{}_fullmzml.json\".format(condition)))\n", + "ResourceWarning: Enable tracemalloc to get the object allocation traceback\n", + "/home/caro/MuMDIA/quantification/test_ipynbs/../../quantification/lfq.py:51: DeprecationWarning: `is_in` with a collection of the same datatype is ambiguous and deprecated.\n", + "Please use `implode` to return to previous behavior.\n", + "\n", + "See https://github.com/pola-rs/polars/issues/22149 for more information.\n", + " df_fragment_mokapot_filtered = df_fragment.filter(pl.col(\"peptide\").is_in(mokapot_peptides[\"Peptide\"]))\n", + "2025-09-26 15:37:47,531 - root - INFO - Quantifying fragments\n", + "177653it [03:53, 759.65it/s]\n", + "/tmp/ipykernel_1475011/136330505.py:4: ResourceWarning: unclosed file <_io.TextIOWrapper name='../../configs/config_conditionB_fullmzml.json' mode='r' encoding='UTF-8'>\n", + " config = json.load(open(\"../../configs/config_condition{}_fullmzml.json\".format(condition)))\n", + "ResourceWarning: Enable tracemalloc to get the object allocation traceback\n", + "2025-09-26 15:41:42,629 - root - INFO - Quantifying fragments\n", + "181325it [03:58, 760.05it/s]\n" + ] + } + ], + "source": [ + "for condition in [\"A\", \"B\"]:\n", + " run_quant(condition)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "269cf063", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Performing label-free protein quantification using directLFQ\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2025-09-26 15:45:43,132 - directlfq.lfq_manager - INFO - Starting directLFQ analysis.\n", + "2025-09-26 15:45:43,664 - directlfq.lfq_manager - INFO - Performing sample normalization.\n", + "2025-09-26 15:45:45,735 - directlfq.lfq_manager - INFO - Estimating lfq intensities.\n", + "2025-09-26 15:45:45,745 - directlfq.protein_intensity_estimation - INFO - 1957 lfq-groups total\n", + "2025-09-26 15:45:53,494 - directlfq.protein_intensity_estimation - INFO - using 60 processes\n", + "2025-09-26 15:45:53,539 - directlfq.protein_intensity_estimation - INFO - lfq-object 0\n", + "2025-09-26 15:45:53,827 - directlfq.protein_intensity_estimation - INFO - lfq-object 100\n", + "2025-09-26 15:45:53,883 - directlfq.protein_intensity_estimation - INFO - lfq-object 200\n", + "2025-09-26 15:45:54,029 - directlfq.protein_intensity_estimation - INFO - lfq-object 300\n", + "2025-09-26 15:45:54,115 - directlfq.protein_intensity_estimation - INFO - lfq-object 400\n", + "2025-09-26 15:45:54,157 - directlfq.protein_intensity_estimation - INFO - lfq-object 600\n", + "2025-09-26 15:45:54,209 - directlfq.protein_intensity_estimation - INFO - lfq-object 500\n", + "2025-09-26 15:45:54,346 - directlfq.protein_intensity_estimation - INFO - lfq-object 700\n", + "2025-09-26 15:45:54,379 - directlfq.protein_intensity_estimation - INFO - lfq-object 900\n", + "2025-09-26 15:45:54,398 - directlfq.protein_intensity_estimation - INFO - lfq-object 800\n", + "2025-09-26 15:45:54,479 - directlfq.protein_intensity_estimation - INFO - lfq-object 1000\n", + "2025-09-26 15:45:54,565 - directlfq.protein_intensity_estimation - INFO - lfq-object 1100\n", + "2025-09-26 15:45:54,669 - directlfq.protein_intensity_estimation - INFO - lfq-object 1200\n", + "2025-09-26 15:45:54,783 - directlfq.protein_intensity_estimation - INFO - lfq-object 1300\n", + "2025-09-26 15:45:54,960 - directlfq.protein_intensity_estimation - INFO - lfq-object 1500\n", + "2025-09-26 15:45:54,977 - directlfq.protein_intensity_estimation - INFO - lfq-object 1400\n", + "2025-09-26 15:45:55,036 - directlfq.protein_intensity_estimation - INFO - lfq-object 1600\n", + "2025-09-26 15:45:55,081 - directlfq.protein_intensity_estimation - INFO - lfq-object 1800\n", + "2025-09-26 15:45:55,159 - directlfq.protein_intensity_estimation - INFO - lfq-object 1700\n", + "2025-09-26 15:45:55,173 - directlfq.protein_intensity_estimation - INFO - lfq-object 1900\n", + "2025-09-26 15:45:56,087 - directlfq.lfq_manager - INFO - Could not add additional columns to protein table, printing without additional columns.\n", + "2025-09-26 15:45:56,088 - directlfq.lfq_manager - INFO - Writing results files.\n", + "2025-09-26 15:45:56,710 - directlfq.lfq_manager - INFO - Analysis finished!\n" + ] + } + ], + "source": [ + "# run directLFQ\n", + "quantify_proteins(\"quant_out/\", output_dir=\"quant_out/\")\n", + "\n", + "# to polars\n", + "directLFQ = pl.read_csv(\"quant_out/combined_fragment_intensities.aq_reformat.tsv.protein_intensities.tsv\", separator=\"\\t\")\n", + "\n", + "# add FCs\n", + "directLFQ = directLFQ.with_columns(\n", + " (pl.col(\"ttSCP_diaPASEF_Condition_A_Sample_Alpha_02_11500_Intensity_integrated\") / pl.col(\"ttSCP_diaPASEF_Condition_B_Sample_Alpha_02_11502_Intensity_integrated\")).log(2).alias(\"log2FC_A_vs_B\")\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "db27ab34", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/tmp/ipykernel_1475011/136330505.py:22: MapWithoutReturnDtypeWarning: 'return_dtype' of function python_udf must be set\n", + "\n", + "A later expression might fail because the output type is not known. Set return_dtype=pl.self_dtype() if the type is unchanged, or set the proper output data type.\n", + " directLFQ_protein_intensities = directLFQ_protein_intensities.with_columns(\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n", + "sys:1: MapWithoutReturnDtypeWarning: Calling `map_elements` without specifying `return_dtype` can lead to unpredictable results. Specify `return_dtype` to silence this warning.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Warning: Protein not assigned to any species: PEDF_BOVIN\n", + "Warning: Protein not assigned to any species: FA5_BOVIN\n", + "Warning: Protein not assigned to any species: FETUA_BOVIN\n", + "Warning: Protein not assigned to any species: G6PI_BOVIN\n", + "Warning: Protein not assigned to any species: HRG_BOVIN\n", + "Warning: Protein not assigned to any species: TRFE_BOVIN\n", + "Warning: Protein not assigned to any species: TRYP_PIG\n", + "Warning: Protein not assigned to any species: A0A3Q1M3L6_BOVIN\n", + "Warning: Protein not assigned to any species: A1AG_BOVIN\n", + "Warning: Protein not assigned to any species: A2MG_BOVIN\n", + "Warning: Protein not assigned to any species: ALBU_BOVIN\n", + "Warning: Protein not assigned to any species: APOA1_BOVIN\n", + "Warning: Protein not assigned to any species: APOE_BOVIN\n", + "Warning: Protein not assigned to any species: Q1RMN8_BOVIN\n", + "Warning: Protein not assigned to any species: Q1WEI2_PSEAI\n", + "Warning: Protein not assigned to any species: HBA_BOVIN\n", + "Warning: Protein not assigned to any species: HBBF_BOVIN\n", + "Warning: Protein not assigned to any species: CAP1_BOVIN\n", + "Warning: Protein not assigned to any species: ITIH4_BOVIN\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAJOCAYAAACqS2TfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAfMpJREFUeJzt3XmcjeX/x/H3OWf2wYx9Zmwz1iwhSxJlZCdbESJLQl8qS1HTZvgWbYpIpa9QkZJSP4USg0oKUSRFlsrYvpYxM2a/f3/Mdw7HzJgx5j73zPF6Ph7343Gd677v6/6cc1/j+Jzruu/bZhiGIQAAAAAAUOjsVgcAAAAAAICnIukGAAAAAMAkJN0AAAAAAJiEpBsAAAAAAJOQdAMAAAAAYBKSbgAAAAAATELSDQAAAACASUi6AQAAAAAwCUk3AAAAAAAmIekGAAtER0fLZrO55ViRkZGKjIx0vo6JiZHNZtNHH33kluMPHTpU4eHhbjlWQcXHx+u+++5TSEiIbDabxo0bl+u24eHhGjp0qNtis9ql/Sc3Wf0qJiamUI9vs9n0wAMPFGqbxc3Bgwdls9n00ksvWR2K5Ww2m6Kjo60OAwCuCEk3AFylhQsXymazORc/Pz+FhYWpU6dOevXVV3Xu3LlCOc6RI0cUHR2tHTt2FEp7hakox5Yf06ZN08KFC/Wvf/1L7777ru655x5L4/n444/Vr18/Va9eXQEBAapTp44efvhhnTlzJl/7R0ZGuvTJi5fffvvN3ODzaf/+/Ro1apSqV68uPz8/lSpVSq1atdKsWbN0/vx5q8O7Zv3yyy/q06ePqlWrJj8/P1WqVEkdOnTQ7NmzrQ4NAIotL6sDAABPMXXqVEVERCg1NVVHjx5VTEyMxo0bp5dfflmfffaZGjZs6Nz2ySef1GOPPXZF7R85ckRTpkxReHi4GjdunO/9vvzyyys6TkFcLra33npLGRkZpsdwNdatW6ebbrpJkydPtjoUSdLIkSMVFhamQYMGqWrVqvrll180Z84cffHFF9q+fbv8/f3zbKNy5cqaPn16tvqwsDAzQr4in3/+ufr27StfX18NHjxYDRo0UEpKir755htNnDhRu3fv1rx586wO85rz3XffqW3btqpatapGjBihkJAQ/fXXX/r+++81a9YsPfjgg1aHqPPnz8vLi/++Aihe+FcLAApJly5d1KxZM+frqKgorVu3Trfffrt69OihPXv2OJMlLy8v0//jmJiYqICAAPn4+Jh6nLx4e3tbevz8OH78uOrVq2d1GE4fffRRtindTZs21ZAhQ7R48WLdd999ebYRFBSkQYMGmRRhwR04cED9+/dXtWrVtG7dOoWGhjrXjRkzRvv27dPnn39uYYTXrmeffVZBQUH68ccfFRwc7LLu+PHj1gR1CT8/P6tDAIArxvRyADDRbbfdpqeeekqHDh3Se++956zP6Zrur776Sq1bt1ZwcLBKlCihOnXq6PHHH5eUeb1s8+bNJUnDhg1zThVeuHChpMzpxA0aNNC2bdt06623KiAgwLlvbtfkpqen6/HHH1dISIgCAwPVo0cP/fXXXy7b5Hb98sVt5hVbTtd0JyQk6OGHH1aVKlXk6+urOnXq6KWXXpJhGC7bZV3Pu2LFCjVo0EC+vr6qX7++Vq9enfMHfonjx49r+PDhqlixovz8/NSoUSMtWrTIuT7rOuQDBw7o888/d8Z+8ODBfLWf5c8//1Tfvn1VpkwZBQQE6KabbsoxcTx06JB69OihwMBAVahQQePHj9eaNWuyXQud0/nq3bu3JGnPnj1XFFtu0tLS9O9//1s1atSQr6+vwsPD9fjjjys5OTnPff/++2/16tXL5X3kZz9JeuGFFxQfH6/58+e7JNxZatasqbFjx2arz6sPHDp0SKNHj1adOnXk7++vsmXLqm/fvtnOZdblIN9++60mTJig8uXLKzAwUL1799aJEydcts3IyFB0dLTCwsIUEBCgtm3b6tdff83x7+LMmTMaN26cs0/XrFlTzz//fLZZHkuXLlXTpk1VsmRJlSpVStdff71mzZqVr89Okl555RVVq1ZN/v7+atOmjXbt2uVct2DBAtlsNv3000/Z9ps2bZocDof++eefXNvev3+/6tevny3hlqQKFSq4vM7621y8eLHq1KkjPz8/NW3aVBs3bsy27z///KN7771XFStWdJ6/t99+O9t2SUlJio6OVu3ateXn56fQ0FDdcccd2r9/v8txL72mO7/tz549W/Xr11dAQIBKly6tZs2aacmSJbl+HgBQWBjpBgCT3XPPPXr88cf15ZdfasSIETlus3v3bt1+++1q2LChpk6dKl9fX+3bt0/ffvutJKlu3bqaOnWqnn76aY0cOVK33HKLJOnmm292tvHf//5XXbp0Uf/+/TVo0CBVrFjxsnE9++yzstlsevTRR3X8+HHNnDlT7du3144dO/I1fTlLfmK7mGEY6tGjh9avX6/hw4ercePGWrNmjSZOnKh//vlHr7zyisv233zzjT7++GONHj1aJUuW1Kuvvqo777xThw8fVtmyZXON6/z584qMjNS+ffv0wAMPKCIiQsuWLdPQoUN15swZjR07VnXr1tW7776r8ePHq3Llynr44YclSeXLl8/3+z927JhuvvlmJSYm6qGHHlLZsmW1aNEi9ejRQx999JEzWU5ISNBtt92m2NhYjR07ViEhIVqyZInWr1+fr+McPXpUklSuXLl8bZ+enq6TJ0+61Pn5+alEiRKSpPvuu0+LFi1Snz599PDDD2vLli2aPn269uzZo08++STXds+fP6927drp8OHDeuihhxQWFqZ3331X69aty1dc//d//6fq1avn2j9ykp8+8OOPP+q7775T//79VblyZR08eFCvv/66IiMj9euvvyogIMClzQcffFClS5fW5MmTdfDgQc2cOVMPPPCAPvjgA+c2UVFReuGFF9S9e3d16tRJO3fuVKdOnZSUlOTSVmJiotq0aaN//vlHo0aNUtWqVfXdd98pKipKsbGxmjlzpqTMH9YGDBigdu3a6fnnn5eU+SPKt99+m+MPDZd65513dO7cOY0ZM0ZJSUmaNWuWbrvtNv3yyy+qWLGi+vTpozFjxmjx4sW64YYbXPZdvHixIiMjValSpVzbr1atmjZv3qxdu3apQYMGecazYcMGffDBB3rooYfk6+uruXPnqnPnzvrhhx+c+x87dkw33XSTM0kvX768Vq1apeHDhysuLs5508L09HTdfvvt+vrrr9W/f3+NHTtW586d01dffaVdu3apRo0aOcaQ3/bfeustPfTQQ+rTp4/Gjh2rpKQk/fzzz9qyZYvuvvvuPN8rAFwVAwBwVRYsWGBIMn788cdctwkKCjJuuOEG5+vJkycbF/8T/MorrxiSjBMnTuTaxo8//mhIMhYsWJBtXZs2bQxJxhtvvJHjujZt2jhfr1+/3pBkVKpUyYiLi3PWf/jhh4YkY9asWc66atWqGUOGDMmzzcvFNmTIEKNatWrO1ytWrDAkGc8884zLdn369DFsNpuxb98+Z50kw8fHx6Vu586dhiRj9uzZ2Y51sZkzZxqSjPfee89Zl5KSYrRs2dIoUaKEy3uvVq2a0a1bt8u2d/G2F38m48aNMyQZmzZtctadO3fOiIiIMMLDw4309HTDMAxjxowZhiRjxYoVzu3Onz9vXHfddYYkY/369Zc97vDhww2Hw2H8/vvvecaY1R8uXbLi3rFjhyHJuO+++1z2e+SRRwxJxrp161zauvhcZ32uH374obMuISHBqFmzZp7v4+zZs4Yko2fPnnm+hyz57QOJiYnZ9t28ebMhyXjnnXecdVl/r+3btzcyMjKc9ePHjzccDodx5swZwzAM4+jRo4aXl5fRq1cvlzajo6NdPkvDMIx///vfRmBgYLZz89hjjxkOh8M4fPiwYRiGMXbsWKNUqVJGWlpavt+/YRjGgQMHDEmGv7+/8ffffzvrt2zZYkgyxo8f76wbMGCAERYW5ux3hmEY27dvz/Xv82Jffvml4XA4DIfDYbRs2dKYNGmSsWbNGiMlJSXbtll9auvWrc66Q4cOGX5+fkbv3r2ddcOHDzdCQ0ONkydPuuzfv39/IygoyHne3n77bUOS8fLLL2c71sXnSZIxefLkK26/Z8+eRv369S/7/gHALEwvBwA3KFGixGXvYp41nfPTTz8t8E3HfH19NWzYsHxvP3jwYJUsWdL5uk+fPgoNDdUXX3xRoOPn1xdffCGHw6GHHnrIpf7hhx+WYRhatWqVS3379u1dRrkaNmyoUqVK6c8//8zzOCEhIRowYICzztvbWw899JDi4+O1YcOGQng3mce58cYb1bp1a2ddiRIlNHLkSB08eFC//vqrJGn16tWqVKmSevTo4dzOz88v19kPF1uyZInmz5+vhx9+WLVq1cpXXOHh4frqq69clkmTJjljlqQJEya47JM10n+5a6q/+OILhYaGqk+fPs66gIAAjRw5Ms+Y4uLiJMml3+VHfvrAxbMzUlNT9d///lc1a9ZUcHCwtm/fnq3NkSNHulziccsttyg9PV2HDh2SJH399ddKS0vT6NGjXfbL6WZiy5Yt0y233KLSpUvr5MmTzqV9+/ZKT093TrkODg5WQkKCvvrqqyt6/1l69erlMlJ94403qkWLFi5/s4MHD9aRI0dcZlAsXrxY/v7+uvPOOy/bfocOHbR582b16NFDO3fu1AsvvKBOnTqpUqVK+uyzz7Jt37JlSzVt2tT5umrVqurZs6fWrFmj9PR0GYah5cuXq3v37jIMw+Wz6dSpk86ePes8N8uXL1e5cuVy/Hxze7zilbQfHBysv//+Wz/++ONlPwMAMANJNwC4QXx8/GUTjX79+qlVq1a67777VLFiRfXv318ffvjhFSXglSpVuqKbpl2avNlsNtWsWfOKr2e+UocOHVJYWFi2z6Nu3brO9RerWrVqtjZKly6t06dP53mcWrVqyW53/arL7TgFdejQIdWpUydb/aXHOXTokGrUqJEtgahZs+Zl29+0aZOGDx+uTp066dlnn813XIGBgWrfvr3LknWzuEOHDslut2c7dkhIiIKDgy/72Rw6dEg1a9bM9j5y+gwuVapUKUm64sfo5acPnD9/Xk8//bTzmupy5cqpfPnyOnPmjM6ePZtnm6VLl5YkZ5tZn8Gln1GZMmWc22b5448/tHr1apUvX95lad++vaQLNyEbPXq0ateurS5duqhy5cq69957831/Ain736wk1a5d2+VvtkOHDgoNDdXixYslZV6X/v7776tnz575+rGjefPm+vjjj3X69Gn98MMPioqK0rlz59SnTx/nD0h5xZOYmKgTJ07oxIkTOnPmjObNm5fts8n6gTDrs9m/f7/q1KlzRTeYvJL2H330UZUoUUI33nijatWqpTFjxjgv3wEAs3FNNwCY7O+//9bZs2cvm1z5+/tr48aNWr9+vT7//HOtXr1aH3zwgW677TZ9+eWXcjgceR7nSq7Dzq/cRpjS09PzFVNhyO04xiU3XfNEO3fuVI8ePdSgQQN99NFHhX7H+9zOr1lKlSqlsLAwl5t/5Ud++sCDDz6oBQsWaNy4cWrZsqWCgoJks9nUv3//HH+8Ksx+lZGRoQ4dOjhnElyqdu3akjJvRrZjxw6tWbNGq1at0qpVq7RgwQINHjzY5QZ/V8PhcOjuu+/WW2+9pblz5+rbb7/VkSNHrvhO9j4+PmrevLmaN2+u2rVra9iwYVq2bNkVPVYv63MfNGiQhgwZkuM2Fz9K8UpdSft169bV3r17tXLlSq1evVrLly/X3Llz9fTTT2vKlCkFjgEA8oOkGwBM9u6770qSOnXqdNnt7Ha72rVrp3bt2unll1/WtGnT9MQTT2j9+vVq3759oSdIf/zxh8trwzC0b98+l/8Ely5dWmfOnMm276FDh1S9enXn6yuJrVq1alq7dq3OnTvnMvL222+/OdcXhmrVqunnn39WRkaGy2i3GcfZu3dvtvpLj1OtWjX9+uuvMgzD5fPat29fju3u379fnTt3VoUKFfTFF184b4BWWDFnZGTojz/+cI7IS5k3pTpz5sxlP5tq1app165d2d5HTp9BTm6//XbNmzdPmzdvVsuWLQv+Ji7x0UcfaciQIZoxY4azLikpKcf+mx9Zn8G+ffsUERHhrP/vf/+bbZZFjRo1FB8f7xzZvhwfHx91795d3bt3V0ZGhkaPHq0333xTTz31VJ6zHi79m5Wk33//PdvTAQYPHqwZM2bo//7v/7Rq1SqVL18+z39/LifrUYixsbH5iicgIMB5M8KSJUsqPT09z8+mRo0a2rJli1JTU/P9mMHy5cvnu30pc/ZHv3791K9fP6WkpOiOO+7Qs88+q6ioKB5FBsBUTC8HABOtW7dO//73vxUREaGBAwfmut2pU6ey1TVu3FiSnI9iCgwMlKQCJxGXyroTcpaPPvpIsbGx6tKli7OuRo0a+v7775WSkuKsW7lyZbZHi11JbF27dlV6errmzJnjUv/KK6/IZrO5HP9qdO3aVUePHnW5G3VaWppmz56tEiVKqE2bNoV2nB9++EGbN2921iUkJGjevHkKDw93Tunu1KmT/vnnH5drY5OSkvTWW29la/Po0aPq2LGj7Ha71qxZc0V3U89vzJKcd9XO8vLLL0uSunXrdtl9jxw5oo8++shZl5iYqHnz5uXr2JMmTVJgYKDuu+8+HTt2LNv6/fv3X9EjtLI4HI5so9SzZ89Wenr6FbclSe3atZOXl5def/11l/pL+60k3XXXXdq8ebPWrFmTbd2ZM2eUlpYmKTNhv5jdbnf+yJWfR66tWLHC5ZFfP/zwg7Zs2ZLtb6Zhw4Zq2LCh/vOf/2j58uXq379/vmZJrF+/PseR/qxrxi+9hGDz5s0u18v/9ddf+vTTT9WxY0c5HA45HA7deeedWr58eY6zGy5+RNudd96pkydP5vj55jb74Erav/Sz9/HxUb169WQYhlJTU3NsHwAKCyPdAFBIVq1apd9++01paWk6duyY1q1bp6+++krVqlXTZ599dtmRlKlTp2rjxo3q1q2bqlWrpuPHj2vu3LmqXLmy8wZdNWrUUHBwsN544w2VLFlSgYGBatGihcso3JUoU6aMWrdurWHDhunYsWOaOXOmatas6XJjr/vuu08fffSROnfurLvuukv79+/Xe++9l+3xPVcSW/fu3dW2bVs98cQTOnjwoBo1aqQvv/xSn376qcaNG5fro4Gu1MiRI/Xmm29q6NCh2rZtm8LDw/XRRx/p22+/1cyZM6/4Zl65eeyxx/T++++rS5cueuihh1SmTBktWrRIBw4c0PLly52j7KNGjdKcOXM0YMAAjR071nndbVa/uHjUuHPnzvrzzz81adIkffPNN/rmm2+c6ypWrKgOHTpcVcyNGjXSkCFDNG/ePJ05c0Zt2rTRDz/8oEWLFqlXr15q27ZtrvuOGDFCc+bM0eDBg7Vt2zaFhobq3XffzfZIrtzUqFFDS5YsUb9+/VS3bl0NHjxYDRo0UEpKir777jvnY92u1O233653331XQUFBqlevnjZv3qy1a9de9rFyl1OxYkWNHTtWM2bMUI8ePdS5c2ft3LlTq1atUrly5VzO18SJE/XZZ5/p9ttv19ChQ9W0aVMlJCTol19+0UcffaSDBw+qXLlyuu+++3Tq1Cnddtttqly5sg4dOqTZs2ercePGLjMOclOzZk21bt1a//rXv5ScnKyZM2eqbNmyOU5rHzx4sB555BFJyvfU8gcffFCJiYnq3bu3rrvuOuc5+eCDDxQeHp7tRo0NGjRQp06dXB4ZJslluvZzzz2n9evXq0WLFhoxYoTq1aunU6dOafv27Vq7dq3zB8fBgwfrnXfe0YQJE/TDDz/olltuUUJCgtauXavRo0erZ8+eOcac3/Y7duyokJAQtWrVShUrVtSePXs0Z84cdevWrdD+LQCAXFlwx3QA8ChZjyDKWnx8fIyQkBCjQ4cOxqxZs1weTZXl0keGff3110bPnj2NsLAww8fHxwgLCzMGDBiQ7RFEn376qVGvXj3Dy8vL5RFAbdq0yfVxOLk9Muz99983oqKijAoVKhj+/v5Gt27djEOHDmXbf8aMGUalSpUMX19fo1WrVsbWrVuztXm52C59ZJhhZD5Sa/z48UZYWJjh7e1t1KpVy3jxxRddHg1kGJmPBxozZky2mHJ7lNmljh07ZgwbNswoV66c4ePjY1x//fU5Pjbpah4ZZhiGsX//fqNPnz5GcHCw4efnZ9x4443GypUrs+37559/Gt26dTP8/f2N8uXLGw8//LCxfPlyQ5Lx/fffO7e7uD9dulz6uefkcv0hS2pqqjFlyhQjIiLC8Pb2NqpUqWJERUUZSUlJ2dq69JiHDh0yevToYQQEBBjlypUzxo4da6xevTpfjz7L8vvvvxsjRowwwsPDDR8fH6NkyZJGq1atjNmzZ7vEkN8+cPr0aee5LlGihNGpUyfjt99+y7Zdbo/4y/q7uDj+tLQ046mnnjJCQkIMf39/47bbbjP27NljlC1b1rj//vtd9j937pwRFRVl1KxZ0/Dx8THKlStn3HzzzcZLL73kfOTWRx99ZHTs2NGoUKGC4ePjY1StWtUYNWqUERsbe9nPKuuRYS+++KIxY8YMo0qVKoavr69xyy23GDt37sxxn9jYWMPhcBi1a9e+bNsXW7VqlXHvvfca1113nVGiRAnDx8fHqFmzpvHggw8ax44dc9k267y89957Rq1atQxfX1/jhhtuyPH8Hzt2zBgzZoxRpUoVw9vb2wgJCTHatWtnzJs3z2W7xMRE44knnnD2yZCQEKNPnz7G/v37XY578SPD8tv+m2++adx6661G2bJlDV9fX6NGjRrGxIkTjbNnz+b78wGAgrIZxjVwJxoAAIqomTNnavz48fr7779dHgeFounMmTMqXbq0nnnmGT3xxBNWh5OrkydPKjQ0VE8//bSeeuqpQm/fZrNpzJgxOU4HBwC44ppuAADc5Pz58y6vk5KS9Oabb6pWrVok3EXQpedLunAdfGRkpHuDuUILFy5Uenq67rnnHqtDAYBrHtd0AwDgJnfccYeqVq2qxo0b6+zZs3rvvff022+/OZ+pjKLlgw8+0MKFC9W1a1eVKFFC33zzjd5//3117NhRrVq1sjq8HK1bt06//vqrnn32WfXq1Svbnc0BAO5H0g0AgJt06tRJ//nPf7R48WKlp6erXr16Wrp0qfr162d1aMhBw4YN5eXlpRdeeEFxcXHOm6s988wzVoeWq6lTp+q7775Tq1atNHv2bKvDAQBI4ppuAAAAAABMwjXdAAAAAACYhKQbAAAAAACTcE23pIyMDB05ckQlS5aUzWazOhwAAAAAQBFnGIbOnTunsLAw2e25j2eTdEs6cuSIqlSpYnUYAAAAAIBi5q+//lLlypVzXU/SLalkyZKSMj+sUqVKWRwNABQ9CQlSWFhm+cgRKTDQ2nhQ+BJSEhQ2I/MkH3n4iAJ9PPwk06kBAFcpLi5OVapUceaTuSHplpxTykuVKkXSDQA5cDgulEuVIj/xRI4Uh+SXWS5VqpTnJ910agBAIcnrEmVupAYAAAAAgElIugEAAAAAMAlJNwAAAAAAJuGabgAAAAAoIjIyMpSSkmJ1GJDk7e0tx8X3ACkgkm4AAAAAKAJSUlJ04MABZWRkWB0K/ic4OFghISF53iztcki6AQB58veXDhy4UIbn8ff214GxB5xlj0enBlDEGIah2NhYORwOValSRXY7VwJbyTAMJSYm6vjx45Kk0NDQArdF0g0AyJPdLoWHWx0FzGS32RUeHG51GO5DpwZQxKSlpSkxMVFhYWEKCAiwOhxI8v/fj7LHjx9XhQoVCjzVnJ9PAAAAAMBi6enpkiQfHx+LI8HFsn4ASU1NLXAbJN0AgDylpEgTJ2Yu3NvFM6Wkp2jilxM18cuJSkm/Bk4ynRpAEXU11w6j8BXG+bAZhmEUQizFWlxcnIKCgnT27FmVKlXK6nAAoMhJSJBKlMgsx8dLgYHWxoPCl5CSoBLTM09yfFS8An08/CTTqQEUMUlJSTpw4IAiIiLk5+dndTj4n8udl/zmkYx0AwAAAABgEm6kBgAAAABFVHRMtHuPF3llxxs6dKjOnDmjFStWuNTHxMSobdu2On36tFasWKFx48bpzJkz2fa32Wz65JNP1KtXLx08eFARERGy2+06fPiwKlWq5NwuNjZWVapUUXp6ug4cOKDwS26G2alTJ61du1bff/+9mjdvni3GRYsWafr06Xrsscec9StWrFDv3r1l9uRvRroBAAAAAEVGpUqV9M4777jULVq0yCUJv9jhw4f13Xff6YEHHtDbb7+d4zZ+fn56/vnndfr06UKPNy8k3QAAAACAImPIkCFasGCBS92CBQs0ZMiQHLdfsGCBbr/9dv3rX//S+++/r/Pnz2fbpn379goJCdH06dNNiflySLoBAAAAAEVGjx49dPr0aX3zzTeSpG+++UanT59W9+7ds21rGIYWLFigQYMG6brrrlPNmjX10UcfZdvO4XBo2rRpmj17tv7++2/T38PFSLoBAAAAAAW2cuVKlShRwmXp0qVLgdvz9vbWoEGDnFPF3377bQ0aNEje3t7Ztl27dq0SExPVqVMnSdKgQYM0f/78HNvt3bu3GjdurMmTJxc4toIg6QYA5MnfX9q1K3Px97c6GpjB39tfu/61S7v+tUv+3tfASaZTA0Chadu2rXbs2OGy/Oc//7mqNu+9914tW7ZMR48e1bJly3TvvffmuN3bb7+tfv36ycsr8x7hAwYM0Lfffqv9+/fnuP3zzz+vRYsWac+ePVcV35Ug6QYA5Mlul+rXz1zsfHN4JLvNrvoV6qt+hfqy266Bk0ynBoBCExgYqJo1a7osF9/0rFSpUkpISFBGRobLfll3Mw8KCsrW5vXXX6/rrrtOAwYMUN26ddWgQYNs25w6dUqffPKJ5s6dKy8vL3l5ealSpUpKS0vL9YZqt956qzp16qSoqKireMdXhm8ZAAAAAIBp6tSpo7S0NO3YscOlfvv27ZKk2rVr57jfvffeq5iYmFxHuRcvXqzKlStr586dLqPsM2bM0MKFC5Wenp7jfs8995z+7//+T5s3by74m7oCPKcbAJCnlBRp2rTM8uOPSz4+1saDwpeSnqJpmzJP8uO3PC4fh4efZDo1ALhN/fr11bFjR917772aMWOGqlevrr1792rcuHHq169fro8CGzFihPr27avg4OAc18+fP199+vTJNgpepUoVRUVFafXq1erWrVu2/a6//noNHDhQr7766lW/t/xgpBsAkKfUVGnKlMwlNdXqaGCG1PRUTdkwRVM2TFFq+jVwkunUAOBWH3zwgdq0aaNRo0apfv36euihh9SzZ8/LXvvt5eWlcuXKOa/Xvti2bdu0c+dO3XnnndnWBQUFqV27drneUE2Spk6dmm26u1lshmEYbjlSERYXF6egoCCdPXtWpUqVsjocAChyEhKkEiUyy/HxUmCgtfEg/6JjonOuj3StT0hJUInpmSc5PipegT4efpLp1ACKmKSkJB04cEARERHy8/OzOhz8z+XOS37zSEa6AQAAAAAwCUk3AAAAAAAmIekGAAAAAMAkJN0AAAAAAJiEpBsAAAAAAJPwnG4AQJ78/KQffrhQhufx8/LTD/f94Cx7PDo1AMBNSLoBAHlyOKTmza2OAmZy2B1qXukaOsl0agCAmzC9HAAAAAAAkzDSDQDIU0qKNGtWZnnsWMnHx9p4UPhS0lM06/vMkzz2prHycXj4SaZTAwDchKQbAJCn1FRp0qTM8ujR5CeeKDU9VZPWZp7k0c1He37STacGALgJSTcAAAAAFFHR0UX7eEOHDtWiRYuy1Xfq1EmrV6+WJP3000+aNm2aNm7cqLNnz6pKlSqKjIzUxIkTVbt2bec+ixYt0pw5c7R79245HA41adJEEydO1O233+7cJiYmRm3bttXp06cVHBycQ/zRWrFihXbs2HFlb8REXNMNAAAAACiwzp07KzY21mV5//33JUkrV67UTTfdpOTkZC1evFh79uzRe++9p6CgID311FPONh555BGNGjVK/fr1088//6wffvhBrVu3Vs+ePTVnzhyr3lqhYKQbAAAAAFBgvr6+CgkJyVafmJioYcOGqWvXrvrkk0+c9REREWrRooXOnDkjSfr+++81Y8YMvfrqq3rwwQed2z377LNKSkrShAkT1LNnT1WpUsX092IGS0e6N27cqO7duyssLEw2m00rVqxwWW+z2XJcXnzxRec24eHh2dY/99xzbn4nAAAAAICLrVmzRidPntSkrHtoXCJrevj777+vEiVKaNSoUdm2efjhh5Wamqrly5ebGaqpLB3pTkhIUKNGjXTvvffqjjvuyLY+NjbW5fWqVas0fPhw3XnnnS71U6dO1YgRI5yvS5YsaU7AAAAAAAAXK1euVIkSJVzqHn/8cXl5Zaab11133WX3//3331WjRg355HBTy7CwMJUqVUq///574QXsZpYm3V26dFGXLl1yXX/pFIVPP/1Ubdu2VfXq1V3qS5YsmeN0BgAAAACAudq2bavXX3/dpa5MmTJ666238t2GYRiFHVaRUWyu6T527Jg+//zzHO+M99xzz+nf//63qlatqrvvvlvjx493/qoCALh6fn7S+vUXyvA8fl5+Wj9kvbPs8ejUAFBoAgMDVbNmzWz1WXcm/+2339SyZctc969du7a++eYbpaSkZBvtPnLkiOLi4lzucl7cFJu7ly9atEglS5bMNg39oYce0tKlS7V+/XqNGjVK06ZNy/WagSzJycmKi4tzWQAAuXM4pMjIzMXhsDoamMFhdygyPFKR4ZFy2K+Bk0ynBgDTdezYUeXKldMLL7yQ4/qsG6n1799f8fHxevPNN7Nt89JLL8nb2zvbJcbFSbEZDn777bc1cOBA+V3ya/SECROc5YYNG8rHx0ejRo3S9OnT5evrm2Nb06dP15QpU0yNFwAAAACuBcnJyTp69KhLnZeXl8qVK6f//Oc/6tu3r3r06KGHHnpINWvW1MmTJ/Xhhx/q8OHDWrp0qVq2bKmxY8dq4sSJSklJUa9evZSamqr33ntPs2bN0syZM7PdufyXX35xuZeXzWZTo0aN3PJ+r1SxSLo3bdqkvXv36oMPPshz2xYtWigtLU0HDx5UnTp1ctwmKirKJVmPi4srtrefBwB3SE2V5s3LLI8cKXl7WxsPCl9qeqrmbcs8ySObjpS3w8NPMp0aAArN6tWrFRoa6lJXp04d/fbbb+rZs6e+++47TZ8+XXfffbcz97rtttv0zDPPOLefOXOmGjZsqLlz5+rJJ5+Uw+FQkyZNtGLFCnXv3j3bMW+99VaX1w6HQ2lpaea8watkM4rIFes2m02ffPKJevXqlW3d0KFDtWvXLm3dujXPdhYvXqzBgwfr5MmTKl26dL6OHRcXp6CgIJ09e1alSpW60tABwOMlJEhZNyWNj5cCA62NB/kXHROdc32ka31CSoJKTM88yfFR8Qr08fCTTKcGUMQkJSXpwIEDioiIyDa7F9a53HnJbx5p6Uh3fHy89u3b53x94MAB7dixQ2XKlFHVqlUlZb6RZcuWacaMGdn237x5s7Zs2aK2bduqZMmS2rx5s8aPH69BgwblO+EGAAAAAMAslibdW7duVdu2bZ2vs6Z8DxkyRAsXLpQkLV26VIZhaMCAAdn29/X11dKlSxUdHa3k5GRFRERo/PjxLlPHAQAAAACwiqVJd2RkZJ7PYxs5cqRGjhyZ47omTZro+++/NyM0AAAAAACuWrF5ZBgAAAAAAMUNSTcAAAAAACYh6QYAAAAAwCTF4jndAABr+fpKK1deKMPz+Hr5auWAlc6yx6NTAwDchKQbAJAnLy+pWzero4CZvOxe6lb7GjrJdGoAgJswvRwAAAAAAJMw0g0AyFNqqrR4cWZ54EDJ29vaeFD4UtNTtfiXzJM88PqB8nZ4+EmmUwMA3ISkGwCQp5QUadiwzHLfvuQnniglPUXDPs08yX3r9fX8pJtODaC4iI4ussczDEMdOnSQw+HQmjVrXNbNnTtXjz/+uObMmaN77rknx/1jY2MVEhLifP3333+revXqql27tnbt2pVt+w0bNmjKlCnasWOHkpKSVKlSJd1888166623NHLkSC1atCjXWKtVq6aDBw/m+70VJqaXAwAAAACumM1m04IFC7Rlyxa9+eabzvoDBw5o0qRJmj17tipXrixJ2rt3r2JjY12WChUquLS3cOFC3XXXXYqLi9OWLVtc1v3666/q3LmzmjVrpo0bN+qXX37R7Nmz5ePjo/T0dM2aNculbUlasGCB8/WPP/5o8qeRO0a6AQAAAAAFUqVKFc2aNUsPPPCAOnbsqPDwcA0fPlwdO3bUPffco5iYGElShQoVFBwcnGs7hmFowYIFmjt3ripXrqz58+erRYsWzvVffvmlQkJC9MILLzjratSooc6dO0uS/P39FRQU5NJmcHCwy0i6VRjpBgAAAAAU2JAhQ9SuXTvde++9mjNnjnbt2uUy8p0f69evV2Jiotq3b69BgwZp6dKlSkhIcK4PCQlRbGysNm7cWNjhm46kGwAAAABwVebNm6ddu3Zp3LhxmjdvnsqXL++yvnLlyipRooRzqV+/vsv6+fPnq3///nI4HGrQoIGqV6+uZcuWOdf37dtXAwYMUJs2bRQaGqrevXtrzpw5iouLc8v7uxok3QAAAACAq1KhQgWNGjVKdevWVa9evbKt37Rpk3bs2OFcvvjiC+e6M2fO6OOPP9agQYOcdYMGDdL8+fOdrx0OhxYsWKC///5bL7zwgipVqqRp06apfv36zmu4iyqu6QYAAAAAXDUvLy95eeWcYkZEROR6TfeSJUuUlJTkcg23YRjKyMjQ77//rtq1azvrK1WqpHvuuUf33HOP/v3vf6t27dp64403NGXKlEJ9L4WJpBsAkCdfX+nDDy+U4Xl8vXz1YZ8PnWWPR6cGgCJj/vz5evjhhzV06FCX+tGjR+vtt9/Wc889l+N+pUuXVmhoqMu130URSTcAIE9eXpmPMobn8rJ7qW/9a+gk06kBwK2OHz+upKQkl7qyZctq9+7d2r59uxYvXqzrrrvOZf2AAQM0depUPfPMM5o/f7527Nih3r17q0aNGkpKStI777yj3bt3a/bs2e58K1eMpBsAAAAAYKo6depkq9u8ebMWL16sevXqZUu4Jal379564IEH9MUXX+jGG2/UN998o/vvv19Hjhxx3oxtxYoVatOmjTveQoHZDMMwrA7CanFxcQoKCtLZs2dVqlQpq8MBgCInLU365JPMcu/emYOEKB6iY6Jzro90rU/LSNMnezJPcu+6veVl9/CTTKcGUMQkJSXpwIEDioiIkJ+fn9Xh4H8ud17ym0fyDQMAyFNysnTXXZnl+HjyE0+UnJasuz7KPMnxUfHy8vHwk0ynBgC4CY8MAwAAAADAJCTdAAAAAACYhKQbAAAAAACTkHQDAAAAAGASkm4AAAAAAExC0g0AAAAAgEl4PgYAIE8+PtKCBRfK8Dw+Dh8t6LnAWfZ4dGoAgJuQdAMA8uTtLQ0danUUMJO3w1tDGw+1Ogz3oVMDANyE6eUAAAAAAJiEkW4AQJ7S0qQ1azLLnTpJXnx7eJy0jDSt2Zd5kjvV7CQvu4efZDo1gGIiJibGrceLjIws0H5//fWXJk+erNWrV+vkyZMKDQ1Vr1699PTTT+vcuXOKiIi47P4LFixQeHi42rZtq9OnTys4ONhlfXh4uMaNG6dx48Y5Xx86dChbO9OnT9djjz2mgwcPuhyzdOnSuv766/XMM8/olltuKdB7LCi+YQAAeUpOlm6/PbMcH09+4omS05J1+/uZJzk+Kl5ePh5+kunUAFBo/vzzT7Vs2VK1a9fW+++/r4iICO3evVsTJ07UqlWrtHnzZsXGxjq3f+mll7R69WqtXbvWWRcUFKQtW7Zc0XGnTp2qESNGuNSVLFnS5fXatWtVv359nTx5Us8++6xuv/12/f7776pYsWIB3mnB8A0DAAAAACiwMWPGyMfHR19++aX8/f0lSVWrVtUNN9ygGjVq6Mknn9Trr7/u3L5EiRLy8vJSSEjIVR23ZMmSebZRtmxZhYSEKCQkRI8//riWLl2qLVu2qEePHld17CvBNd0AAAAAgAI5deqU1qxZo9GjRzsT7iwhISEaOHCgPvjgAxmGYVGEmc6fP6933nlHkuTj5qdWMNINAAAAACiQP/74Q4ZhqG7dujmur1u3rk6fPq0TJ06oQoUK+WqzcuXK2eoSExOz1T366KN68sknXepWrVrlcs32zTffLLvdrsTERBmGoaZNm6pdu3b5iqOwkHQDAAAAAK5KYY5kb9q0Kdu12Tnd4G3ixIkaesnjHytVquTy+oMPPtB1112nXbt2adKkSVq4cKG8vb0LLdb8IOkGAAAAABRIzZo1ZbPZtGfPHvXu3Tvb+j179qh06dIqX758vtuMiIjIdvdyrxxueFmuXDnVrFnzsm1VqVJFtWrVUq1atZSWlqbevXtr165d8vX1zXc8V4trugEAAAAABVK2bFl16NBBc+fO1fnz513WHT16VIsXL1a/fv1ks9ksivCCPn36yMvLS3PnznXrcRnpBgDkycdHmjPnQhmex8fhozld5jjLHo9ODQCFZs6cObr55pvVqVMnPfPMMy6PDKtUqZKeffZZU4577tw5HT161KUuICBApUqVynF7m82mhx56SNHR0Ro1apQCAgJMietSjHQDAPLk7S2NGZO5uPkyKLiJt8NbY24cozE3jpG34xo4yXRqACg0tWrV0tatW1W9enXdddddqlGjhkaOHKm2bdtq8+bNKlOmjCnHffrppxUaGuqyTJo06bL7DBkyRKmpqZqT9cOrG9gMq+/dXgTExcUpKChIZ8+ezfVXEQAAiqPomOic6yNzrgcAWCMpKUkHDhxQRESE/Pz8rA4H/3O585LfPJLp5QCAPKWnS5s2ZZZvuUVyOKyNB4UvPSNdmw5nnuRbqt4ih93DTzKdGgDgJiTdAIA8JSVJbdtmluPjpcBAa+NB4UtKS1LbRZknOT4qXoE+Hn6S6dQAADch6QYA4Bp06bTzlPQUawIBAMDDcSM1AAAAAABMQtINAAAAAEUE97kuWgrjfJB0AwAAAIDFHP+7oWNKCpf7FCWJiYmSJO+reLwk13QDAAAAgMW8vLwUEBCgEydOyNvbW3Y746NWMgxDiYmJOn78uIKDg50/ihQESTcAAAAAWMxmsyk0NFQHDhzQoUOHrA4H/xMcHKyQkJCraoOkGwCQJ29v6YUXLpTheRw2h15on3mSvR3XwEmmUwMognx8fFSrVi2mmBcR3t7eVzXCnYWkGwCQJx8faeJEq6OAmRx2hya2uoZOMp0aQBFlt9vl5+dndRgoRFwoAAAAAACASRjpBgDkKT1d2r49s9ykiVQIM61QxGQYGfrxnx8lSU1Cm8hh9/CTTKcGALgJSTcAIE9JSdKNN2aW4+OlwEBr40HhS8tI043/yTzJ8VHxCvTx8JNMpwYAuAnTywEAAAAAMAlJNwAAAAAAJiHpBgAAAADAJJYm3Rs3blT37t0VFhYmm82mFStWuKwfOnSobDaby9K5c2eXbU6dOqWBAweqVKlSCg4O1vDhwxUfH+/GdwEAAAAAQM4sTboTEhLUqFEjvfbaa7lu07lzZ8XGxjqX999/32X9wIEDtXv3bn311VdauXKlNm7cqJEjR5odOgAAAAAAebL07uVdunRRly5dLruNr6+vQkJCcly3Z88erV69Wj/++KOaNWsmSZo9e7a6du2ql156SWFhYYUeMwAAnu7ZTc/Kx+HjfB0dGW1dMAAAFHNF/pFhMTExqlChgkqXLq3bbrtNzzzzjMqWLStJ2rx5s4KDg50JtyS1b99edrtdW7ZsUe/eva0KGwA8ire3NHnyhTI8j8PmUJtqbZxlj0enBgC4SZFOujt37qw77rhDERER2r9/vx5//HF16dJFmzdvlsPh0NGjR1WhQgWXfby8vFSmTBkdPXo013aTk5OVnJzsfB0XF2faewAAT+DjI0VHWx0FzOSwOxQZHml1GO5DpwYAuEmRTrr79+/vLF9//fVq2LChatSooZiYGLVr167A7U6fPl1TpkwpjBABAAAAAMhVsXpkWPXq1VWuXDnt27dPkhQSEqLjx4+7bJOWlqZTp07leh24JEVFRens2bPO5a+//jI1bgAo7jIypN27M5eMDKujgRkMw9DxhOM6nnBchmFYHY756NQAADcp0iPdl/r777/13//+V6GhoZKkli1b6syZM9q2bZuaNm0qSVq3bp0yMjLUokWLXNvx9fWVr6+vW2IGAE9w/rzUoEFmOT5eCgy0Nh4UvtSMVL2+9XVJUlTrKJcbqXkkOjUAwE0sTbrj4+Odo9aSdODAAe3YsUNlypRRmTJlNGXKFN15550KCQnR/v37NWnSJNWsWVOdOnWSJNWtW1edO3fWiBEj9MYbbyg1NVUPPPCA+vfvz53LAQAAAACWs3R6+datW3XDDTfohhtukCRNmDBBN9xwg55++mk5HA79/PPP6tGjh2rXrq3hw4eradOm2rRpk8so9eLFi3XdddepXbt26tq1q1q3bq158+ZZ9ZYAAAAAAHCydKQ7MjLysteNrVmzJs82ypQpoyVLlhRmWAAAAAAAFIpidSM1AAAAAACKE5JuAAAAAABMQtINAAAAAIBJitUjwwAA1vD2lh555EIZnsdhc6hl5ZbOssejUwMA3ISkGwCQJx8f6cUXrY4CZnLYHepYo6PVYbgPnRoA4CZMLwcAAAAAwCSMdAMA8pSRIR0+nFmuWlWy85OtxzEMQ2eTz0qSgnyDZLPZLI7IZHRqAICbkHQDAPJ0/rwUEZFZjo+XAgOtjQeFLzUjVbO2zJIkRbWOko/Dx+KITEanBgC4CT/rAgAAAABgEpJuAAAAAABMQtINAAAAAIBJSLoBAAAAADAJSTcAAAAAACYh6QYAAAAAwCQ8MgwAkCcvL2n06AtleB67za5mYc2cZY9HpwYAuAnfMgCAPPn6Sq+9ZnUUMJOX3UvdanWzOgz3oVMDANzkGvgpGwAAAAAAazDSDQDIk2FIJ09mlsuVk2w2a+NB4TMMQ4mpiZKkAO8A2Tz9JNOpAQBuQtINAMhTYqJUoUJmOT5eCgy0Nh4UvtSMVL20+SVJUlTrKPk4fCyOyGR0agCAmzC9HAAAAAAAk5B0AwAAAABgEpJuAAAAAABMQtINAAAAAIBJSLoBAAAAADAJSTcAAAAAACbhkWEAgDx5eUlDhlwow/PYbXY1qtjIWfZ4dGoAgJvwLQMAyJOvr7RwodVRwExedi/1uq6X1WG4D50aAOAm18BP2QAAAAAAWIORbgBAngxDSkzMLAcESDabtfGg8BmGodSMVEmSt91bNk8/yXRqAICbMNINAMhTYqJUokTmkpWnwLOkZqRq+jfTNf2b6c7k26PRqQEAbkLSDQAAAACASUi6AQAAAAAwCUk3AAAAAAAmIekGAAAAAMAkJN0AAAAAAJiEpBsAAAAAAJPwnG4AQJ4cDqlPnwtleB67za565eo5yx6PTg0AcBOSbgBAnvz8pGXLrI4CZvKye6lv/b5Wh+E+dGoAgJtcAz9lAwAAAABgDZJuAAAAAABMQtINAMhTQoJks2UuCQlWRwMzpKSnaMqGKZqyYYpS0lOsDsd8dGoAgJuQdAMAAAAAYBKSbgAAAAAATELSDQAAAACASUi6AQAAAAAwCUk3AAAAAAAmIekGAAAAAMAkXlYHAAAo+hwOqWvXC2V4HrvNrlplajnLHo9ODQBwE5JuAECe/Pykzz+3OgqYycvupbuvv9vqMNyHTg0AcJNr4KdsAAAAAACsQdINAAAAAIBJSLoBAHlKSJACAzOXhASro4EZUtJTNG3TNE3bNE0p6SlWh2M+OjUAwE24phsAkC+JiVZHALOlZqRaHYJ70akBAG7ASDcAAAAAACYh6QYAAAAAwCSWJt0bN25U9+7dFRYWJpvNphUrVjjXpaam6tFHH9X111+vwMBAhYWFafDgwTpy5IhLG+Hh4bLZbC7Lc8895+Z3AgCA9aJjorMtAADAWpYm3QkJCWrUqJFee+21bOsSExO1fft2PfXUU9q+fbs+/vhj7d27Vz169Mi27dSpUxUbG+tcHnzwQXeEDwAAAADAZVl6I7UuXbqoS5cuOa4LCgrSV1995VI3Z84c3XjjjTp8+LCqVq3qrC9ZsqRCQkJMjRUAAAAAgCtVrK7pPnv2rGw2m4KDg13qn3vuOZUtW1Y33HCDXnzxRaWlpVkTIAB4KLtdatMmc7EXq28O5JdNNlULqqZqQdVkk83qcMxHpwYAuEmxeWRYUlKSHn30UQ0YMEClSpVy1j/00ENq0qSJypQpo++++05RUVGKjY3Vyy+/nGtbycnJSk5Odr6Oi4szNXYAKO78/aWYGKujgJm8Hd4a2nio1WG4D50aAOAmxSLpTk1N1V133SXDMPT666+7rJswYYKz3LBhQ/n4+GjUqFGaPn26fH19c2xv+vTpmjJliqkxAwAAAABQ5OdTZSXchw4d0ldffeUyyp2TFi1aKC0tTQcPHsx1m6ioKJ09e9a5/PXXX4UcNQAAAAAARXykOyvh/uOPP7R+/XqVLVs2z3127Nghu92uChUq5LqNr69vrqPgAIDsEhKk8PDM8sGDUmCgldHADCnpKZq1ZZYkaWyLsfJx+Fgckcno1AAAN7E06Y6Pj9e+ffucrw8cOKAdO3aoTJkyCg0NVZ8+fbR9+3atXLlS6enpOnr0qCSpTJky8vHx0ebNm7Vlyxa1bdtWJUuW1ObNmzV+/HgNGjRIpUuXtuptAYBHOnnS6ghgtsTURKtDcC86NQDADSxNurdu3aq2bds6X2ddnz1kyBBFR0frs88+kyQ1btzYZb/169crMjJSvr6+Wrp0qaKjo5WcnKyIiAiNHz/e5TpvAAAAAACsYmnSHRkZKcMwcl1/uXWS1KRJE33//feFHRYAAAAAAIWiyN9IDQAAAACA4oqkGwAAAAAAk5B0AwAAAABgkiL9yDAAQNFgt0vNml0ow/PYZFNYyTBn2ePRqQEAbkLSDQDIk7+/9OOPVkcBM3k7vDWiyQirw3AfOjUAwE34aRcAAAAAAJOQdAMAAAAAYBKmlwMA8pSYKNWrl1n+9VcpIMDaeFD4UtNT9dqPr0mSxjQfI2+Ht3NddEx0tu2jI7PXFSt0agCAm5B0AwDyZBjSoUMXyvA8hgydTT7rLHs8OjUAwE2YXg4AAAAAgElIugEAAAAAMAlJNwAAAAAAJiHpBgAAAADAJCTdAAAAAACYhLuXAwDyZLNdeLqSzWZtLDCHTTaVDyjvLHs8OjUAwE1IugEAeQoIkHbvtjoKmMnb4a3RzUdbHYb70KkBAG7C9HIAAAAAAExC0g0AAAAAgElIugEAeUpMlOrXz1wSE62OBmZITU/V3B/nau6Pc5Wanmp1OOajUwMA3IRrugEAeTIM6ddfL5TheQwZOpF4wln2eHRqAICbMNINAAAAAIBJSLoBAAAAADAJSTcAAAAAACYh6QYAAAAAwCQk3QAAAAAAmIS7lwMA8mSzSdWqXSjD89hkU5BvkLPs8ejUAAA3IekGAOQpIEA6eNDqKGAmb4e3xt00zuow3IdODQBwE6aXAwAAAABgEpJuAAAAAABMQtINAMjT+fNS8+aZy/nzVkcDM6Smp+qt7W/pre1vKTU91epwzEenBgC4Cdd0AwDylJEhbd16oQzPY8jQkXNHnGWPR6cGALgJI90AAAAAAJiEpBsAAAAAAJOQdAMAAAAAYBKSbgAAAAAATELSDQAAAACASbh7OQAgX8qVszoCmC3AO8DqENyLTg0AcAOSbgBAngIDpRMnrI4CZvJx+GjizROtDsN96NQAADcp0PTyP//8s7DjAAAAAADA4xQo6a5Zs6batm2r9957T0lJSYUdEwAAAAAAHqFASff27dvVsGFDTZgwQSEhIRo1apR++OGHwo4NAFBEnD8vRUZmLufPWx0NzJCanqqFOxZq4Y6FSk1PtToc89GpAQBuUqCku3Hjxpo1a5aOHDmit99+W7GxsWrdurUaNGigl19+WSe4RgoAPEpGhrRhQ+aSkWF1NDCDIUOHzh7SobOHZMiwOhzz0akBAG5yVY8M8/Ly0h133KFly5bp+eef1759+/TII4+oSpUqGjx4sGJjYwsrTgAAAAAAip2rSrq3bt2q0aNHKzQ0VC+//LIeeeQR7d+/X1999ZWOHDminj17FlacAAAAAAAUOwV6ZNjLL7+sBQsWaO/everataveeecdde3aVXZ7Zg4fERGhhQsXKjw8vDBjBQAAAACgWClQ0v3666/r3nvv1dChQxUaGprjNhUqVND8+fOvKjgAAAAAAIqzAiXdf/zxR57b+Pj4aMiQIQVpHgAAAAAAj1CgpHvBggUqUaKE+vbt61K/bNkyJSYmkmwDgAcKCLA6ApjN2+6d722jY6Jzro/Mub5IolMDANygQDdSmz59usqVK5etvkKFCpo2bdpVBwUAKFoCA6WEhMwlMNDqaGAGH4ePHr/lcT1+y+PycfhYHY756NQAADcp0Ej34cOHFRERka2+WrVqOnz48FUHBQAAcpfbKDMAACh6CjTSXaFCBf3888/Z6nfu3KmyZctedVAAAAAAAHiCAiXdAwYM0EMPPaT169crPT1d6enpWrduncaOHav+/fsXdowAAIslJUndumUuSUlWRwMzpGWkackvS7TklyVKy0izOhzz0akBAG5SoOnl//73v3Xw4EG1a9dOXl6ZTWRkZGjw4MFc0w0AHig9XfriiwtleJ4MI0N/nPrDWfZ4dGoAgJsUKOn28fHRBx98oH//+9/auXOn/P39df3116tatWqFHR8AAAAAAMVWgZLuLLVr11bt2rULKxYAAAAAADxKga7pTk9P1/z583X33Xerffv2uu2221yW/Nq4caO6d++usLAw2Ww2rVixwmW9YRh6+umnFRoaKn9/f7Vv315//PGHyzanTp3SwIEDVapUKQUHB2v48OGKj48vyNsCAAAAAKBQFSjpHjt2rMaOHav09HQ1aNBAjRo1clnyKyEhQY0aNdJrr72W4/oXXnhBr776qt544w1t2bJFgYGB6tSpk5IuuuHJwIEDtXv3bn311VdauXKlNm7cqJEjRxbkbQEAAAAAUKgKNL186dKl+vDDD9W1a9erOniXLl3UpUuXHNcZhqGZM2fqySefVM+ePSVJ77zzjipWrKgVK1aof//+2rNnj1avXq0ff/xRzZo1kyTNnj1bXbt21UsvvaSwsLCrig8AAAAAgKtRoJFuHx8f1axZs7BjcXHgwAEdPXpU7du3d9YFBQWpRYsW2rx5syRp8+bNCg4OdibcktS+fXvZ7XZt2bLF1PgAAAAAAMhLgUa6H374Yc2aNUtz5syRzWYr7JgkSUePHpUkVaxY0aW+YsWKznVHjx5VhQoVXNZ7eXmpTJkyzm1ykpycrOTkZOfruLi4wgobADxSYKBkGFZHATP5OHw0uc1kq8NwHzo1AMBNCpR0f/PNN1q/fr1WrVql+vXry9vb22X9xx9/XCjBmWX69OmaMmWK1WEAAAAAADxcgZLu4OBg9e7du7BjcRESEiJJOnbsmEJDQ531x44dU+PGjZ3bHD9+3GW/tLQ0nTp1yrl/TqKiojRhwgTn67i4OFWpUqUQowcAAAAAoIBJ94IFCwo7jmwiIiIUEhKir7/+2plkx8XFacuWLfrXv/4lSWrZsqXOnDmjbdu2qWnTppKkdevWKSMjQy1atMi1bV9fX/n6+pr+HgDAUyQlSffck1l+913Jz8/aeFD40jLS9MmeTyRJvev2lpe9QP9FKD7o1AAANynwN2paWppiYmK0f/9+3X333SpZsqSOHDmiUqVKqUSJEvlqIz4+Xvv27XO+PnDggHbs2KEyZcqoatWqGjdunJ555hnVqlVLEREReuqppxQWFqZevXpJkurWravOnTtrxIgReuONN5SamqoHHnhA/fv3587lAFCI0tOljz7KLC9caGkoMEmGkaFfT/4qSepp9LQ4GjegUwMA3KRASfehQ4fUuXNnHT58WMnJyerQoYNKliyp559/XsnJyXrjjTfy1c7WrVvVtm1b5+usKd9DhgzRwoULNWnSJCUkJGjkyJE6c+aMWrdurdWrV8vvol+jFy9erAceeEDt2rWT3W7XnXfeqVdffbUgbwsAAAAAgEJVoKR77NixatasmXbu3KmyZcs663v37q0RI0bku53IyEgZl7lzqM1m09SpUzV16tRctylTpoyWLFmS72MCAAAAAOAuBUq6N23apO+++04+Pj4u9eHh4frnn38KJTAAAAAAAIo7e0F2ysjIUHp6erb6v//+WyVLlrzqoAAAAAAA8AQFSro7duyomTNnOl/bbDbFx8dr8uTJ6tq1a2HFBgAAAABAsVag6eUzZsxQp06dVK9ePSUlJenuu+/WH3/8oXLlyun9998v7BgBAAAAACiWCpR0V65cWTt37tTSpUv1888/Kz4+XsOHD9fAgQPl7+9f2DECACwWECDFx18ow/N4270V1TrKWfZ4dGoAgJsU+DndXl5eGjRoUGHGAgAoomw2KTDQ6ihgJpvNJh+HT94bego6NQDATQqUdL/zzjuXXT948OACBQMAAAAAgCcp8HO6L5aamqrExET5+PgoICCApBsAPExysjRqVGb5zTclX19r40HhS8tI08rfV0qSbq99u7zsBZ4MVzzQqQEAblKgu5efPn3aZYmPj9fevXvVunVrbqQGAB4oLU1atChzSUuzOhqYIcPI0M5jO7Xz2E5lGBlWh2M+OjUAwE0KlHTnpFatWnruueeyjYIDAAAAAHCtKrSkW8q8udqRI0cKs0kAAAAAAIqtAl2w9dlnn7m8NgxDsbGxmjNnjlq1alUogQEAAAAAUNwVKOnu1auXy2ubzaby5cvrtttu04wZMwojLgAAAAAAir0CJd0ZGdfADVYAAAAAALhKhXpNNwAAAAAAuKBAI90TJkzI97Yvv/xyQQ4BAChCAgKk48cvlOF5vO3eeqTlI86yx6NTAwDcpEBJ908//aSffvpJqampqlOnjiTp999/l8PhUJMmTZzb2Wy2wokSAGApm00qX97qKGAmm82mQJ9Aq8NwHzo1AMBNCpR0d+/eXSVLltSiRYtUunRpSdLp06c1bNgw3XLLLXr44YcLNUgAAAAAAIqjAl3TPWPGDE2fPt2ZcEtS6dKl9cwzz3D3cgDwQMnJ0pgxmUtystXRwAxpGWn6/I/P9fkfnystI83qcMxHpwYAuEmBku64uDidOHEiW/2JEyd07ty5qw4KAFC0pKVJc+dmLmnXQD52LcowMrT1yFZtPbJVGcY18JQSOjUAwE0KlHT37t1bw4YN08cff6y///5bf//9t5YvX67hw4frjjvuKOwYAQAAAAAolgp0Tfcbb7yhRx55RHfffbdSU1MzG/Ly0vDhw/Xiiy8WaoAAAAAAABRXBUq6AwICNHfuXL344ovav3+/JKlGjRoKDLyG7noKAAAAAEAeCjS9PEtsbKxiY2NVq1YtBQYGyjCMwooLAAAAAIBir0BJ93//+1+1a9dOtWvXVteuXRUbGytJGj58OI8LAwAAAADgfwqUdI8fP17e3t46fPiwAgICnPX9+vXT6tWrCy04AAAAAACKswJd0/3ll19qzZo1qly5skt9rVq1dOjQoUIJDABQdPj7SwcOXCjD83jbvTW2xVhn2ePRqQEAblKgpDshIcFlhDvLqVOn5Ovre9VBAQCKFrtdCg+3OgqYyWazKdgv2Oow3IdODQBwkwJNL7/lllv0zjvvOF/bbDZlZGTohRdeUNu2bQstOAAAAAAAirMCjXS/8MILateunbZu3aqUlBRNmjRJu3fv1qlTp/Ttt98WdowAAIulpEhPPJFZfvZZycfH2nhQ+NIz0vX1ga8lSe0i2slhd1gckcno1AAANynQSHeDBg30+++/q3Xr1urZs6cSEhJ0xx136KefflKNGjUKO0YAgMVSU6WXXspcUlOtjgZmSDfStfnvzdr892alG+lWh2M+OjUAwE2ueKQ7NTVVnTt31htvvKEnsn4hBgAAAAAA2VzxSLe3t7d+/vlnM2IBAAAAAMCjFGh6+aBBgzR//vzCjgUAAAAAAI9SoBuppaWl6e2339batWvVtGlTBQYGuqx/+eWXCyU4AAAAAACKsytKuv/880+Fh4dr165datKkiSTp999/d9nGZrMVXnQAAAAAABRjV5R016pVS7GxsVq/fr0kqV+/fnr11VdVsWJFU4IDAAAAAKA4u6Kk2zAMl9erVq1SQkJCoQYEACh6/P2lXbsulOF5vO3e+lezfznLHo9ODQBwkwJd053l0iQcAOCZ7Hapfn2ro4CZbDabKgRWsDoM96FTAwDc5IruXm6z2bJds8013AAAAAAA5OyKp5cPHTpUvr6+kqSkpCTdf//92e5e/vHHHxdehAAAy6WkSNOmZZYff1zy8bE2HhS+9Ix0bTq8SZJ0S9Vb5LA7LI7IZHRqAICbXFHSPWTIEJfXgwYNKtRgAABFU2qqNGVKZnniRPITT5RupGvDoQ2SpJur3CyHPDzpplMDANzkipLuBQsWmBUHAAAAAAAe54qu6QYAAAAAAPlH0g0AAAAAgElIugEAAAAAMAlJNwAAAAAAJiHpBgAAAADAJFd093IAwLXJz0/64YcLZXgeL7uX7rvhPmfZ49GpAQBucg18qwIArpbDITVvbnUUMJPdZlelUpWsDsN96NQAADdhejkAAAAAACZhpBsAkKeUFGnWrMzy2LGSj4+18aDwpWek6/t/vpck3VTpJjnsDosjMhmdGgDgJiTdAIA8paZKkyZllkePJj/xROlGutb+uVaS1DysuRzy8KSbTg0AcBOmlwMAAAAAYBKSbgAAAAAATFLkk+7w8HDZbLZsy5gxYyRJkZGR2dbdf//9FkcNAAAAAEAxuKb7xx9/VHp6uvP1rl271KFDB/Xt29dZN2LECE2dOtX5OiAgwK0xAgAAAACQkyKfdJcvX97l9XPPPacaNWqoTZs2zrqAgACFhIS4OzQAAHCJ6Jjo7HWR2esAALhWFPnp5RdLSUnRe++9p3vvvVc2m81Zv3jxYpUrV04NGjRQVFSUEhMTLYwSAAAAAIBMRX6k+2IrVqzQmTNnNHToUGfd3XffrWrVqiksLEw///yzHn30Ue3du1cff/xxru0kJycrOTnZ+TouLs7MsAGg2PPzk9avv1CG5/Gye2lIoyHOssejUwMA3KRYfavOnz9fXbp0UVhYmLNu5MiRzvL111+v0NBQtWvXTvv371eNGjVybGf69OmaMmWK6fECgKdwOKTISKujgJnsNrvCg8OtDsN96NQAADcpNtPLDx06pLVr1+q+++677HYtWrSQJO3bty/XbaKionT27Fnn8tdffxVqrAAAAAAASMVopHvBggWqUKGCunXrdtntduzYIUkKDQ3NdRtfX1/5+voWZngA4NFSU6V58zLLI0dK3t7WxoPCl56Rrm2x2yRJTUObymF3WByRyejUAAA3KRZJd0ZGhhYsWKAhQ4bIy+tCyPv379eSJUvUtWtXlS1bVj///LPGjx+vW2+9VQ0bNrQwYgDwLCkp0gMPZJaHDiU/8UTpRrpW7VslSWoc0lgOeXjSTacGALhJsUi6165dq8OHD+vee+91qffx8dHatWs1c+ZMJSQkqEqVKrrzzjv15JNPWhQpAAAAAAAXFIuku2PHjjIMI1t9lSpVtGHDBgsiAgAAAAAgb8XmRmoAAAAAABQ3JN0AAAAAAJiEpBsAAAAAAJOQdAMAAAAAYJJicSM1AIC1fH2llSsvlOF5vOxeGtBggLPs8ejUAAA3uQa+VQEAV8vLS+rWzeooYCa7za7aZWtbHYb70KkBAG7C9HIAAAAAAEzCSDcAIE+pqdLixZnlgQMlb29r40HhS89I1y/Hf5EkXV/hejnsDosjMhmdGgDgJiTdAIA8paRIw4Zllvv2JT/xROlGuj7d+6kkqV75enLIw5NuOjUAwE2YXg4AAAAAgElIugEAAAAAMAlJNwAAAAAAJiHpBgAAAADAJCTdAAAAAACYhKQbAAAAAACT8MgwAECefH2lDz+8UIb7RMdEu+U4XnYv9anXx1n2eHRqAICbXAPfqgCAq+XllfkoY3guu82u+uXrWx2G+9CpAQBuwvRyAAAAAABMwkg3ACBPaWnSJ59klnv3zhwkhGfJMDK05+QeSVLdcnVlt3n47/J0agCAm/ANAwDIU3KydNddmeX4ePITT5SWkaaPfv1IkhTVOko+Dh+LIzIZnRoA4CYe/jM2AAAAAADWIekGAAAAAMAkJN0AAAAAAJiEpBsAAAAAAJOQdAMAAAAAYBKSbgAAAAAATMLzMQAAefLxkRYsuFCG53HYHOpZp6ez7PHo1AAANyHpBgDkydtbGjrU6ihgJofdocYhja0Ow33o1AAAN2F6OQAAAAAAJmGkGwCQp7Q0ac2azHKnTpIX3x4eJ8PI0L5T+yRJNcvUlN3m4b/L06kBAG7CNwwAIE/JydLtt2eW4+PJTzxRWkaa3t/1viQpqnWUfBwefp0znRoA4CYe/jM2AAAAAADWIekGAAAAAMAkJN0AAAAAAJiEpBsAAAAAAJOQdAMAAAAAYBKSbgAAAAAATMLzMQAAefLxkebMuVCG53HYHOpSs4uz7PHo1AAANyHpBgDkydtbGjPG6ihgJofdoRsr3Wh1GO5DpwYAuAnTywEAAAAAMAkj3QCAPKWnS5s2ZZZvuUVyXAOzj681GUaGDp89LEmqGlRVdpuH/y5PpwYAuAlJNwAgT0lJUtu2meX4eCkw0Np4UPjSMtK0aOciSVJU6yj5ODz8Omc6NQDATTz8Z2wAAAAAAKxD0g0AAAAAgEmYXg4AQBEQHRNtdQgAAMAEjHQDAAAAAGASkm4AAAAAAExC0g0AAAAAgEm4phsAkCdvb+mFFy6U4XkcNofaV2/vLHs8OjUAwE1IugEAefLxkSZOtDoKmMlhd6hVlVZWh+E+dGoAgJswvRwAAAAAAJMw0g0AyFN6urR9e2a5SRPJcQ3MPr7WZBgZij0XK0kKLRkqu83Df5enUwMA3ISkGwCQp6Qk6cYbM8vx8VJgoLXxoPClZaTpPz/9R5IU1TpKPg4fiyMyGZ0aAOAmHv4zNgAAAAAA1iHpBgAAAADAJEU66Y6OjpbNZnNZrrvuOuf6pKQkjRkzRmXLllWJEiV055136tixYxZGDAAAAADABUU66Zak+vXrKzY21rl88803znXjx4/X//3f/2nZsmXasGGDjhw5ojvuuMPCaAEAAAAAuKDI30jNy8tLISEh2erPnj2r+fPna8mSJbrtttskSQsWLFDdunX1/fff66abbnJ3qAAAAAAAuCjyI91//PGHwsLCVL16dQ0cOFCHDx+WJG3btk2pqalq3769c9vrrrtOVatW1ebNm60KFwAAAAAApyI90t2iRQstXLhQderUUWxsrKZMmaJbbrlFu3bt0tGjR+Xj46Pg4GCXfSpWrKijR49ett3k5GQlJyc7X8fFxZkRPgB4DG9vafLkC2V4HofNoTbV2jjLHo9ODQBwkyKddHfp0sVZbtiwoVq0aKFq1arpww8/lL+/f4HbnT59uqZMmVIYIQLANcHHR4qOtjoKmMlhdygyPNLqMNyHTg0AcJMiP738YsHBwapdu7b27dunkJAQpaSk6MyZMy7bHDt2LMdrwC8WFRWls2fPOpe//vrLxKgBAAAAANeqYpV0x8fHa//+/QoNDVXTpk3l7e2tr7/+2rl+7969Onz4sFq2bHnZdnx9fVWqVCmXBQCQu4wMaffuzCUjw+poYAbDMHQ84biOJxyXYRhWh2M+OjUAwE2K9PTyRx55RN27d1e1atV05MgRTZ48WQ6HQwMGDFBQUJCGDx+uCRMmqEyZMipVqpQefPBBtWzZkjuXA0AhO39eatAgsxwfLwUGWhsPCl9qRqpe3/q6JCmqdZR8HD4WR2QyOjUAwE2KdNL9999/a8CAAfrvf/+r8uXLq3Xr1vr+++9Vvnx5SdIrr7wiu92uO++8U8nJyerUqZPmzp1rcdQAAAAAAGQq0kn30qVLL7vez89Pr732ml577TU3RQQAAK5UdEx0zvWROdcDAOBJitU13QAAAAAAFCck3QAAAAAAmISkGwAAAAAAk5B0AwAAAABgkiJ9IzUAQNHg7S098siFMjyPw+ZQy8otnWWPR6cGALgJSTcAIE8+PtKLL1odBczksDvUsUZHq8NwHzo1AMBNmF4OAAAAAIBJGOkGAOQpI0M6fDizXLWqZOcnW49jGIbOJp+VJAX5Bslms1kckcno1AAANyHpBgDk6fx5KSIisxwfLwUGWhsPCl9qRqpmbZklSYpqHSUfh4/FEZmMTg0AcBN+1gUAAAAAwCQk3QAAAAAAmISkGwAAAAAAk5B0AwAAAABgEpJuAAAAAABMQtINAAAAAIBJeGQYACBPXl7S6NEXyvA8dptdzcKaOcsej04NAHATvmUAAHny9ZVee83qKGAmL7uXutXqZnUY7kOnBgC4yTXwUzYAAAAAANZgpBsAkCfDkE6ezCyXKyfZbNbGg8JnGIYSUxMlSQHeAbJ5+kmmUwMA3ISkGwCQp8REqUKFzHJ8vBQYaG08KHypGal6afNLkqSo1lHycfhYHJHJ6NQAADdhejkAAAAAACYh6QYAAAAAwCQk3QAAAAAAmISkGwAAAAAAk5B0AwAAAABgEpJuAAAAAABMwiPDAAB58vKShgy5UIbnsdvsalSxkbPs8ejUAAA34VsGAJAnX19p4UKro4CZvOxe6nVdL6vDcB86NQDATUi6AQBwo+iYaKtDAAAAbkTSDQDIk2FIiYmZ5YAAyWazNh4UPsMwlJqRKknytnvL5uknmU4NAHCTa+CiLQDA1UpMlEqUyFyy8hR4ltSMVE3/ZrqmfzPdmXx7NDo1AMBNSLoBAAAAADAJSTcAAAAAACbhmm4AAGCJ3G4qFx2Zcz0AAMURI90AAAAAAJiEpBsAAAAAAJMwvRwAAJPwTG4AAEDSDQDIk8Mh9elzoQzPY7fZVa9cPWfZ49GpAQBuQtINAMiTn5+0bJnVUcBMXnYv9a3f1+ow3IdODQBwk2vgp2wAAAAAAKxB0g0AAAAAgElIugEAeUpIkGy2zCUhwepoYIaU9BRN2TBFUzZMUUp6itXhmI9ODQBwE5JuAAAAAABMQtINAAAAAIBJSLoBAAAAADAJSTcAAAAAACYh6QYAAAAAwCQk3QAAAAAAmMTL6gAAAEWfwyF17XqhDM9jt9lVq0wtZ9nj0akBAG5C0g0AyJOfn/T551ZHATN52b109/V3Wx2G+9CpAQBucg38lA0AAAAAgDVIugEAAAAAMAlJNwAgTwkJUmBg5pKQYHU0MENKeoqmbZqmaZumKSU9xepwzEenBgC4Cdd0AwDyJTHR6ghgttSMVKtDcC86NQDADRjpBgAAAADAJEU66Z4+fbqaN2+ukiVLqkKFCurVq5f27t3rsk1kZKRsNpvLcv/991sUMQAAAAAAFxTppHvDhg0aM2aMvv/+e3311VdKTU1Vx44dlXDJtVcjRoxQbGysc3nhhRcsihgAAAAAgAuK9DXdq1evdnm9cOFCVahQQdu2bdOtt97qrA8ICFBISIi7wwMAAAAA4LKK9Ej3pc6ePStJKlOmjEv94sWLVa5cOTVo0EBRUVFK5MYoAAAAAIAioEiPdF8sIyND48aNU6tWrdSgQQNn/d13361q1aopLCxMP//8sx599FHt3btXH3/8ca5tJScnKzk52fk6Li7O1NgBoLiz26U2bS6U4XlssqlaUDVn2ePRqQEAblJsku4xY8Zo165d+uabb1zqR44c6Sxff/31Cg0NVbt27bR//37VqFEjx7amT5+uKVOmmBovAHgSf38pJsbqKGAmb4e3hjYeanUY7kOnBgC4SbH4afeBBx7QypUrtX79elWuXPmy27Zo0UKStG/fvly3iYqK0tmzZ53LX3/9VajxAgAAAAAgFfGRbsMw9OCDD+qTTz5RTEyMIiIi8txnx44dkqTQ0NBct/H19ZWvr29hhQkAuMZFx0RbHQIAACiiinTSPWbMGC1ZskSffvqpSpYsqaNHj0qSgoKC5O/vr/3792vJkiXq2rWrypYtq59//lnjx4/XrbfeqoYNG1ocPQB4joQEKTw8s3zwoBQYaGU0MENKeopmbZklSRrbYqx8HD4WR2QyOjUAwE2KdNL9+uuvS5IiIyNd6hcsWKChQ4fKx8dHa9eu1cyZM5WQkKAqVarozjvv1JNPPmlBtADg2U6etDoCmC0x9Rp7+gedGgDgBkU66TYM47Lrq1Spog0bNrgpGgDAtSS3KePRkTnXAwAA5KRY3EgNAAAAAIDiiKQbAAAAAACTkHQDAAAAAGASkm4AAAAAAExSpG+kBgAoGux2qVmzC2V4HptsCisZ5ix7PDo1AMBNSLoBAHny95d+/NHqKGAmb4e3RjQZYXUY7kOnBgC4CT/tAgAAAABgEka6AQC4Ark9vxuFJ6fPmOejAwCKK0a6AQB5SkyUwsMzl8REq6OBGVLTUzXz+5ma+f1MpaanWh2O+ejUAAA3YaQbAJAnw5AOHbpQhucxZOhs8lln2ePRqQEAbsJINwAAAAAAJiHpBgAAAADAJEwvBwAARV5uN7DjBmsAgKKOkW4AAAAAAExC0g0AAAAAgEmYXg4AyJPNJtWrd6EMz2OTTeUDyjvLHo9ODQBwE5JuAECeAgKk3butjgJm8nZ4a3Tz0VaH4T50agCAm5B0AwAAzxcdnXMZAACTcU03AAAAAAAmIekGAOQpMVGqXz9zSUy0OhqYITU9VXN/nKu5P85Vanqq1eGYj04NAHATppcDAPJkGNKvv14ow/MYMnQi8YSz7PEu7tTPPCP5+GSWmXoOAChkjHQDAAAAAGASkm4AAAAAAEzC9HIAAFB8XW46OFPFAQBFACPdAAAAAACYhJFuAABQbMUcjMlWFxke6fY4AADIDUk3ACBPNptUrdqFMjyPTTYF+QY5yx6PTg0AcBOSbgBAngICpIMHrY4CZvJ2eGvcTeOsDsN9Lu7UXPsNADARSTcAAPBMJNMAgCKAG6kBAAAAAGASkm4AQJ7On5eaN89czp+3OhqYITU9VW9tf0tvbX9LqempVodjvos7deo18H4BAJZhejkAIE8ZGdLWrRfK8DyGDB05d8RZ9ngXd+oOHS7UXzwlnenpAIBCwEg3AAAAAAAmYaQbAAAUOZELY1xexwyNzHVdfmU90zsmJlre51P0RMFCy8SIOAAgnxjpBgAAAADAJCTdAAAAAACYhOnlAACPFB0Tnf9tI/O/La7cxdPBL54mbpasaeS5iVwYI3tqet4NXTptnGnkAIACIOkGAORLuXJWRwCzBXgHWB2CW6X4esvH4W11GAAAD0fSDQDIU2CgdOKE1VHATD4OH028eaLVYbhNhrdD3/W/WZHhkVaHAgDwcCTdAACgyCvoHcsLVX6mlzMlHQBwCW6kBgAAAACASRjpBgDk6fx5qUuXzPKqVZK/v7XxoPClpqdq8S+LJUkDrx8o73xc63y5Z2kXdfa0dDVc+4vkd1AaOFDy5tpuAIA5SLoBAHnKyJA2bLhQhucxZOjQ2UPOssczpOBjZyWdlYxr4P0CACzD9HIAAAAAAEzCSDcA4IpdfG8o7hOFa0Vuz/+OdGsUAIDihpFuAAAAAABMQtINAAAAAIBJmF4OAB4mt6nfnjolPDom2uoQPJ67npFt1bO4Nx7apAxvR8EbyO8flKf+EQIeJKfvlOjI7HXAlSDpBgDkS0CA1RHAbN72a+uxWeleTPgDAJiPpBsArnGXDrjlNAAXGCglJLgjmssza1Sb0XLJx+Gjx295/KrauHikOrdndls1mn2pDG+HNg28xfwDFWREmxFxwCm3f5+L6uhzcYsX7sFPvAAAAAAAmISkGwAAAAAAkzC9HAAscLkZo1bPJs1pZmtSknTnnZnl5cuvfH8UjvxM376iNiL/10a4lJaRpg93fyhJuqv+XfKye/Z/EezpGaq/frckaXfb+spwuHkcIj/XdRS0Pf7wiiSmHRctV3JZEecOV8uzv1EBAIUiPV364osLZXieDCNDf5z6w1n2eBmGyv5zylnWVdy8HACAyyHphsfL16MfCnGEoCg8aqIoxFAUFPSXaasHjK70+IU9ah5zMObC/jGZ5YnNc24o5mCMIoe61kWGR+YZy8XlyKExLttFRsY424jMGok1w8GLyuHmHaaoORgcnmP9Le9t0vTwC2V/4+qy0KJywzR3uPhvJi+5/X3kKoc/nKx/2y7+jK+43csdJ4djZds8r++Uwh7Jd4OsELPO56X/NvE9+r+6Ivw5uHtEujBGy81SlM/Ttchjrul+7bXXFB4eLj8/P7Vo0UI//PCD1SEBAAAAAK5xHpF0f/DBB5owYYImT56s7du3q1GjRurUqZOOHz9udWgAAAAAgGuYR0wvf/nllzVixAgNGzZMkvTGG2/o888/19tvv63HHnvM4ugKj6fcxOFK3ke+ti3A3NpL243MmkqWw/S8K50OdOmUypihkVf1nmMWXogp21S3i1/m8l6z9o9cGKPIoTHZbsRU4M/9Svwvtv/NVFZMZHSBH117NdP+sj4Ll89NUrRc94383/rc4ow5GKOhOxY6Xy9sPDRzv//1n9ymUefnPV86TfXSqdsu6y7qr5GX9ptLPo9L242Oibmizyzlywt1zz4r+fhkli/+HJx2LLzwmQzNud3LvS/pQrzRF2136TnPdd+L/mZcj3lh//BxB3XmumBJ0pnOwbm2ldsUaedNzA5mNZiv0JztHdwRfuEzioxR+I6DLtuFn7nwOutYMQsjnduFN3bd/tJYc5pCntVmbtPLL3Y4qKrqnPknz+2uRlYcF7/XS+sujtVlux0X1efwWRRluU5Fz+HGeJe+z5jLfB/FHIxxbr/won8noyOjs/3jk1sMLt+B/9sn8mCMy037svpZ1t9Zrs9jv+gYkeGR2f4xzOl7Jre/Xcn179f5vRYeedl/Y13W5eO7POZg7sfPjbunDJslx/dx8KJyeB7b5tZuLt8zhTFtvShP7S4KrvY9c2lk4Sr2I90pKSnatm2b2rdv76yz2+1q3769Nm/ebGFkAAAAAIBrXbEf6T558qTS09NVsWJFl/qKFSvqt99+y3Gf5ORkJScnO1+fPXtWkhQXF2deoIUgOSE5x/qiHvelruR95Gvb5Jy3+d+GubZxsYSUtMzNk5OlS+LIz/45tZWf/fPzntNSEnJdF3fxy4vauni7S/e/OL7khOSCf+5X4n/nKCEt62XcpR9zvpvJej/ZPosr+CyTk123jZPrdhfHmZO0lASdT0/Jtd3cumd+uu3F5ysvF8eXkJac67qc2r343F8aS04xXNxecrJkGJnliz+HvNrIr+Tzufff/Mjt2Bfvfz49RUmpyReOl0u4l/49Z2vrfNaG+QrN2d759JQL/eZ8ss4np+S43cXHurjfJaSkZfs8Lt7n0vZcjn3ROudnnSAlpqRJSVn7p+b63gtLVhw5xZ1TrC7bpedcfyXsaenK6tUJqWnKyOrURcil7/Ny5/zi7S/+dzIuLi7bH3lun1lcDv8wXXrc/H7HXbxdtnZz+W6+3L8bOX2vXfpdksNh8ozzcscvyPdMQRSF/8fl+D7OX1Qu4D/pub23nI53JdvCfO7ul1fSJ4qSrBiNPL5DbEZeWxRxR44cUaVKlfTdd9+pZcuWzvpJkyZpw4YN2rJlS7Z9oqOjNWXKFHeGCQAAAADwQH/99ZcqV66c6/piP9Jdrlw5ORwOHTt2zKX+2LFjCgkJyXGfqKgoTZgwwfk6IyNDp06dUtmyZWWz2UyN91oVFxenKlWq6K+//lKpUqWsDgcm4Bx7Ps6x5+McezbOr+fjHHs+znHRYhiGzp07p7CwsMtuV+yTbh8fHzVt2lRff/21evXqJSkzif7666/1wAMP5LiPr6+vfH19XeqCg4NNjhSSVKpUKf6B8HCcY8/HOfZ8nGPPxvn1fJxjz8c5LjqCgoLy3KbYJ92SNGHCBA0ZMkTNmjXTjTfeqJkzZyohIcF5N3MAAAAAAKzgEUl3v379dOLECT399NM6evSoGjdurNWrV2e7uRoAAAAAAO7kEUm3JD3wwAO5TieH9Xx9fTV58uRs0/rhOTjHno9z7Pk4x56N8+v5OMeej3NcPBX7u5cDAAAAAFBU2a0OAAAAAAAAT0XSDQAAAACASUi6AQAAAAAwCUk3LPP555+rRYsW8vf3V+nSpZ3PWYfnSE5OVuPGjWWz2bRjxw6rw0EhOXjwoIYPH66IiAj5+/urRo0amjx5slJSUqwODVfhtddeU3h4uPz8/NSiRQv98MMPVoeEQjJ9+nQ1b95cJUuWVIUKFdSrVy/t3bvX6rBgkueee042m03jxo2zOhQUon/++UeDBg1S2bJl5e/vr+uvv15bt261OizkE0k3LLF8+XLdc889GjZsmHbu3Klvv/1Wd999t9VhoZBNmjRJYWFhVoeBQvbbb78pIyNDb775pnbv3q1XXnlFb7zxhh5//HGrQ0MBffDBB5owYYImT56s7du3q1GjRurUqZOOHz9udWgoBBs2bNCYMWP0/fff66uvvlJqaqo6duyohIQEq0NDIfvxxx/15ptvqmHDhlaHgkJ0+vRptWrVSt7e3lq1apV+/fVXzZgxQ6VLl7Y6NOQTdy+H26WlpSk8PFxTpkzR8OHDrQ4HJlm1apUmTJig5cuXq379+vrpp5/UuHFjq8OCSV588UW9/vrr+vPPP60OBQXQokULNW/eXHPmzJEkZWRkqEqVKnrwwQf12GOPWRwdCtuJEydUoUIFbdiwQbfeeqvV4aCQxMfHq0mTJpo7d66eeeYZNW7cWDNnzrQ6LBSCxx57TN9++602bdpkdSgoIEa64Xbbt2/XP//8I7vdrhtuuEGhoaHq0qWLdu3aZXVoKCTHjh3TiBEj9O677yogIMDqcOAGZ8+eVZkyZawOAwWQkpKibdu2qX379s46u92u9u3ba/PmzRZGBrOcPXtWkvib9TBjxoxRt27dXP6W4Rk+++wzNWvWTH379lWFChV0ww036K233rI6LFwBkm64XdZIWHR0tJ588kmtXLlSpUuXVmRkpE6dOmVxdLhahmFo6NChuv/++9WsWTOrw4Eb7Nu3T7Nnz9aoUaOsDgUFcPLkSaWnp6tixYou9RUrVtTRo0ctigpmycjI0Lhx49SqVSs1aNDA6nBQSJYuXart27dr+vTpVocCE/z55596/fXXVatWLa1Zs0b/+te/9NBDD2nRokVWh4Z8IulGoXnsscdks9kuu2RdCypJTzzxhO688041bdpUCxYskM1m07Jlyyx+F8hNfs/v7Nmzde7cOUVFRVkdMq5Qfs/xxf755x917txZffv21YgRIyyKHEB+jRkzRrt27dLSpUutDgWF5K+//tLYsWO1ePFi+fn5WR0OTJCRkaEmTZpo2rRpuuGGGzRy5EiNGDFCb7zxhtWhIZ+8rA4AnuPhhx/W0KFDL7tN9erVFRsbK0mqV6+es97X11fVq1fX4cOHzQwRVyG/53fdunXavHmzfH19XdY1a9ZMAwcO5FfZIiy/5zjLkSNH1LZtW918882aN2+eydHBLOXKlZPD4dCxY8dc6o8dO6aQkBCLooIZHnjgAa1cuVIbN25U5cqVrQ4HhWTbtm06fvy4mjRp4qxLT0/Xxo0bNWfOHCUnJ8vhcFgYIa5WaGioy/+bJalu3bpavny5RRHhSpF0o9CUL19e5cuXz3O7pk2bytfXV3v37lXr1q0lSampqTp48KCqVatmdpgooPye31dffVXPPPOM8/WRI0fUqVMnffDBB2rRooWZIeIq5fccS5kj3G3btnXOVLHbmThVXPn4+Khp06b6+uuvnY9uzMjI0Ndff60HHnjA2uBQKAzD0IMPPqhPPvlEMTExioiIsDokFKJ27drpl19+cakbNmyYrrvuOj366KMk3B6gVatW2R7z9/vvv/P/5mKEpBtuV6pUKd1///2aPHmyqlSpomrVqunFF1+UJPXt29fi6HC1qlat6vK6RIkSkqQaNWowsuIh/vnnH0VGRqpatWp66aWXdOLECec6RkaLpwkTJmjIkCFq1qyZbrzxRs2cOVMJCQkaNmyY1aGhEIwZM0ZLlizRp59+qpIlSzqv1Q8KCpK/v7/F0eFqlSxZMtv1+YGBgSpbtizX7XuI8ePH6+abb9a0adN011136YcfftC8efOYZVaMkHTDEi+++KK8vLx0zz336Pz582rRooXWrVvH8waBYuCrr77Svn37tG/fvmw/pPAUyuKpX79+OnHihJ5++mkdPXpUjRs31urVq7PdXA3F0+uvvy5JioyMdKlfsGBBnpeUALBe8+bN9cknnygqKkpTp05VRESEZs6cqYEDB1odGvKJ53QDAAAAAGASLsIDAAAAAMAkJN0AAAAAAJiEpBsAAAAAAJOQdAMAAAAAYBKSbgAAAAAATELSDQAAAACASUi6AQAAAAAwCUk3AAAAAAAmIekGABRrkZGRGjdunNVhXLWFCxcqODj4sttER0ercePGV32swmrHSikpKapZs6a+++47q0MpkJtuuknLly+3OgwAgBuQdAMAcImPP/5YHTp0UPny5VWqVCm1bNlSa9asuew+MTExstls2ZYnn3zSTVFfsHz5ckVGRiooKEglSpRQw4YNNXXqVJ06dcrtsZjljTfeUEREhG6++eZs60aNGiWHw6Fly5a5Pa6FCxe6nP8SJUqoadOm+vjjj122e/LJJ/XYY48pIyPD7TECANyLpBsAgEts3LhRHTp00BdffKFt27apbdu26t69u3766ac89927d69iY2Ody2OPPeaGiC944okn1K9fPzVv3lyrVq3Srl27NGPGDO3cuVPvvvuuW2Mxi2EYmjNnjoYPH55tXWJiopYuXapJkybp7bfftiA6qVSpUs7z/9NPP6lTp0666667tHfvXuc2Xbp00blz57Rq1SpLYgQAuA9JNwDAo5w+fVqDBw9W6dKlFRAQoC5duuiPP/5w2eatt95SlSpVFBAQoN69e+vll192mdo9c+ZMTZo0Sc2bN1etWrU0bdo01apVS//3f/+X5/ErVKigkJAQ51KiRIl8x3Wp5557ThUrVlTJkiU1fPhwJSUlXXb7H374QdOmTdOMGTP04osv6uabb1Z4eLg6dOig5cuXa8iQIS7bv/vuuwoPD1dQUJD69++vc+fOOdetXr1arVu3VnBwsMqWLavbb79d+/fvd64/ePCgbDabPv74Y7Vt21YBAQFq1KiRNm/e7HKMvD5rSfr000/VpEkT+fn5qXr16poyZYrS0tJyfZ/btm3T/v371a1bt2zrli1bpnr16umxxx7Txo0b9ddff+XaTkZGhipXrqzXX3/dpf6nn36S3W7XoUOHZBiGoqOjVbVqVfn6+iosLEwPPfRQrm1Kks1mc57/WrVq6ZlnnpHdbtfPP//s3MbhcKhr165aunTpZdsCABR/JN0AAI8ydOhQbd26VZ999pk2b94swzDUtWtXpaamSpK+/fZb3X///Ro7dqx27NihDh066Nlnn71smxkZGTp37pzKlCljWlyX+vDDDxUdHa1p06Zp69atCg0N1dy5cy97jMWLF6tEiRIaPXp0jusvTnb379+vFStWaOXKlVq5cqU2bNig5557zrk+ISFBEyZM0NatW/X111/Lbrerd+/e2aZDP/HEE3rkkUe0Y8cO1a5dWwMGDHAmzPn5rDdt2qTBgwdr7Nix+vXXX/Xmm29q4cKFlz0nmzZtUu3atVWyZMls6+bPn69BgwYpKChIXbp00cKFC3Ntx263a8CAAVqyZIlL/eLFi9WqVStVq1ZNy5cv1yuvvKI333xTf/zxh1asWKHrr78+1zYvlZ6erkWLFkmSmjRp4rLuxhtv1KZNm/LdFgCgmDIAACjG2rRpY4wdO9YwDMP4/fffDUnGt99+61x/8uRJw9/f3/jwww8NwzCMfv36Gd26dXNpY+DAgUZQUFCux3j++eeN0qVLG8eOHct1m/Xr1xuSjMDAQJfl5MmT+YprwYIFLjG0bNnSGD16tMsxWrRoYTRq1CjXGLp06WI0bNgw1/VZJk+ebAQEBBhxcXHOuokTJxotWrTIdZ8TJ04YkoxffvnFMAzDOHDggCHJ+M9//uPcZvfu3YYkY8+ePYZh5O+zbteunTFt2jSXbd59910jNDQ011jGjh1r3Hbbbdnqf//9d8Pb29s4ceKEYRiG8cknnxgRERFGRkZGrm399NNPhs1mMw4dOmQYhmGkp6cblSpVMl5//XXDMAxjxowZRu3atY2UlJRc27jYggULXPqB3W43fH19jQULFmTb9tNPPzXsdruRnp6er7YBAMUTI90AAI+xZ88eeXl5qUWLFs66smXLqk6dOtqzZ4+kzGuub7zxRpf9Ln19sSVLlmjKlCn68MMPVaFChTxj2LRpk3bs2OFcSpcuna+4cnovF28vSS1btrzssQ3DyDO+LOHh4S4jxaGhoTp+/Ljz9R9//KEBAwaoevXqKlWqlMLDwyVJhw8fdmmnYcOGLm1IcraTn896586dmjp1qkqUKOFcRowYodjYWCUmJuYY+/nz5+Xn55et/u2331anTp1Urlw5SVLXrl119uxZrVu3LtfPoXHjxqpbt65ztHvDhg06fvy4+vbtK0nq27evzp8/r+rVq2vEiBH65JNPLjv1XZJKlizpPP8//fSTpk2bpvvvvz/b5Qn+/v7KyMhQcnLyZdsDABRvXlYHAABAUbV06VLdd999WrZsmdq3b5+vfSIiIvJ89JdZateurW+++Uapqany9va+7LaXrrfZbC5Tx7t3765q1arprbfeUlhYmDIyMtSgQQOlpKTk2o7NZpOkK7ojd3x8vKZMmaI77rgj27qcEmtJKleunH755ReXuqxp3EePHpWXl5dL/dtvv6127drlGsPAgQO1ZMkSPfbYY1qyZIk6d+6ssmXLSpKqVKmivXv3au3atfrqq680evRovfjii9qwYUOun7HdblfNmjWdrxs2bKgvv/xSzz//vLp37+6sP3XqlAIDA+Xv759rbACA4o+RbgCAx6hbt67S0tK0ZcsWZ91///tf7d27V/Xq1ZMk1alTRz/++KPLfpe+lqT3339fw4YN0/vvv5/jDbsKO66c9rl4e0n6/vvvL3ucu+++W/Hx8ble+33mzJl8xZsV25NPPql27dqpbt26On36dL72vVh+PusmTZpo7969qlmzZrbFbs/5vyk33HCDfvvtN5eR/S+++ELnzp3TTz/95DLT4P3339fHH3982fd+9913a9euXdq2bZs++ugjDRw40GW9v7+/unfvrldffVUxMTHavHlztqQ/Lw6HQ+fPn3ep27Vrl2644YYragcAUPww0g0A8Bi1atVSz549NWLECL355psqWbKkHnvsMVWqVEk9e/aUJD344IO69dZb9fLLL6t79+5at26dVq1a5RyllTKnlA8ZMkSzZs1SixYtdPToUUmZyVdQUJApcV1q7NixGjp0qJo1a6ZWrVpp8eLF2r17t6pXr57rcVq0aKFJkybp4Ycf1j///KPevXsrLCxM+/bt0xtvvKHWrVtr7NixecZbunRplS1bVvPmzVNoaKgOHz5coEef5eezfvrpp3X77beratWq6tOnj+x2u3bu3Kldu3bpmWeeybHdtm3bKj4+Xrt371aDBg0kZd5ArVu3bmrUqJHLtvXq1dP48eO1ePFijRkzJsf2wsPDdfPNN2v48OFKT09Xjx49nOsWLlyo9PR0tWjRQgEBAXrvvffk7++vatWq5fq+DcNw9pnz58/rq6++0po1a/T000+7bLdp0yZ17NjxMp8gAMATMNINAPAoCxYsUNOmTXX77berZcuWMgxDX3zxhXMqcKtWrfTGG2/o5ZdfVqNGjbR69WqNHz/eZSrzvHnzlJaWpjFjxig0NNS55CdhLWhcl+rXr5+eeuopTZo0SU2bNtWhQ4f0r3/9K8/jPP/881qyZIm2bNmiTp06qX79+powYYIaNmyY7ZFhubHb7Vq6dKm2bdumBg0aaPz48XrxxRev6P1K+fusO3XqpJUrV+rLL79U8+bNddNNN+mVV165bFJbtmxZ9e7dW4sXL5YkHTt2TJ9//rnuvPPOHN9L7969NX/+/MvGOnDgQO3cuVO9e/d2me4dHByst956S61atVLDhg21du1a/d///Z9z+nlO4uLinH2mbt26mjFjhqZOnaonnnjCuc0///yj7777TsOGDbtsXACA4s9mXMldVwAA8EAjRozQb7/9xuOb3KCwPuuff/5ZHTp00P79+53PQi9OHn30UZ0+fVrz5s2zOhQAgMmYXg4AuOa89NJL6tChgwIDA7Vq1SotWrQoz2dgo2DM+qwbNmyo559/XgcOHLii52YXFRUqVNCECROsDgMA4AaMdAMArjl33XWXYmJidO7cOVWvXl0PPvig7r//fqvD8kh81gCAax1JNwAAAAAAJuFGagAAAAAAmISkGwAAAAAAk5B0AwAAAABgEpJuAAAAAABMQtINAAAAAIBJSLoBAAAAADAJSTcAAAAAACYh6QYAAAAAwCQk3QAAAAAAmOT/AU6lzRXAOrUuAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plot_logfc_histograms(directLFQ)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/run.py b/run.py index 7416ec9..75961ad 100644 --- a/run.py +++ b/run.py @@ -471,11 +471,23 @@ def main() -> str: df_fragment_max_peptide, ) = retention_window_searches(mzml_dict, peptide_df, config, perc_95) - log_info("Adding the PSM identifier to fragments...") + log_info("Adding columns ['psm_id', 'scannr', 'stripped_peptide', 'proteins'] from PSMs to fragment df...") df_fragment = df_fragment.join( - df_psms.select(["psm_id", "scannr"]), on="psm_id", how="left" + df_psms.select(["psm_id", "scannr", "stripped_peptide", "proteins"]), on="psm_id", how="left" ) + log_info("Adding fragment names to fragments...") + df_fragment = df_fragment.with_columns( + pl.Series( + "fragment_name", + df_fragment["fragment_type"] + + df_fragment["fragment_ordinals"] + + "/" + + df_fragment["fragment_charge"], + ) + ) + + # Narrow types for static analysis assert isinstance(df_fragment, pl.DataFrame) assert isinstance(df_psms, pl.DataFrame) @@ -517,6 +529,14 @@ def main() -> str: config["sage_basic"]["mzml_paths"][0] # TODO: should be for all mzml files ) + # Add retention time margins to precursor ions by going right and left from + # apex retention time in MS2 scans until fragment intensity drops below 1% of apex + log_info("Adding retention time margins to precursor ions and fragments...") + df_psms, df_fragment = mumdia.add_retention_time_margins_loop( + df_psms=df_psms, df_fragment=df_fragment, top_n=100, intensity_threshold=0.01 + ) + + import pickle with open("debug/ms1_dict.pkl", "wb") as f: diff --git a/utilities/plotting.py b/utilities/plotting.py index deae7e2..af734c5 100644 --- a/utilities/plotting.py +++ b/utilities/plotting.py @@ -107,3 +107,156 @@ def plot_XIC(df: pl.DataFrame, output_dir: str = "xics"): plt.tight_layout() plt.savefig(f"{output_dir}/{precursor}_XIC.svg") + + +def plot_XIC_with_margins(df: pl.DataFrame, output_dir: str = "debug/calibration_xics", adapted_interval=None, min_interval=None, max_interval=None, apex_rt=None, cutoff=None): + """ + Plots fragment_intensity vs rt for each unique fragment_name. + Colors by fragment_name, lines connect fragments, marker shape by psm_id. + Adds two separate legends: one for fragment_name (colors), one for psm_id (shapes). + Works with a Polars DataFrame. + In addition to the normal plot_XIC, this function adds vertical lines for the provided intervals: + - adapted_interval: tuple (left, right) for the adapted margins + - min_interval: tuple (left, right) for the minimum RT interval + - max_interval: tuple (left, right) for the maximum RT interval + - apex_rt: float for the apex retention time + """ + # Convert to pandas first + pdf = df.to_pandas() + + precursor = (pdf["peptide"] + "_" + pdf["charge"].astype(str)).unique()[0] + + # Validate required columns + required_cols = {"fragment_intensity", "rt", "fragment_name", "psm_id"} + missing = required_cols - set(pdf.columns) + if missing: + raise ValueError(f"DataFrame missing required columns: {missing}") + + # Unique colors for fragment_name + fragment_names = pdf["fragment_name"].unique() + colors = plt.cm.get_cmap("tab20", len(fragment_names)) + fragment_color_map = {frag: colors(i) for i, frag in enumerate(fragment_names)} + + # Unique marker styles for psm_id (repeat if there are many) + marker_styles = ["o", "s", "D", "^", "v", "p", "*", "X", "H", "<", ">"] + psm_ids = pdf["psm_id"].unique() + psm_marker_map = { + psm: marker_styles[i % len(marker_styles)] for i, psm in enumerate(psm_ids) + } + + plt.figure(figsize=(10, 6)) + + # Step 1: Plot continuous lines by fragment_name + for frag in fragment_names: + frag_df = pdf[pdf["fragment_name"] == frag].sort_values("rt") + plt.plot( + frag_df["rt"], + frag_df["fragment_intensity"], + color=fragment_color_map[frag], + linestyle="-", + linewidth=1, + ) + + # Step 2: Overlay markers by fragment_name + psm_id + for frag in fragment_names: + frag_df = pdf[pdf["fragment_name"] == frag] + for psm in psm_ids: + psm_df = frag_df[frag_df["psm_id"] == psm] + if psm_df.empty: + continue + plt.scatter( + psm_df["rt"], + psm_df["fragment_intensity"], + color=fragment_color_map[frag], + marker=psm_marker_map[psm], + edgecolors="black", + linewidths=0.5, + ) + + plt.xlabel("Retention Time (RT)") + plt.ylabel("Fragment Intensity") + plt.title("Extracted Ion Chromatogram by Fragment") + + # --- Create two legends manually --- + # Legend 1: fragment_name (color lines) + frag_legend_elements = [ + Line2D([0], [0], color=fragment_color_map[frag], lw=2, label=frag) + for frag in fragment_names + ] + legend1 = plt.legend( + handles=frag_legend_elements, + title="Fragment", + bbox_to_anchor=(1.05, 1), + loc="upper left", + ) + plt.gca().add_artist(legend1) # add first legend manually + + # Legend 2: psm_id (marker shapes) + psm_legend_elements = [ + Line2D( + [0], + [0], + marker=psm_marker_map[psm], + color="w", + markerfacecolor="gray", + markeredgecolor="black", + markersize=8, + label=str(psm), + ) + for psm in psm_ids + ] + + + # add vertical lines for intervals + if adapted_interval: + plt.axvline(x=adapted_interval[0], color='gray', linestyle='--', label='Adapted Interval Start') + plt.axvline(x=adapted_interval[1], color='gray', linestyle='--', label='Adapted Interval End') + + if min_interval: + plt.axvline(x=min_interval[0], color='blue', linestyle='--', label='Min Interval Start') + plt.axvline(x=min_interval[1], color='blue', linestyle='--', label='Min Interval End') + + if max_interval: + plt.axvline(x=max_interval[0], color='red', linestyle='--', label='Max Interval Start') + plt.axvline(x=max_interval[1], color='red', linestyle='--', label='Max Interval End') + + if apex_rt: + plt.axvline(x=apex_rt, color='green', linestyle='-', label='Apex RT') + + if cutoff: + plt.axhline(y=cutoff, color='purple', linestyle='-.', label='Cutoff Intensity') + + plt.legend( + handles=psm_legend_elements, + title="PSM ID", + bbox_to_anchor=(1.05, 0), + loc="lower left", + ) + + plt.title(f"{precursor}") + + plt.tight_layout() + plt.savefig(f"{output_dir}/{precursor}_XIC.svg") + +def plot_rt_margin_histogram(rt_margins, output_dir: str = "debug/calibration_xics", min_diff=None, max_diff=None): + """ + Plots a histogram of the rt_margins from the PSM DataFrame. + Expects rt_margins to be a list of tuples (left_margin, right_margin). + """ + if not rt_margins: + raise ValueError("rt_margins list is empty.") + + plt.figure(figsize=(10, 6)) + plt.hist(rt_margins, bins=100, alpha=0.5, label='Margins', color='orange') + + if min_diff is not None: + plt.axvline(x=min_diff, color='red', linestyle='--', label='Min RT Margin') + if max_diff is not None: + plt.axvline(x=max_diff, color='green', linestyle='--', label='Max RT Margin') + + plt.xlabel('Retention Time Margin') + plt.ylabel('Frequency') + plt.title('Histogram of Retention Time Margins') + plt.legend() + plt.tight_layout() + plt.savefig(f"{output_dir}/rt_margin_histogram.svg") \ No newline at end of file