Merge pull request #524 from bigbio/dev

Dev

Author: ypriverol

pmultiqc/modules/common/common_utils.py

@@ -329,16 +329,16 @@ def hist_compute(rt_list, rt_range):
 
 def evidence_calibrated_mass_error(
     evidence_data,
-    recommpute=False,
+    recompute=False,
     filter_outliers_ppm: bool = False
 ):
     # filter_outliers_ppm (if True): Remove rows with mass error [ppm] greater than 1000 (Default: False)
 
     if "potential contaminant" in evidence_data.columns:
         evidence_data = evidence_data[evidence_data["potential contaminant"] != "+"].copy()
 
-    if recommpute:
-        evd_df = recommpute_mass_error(evidence_data)
+    if recompute:
+        evd_df = recompute_mass_error(evidence_data)
     else:
         evd_df = evidence_data.copy()
 
@@ -395,7 +395,7 @@ def evidence_calibrated_mass_error(
     return result_dict
 
 # re-compute mass error
-def recommpute_mass_error(evidence_df):
+def recompute_mass_error(evidence_df):
     required_cols = [
         "mass error [ppm]",
         "uncalibrated mass error [ppm]",
@@ -407,7 +407,7 @@ def recommpute_mass_error(evidence_df):
     ]
 
     if not all(col in evidence_df.columns for col in required_cols):
-        log.info("Evidence is missing one or more required columns in recommpute_mass_error.")
+        log.info("Evidence is missing one or more required columns in recompute_mass_error.")
         return None
 
     df = evidence_df[required_cols].copy()

pmultiqc/modules/common/dia_utils.py

@@ -705,15 +705,62 @@ def cal_rt_irt_loess(report_df, frac=0.3, data_bins: int = DEFAULT_BINS):
     return plot_dict
 
 
-# DIA-NN: Peptides Quantification Table
-def create_peptides_table(report_df, sample_df, file_df):
-    # Validation: remove rows with 0 or NA Precursor.Normalised values
+def _prepare_quant_table_data(report_df):
+    """
+    Common preprocessing for quantification table creation.
+
+    Returns:
+        pd.DataFrame: Preprocessed report data with positive Precursor.Normalised values.
+    """
     report_data = report_df[report_df["Precursor.Normalised"] > 0].copy()
-    report_data = drop_empty_row(report_data, ["Protein.Names", "Stripped.Sequence"])
+    return drop_empty_row(report_data, ["Protein.Names", "Stripped.Sequence"])
+
+
+def _merge_condition_data(report_data, sample_df, file_df):
+    """
+    Merge report data with condition information from sample/file DataFrames.
+
+    Returns:
+        tuple: (merged DataFrame with condition info, list of unique conditions) or (None, [])
+    """
+    if sample_df.empty or file_df.empty:
+        return None, []
+
+    sample_cond_df = pd.merge(
+        sample_df[["Sample", "MSstats_Condition"]],
+        file_df[["Sample", "Spectra_Filepath"]],
+        on="Sample",
+    )
+    # Vectorized path splitting (more efficient than apply with lambda)
+    sample_cond_df["Run"] = sample_cond_df["Spectra_Filepath"].str.rsplit(".", n=1).str[0]
 
+    cond_report_data = pd.merge(
+        report_data[["Stripped.Sequence", "Protein.Names", "Precursor.Normalised", "Run"]],
+        sample_cond_df[["Run", "MSstats_Condition"]].drop_duplicates(),
+        on="Run",
+    )
+
+    unique_conditions = sample_df["MSstats_Condition"].drop_duplicates().tolist()
+    return cond_report_data, unique_conditions
+
+
+def _add_condition_headers(headers, conditions):
+    """Add condition-based headers to the headers dictionary."""
+    for exp_condition in conditions:
+        headers[str(exp_condition)] = {
+            "title": str(exp_condition),
+            "description": "MSstats Condition",
+            "format": "{:,.4f}",
+        }
+
+
+# DIA-NN: Peptides Quantification Table
+def create_peptides_table(report_df, sample_df, file_df):
+    """Create peptides quantification table from DIA-NN report."""
+    report_data = _prepare_quant_table_data(report_df)
     report_data["BestSearchScore"] = 1 - report_data["Q.Value"]
 
-    table_dict = dict()
+    table_dict = {}
     for sequence_protein, group in report_data.groupby(["Stripped.Sequence", "Protein.Names"]):
         table_dict[sequence_protein] = {
             "ProteinName": sequence_protein[1],
@@ -737,52 +784,29 @@ def create_peptides_table(report_df, sample_df, file_df):
         },
     }
 
-    if not sample_df.empty and not file_df.empty:
-
-        sample_cond_df = pd.merge(
-            sample_df[["Sample", "MSstats_Condition"]],
-            file_df[["Sample", "Spectra_Filepath"]],
-            on="Sample",
-        )
-        sample_cond_df["Run"] = sample_cond_df["Spectra_Filepath"].apply(
-            lambda x: os.path.splitext(x)[0]
-        )
-
-        cond_report_data = pd.merge(
-            report_data[["Stripped.Sequence", "Protein.Names", "Precursor.Normalised", "Run"]],
-            sample_cond_df[["Run", "MSstats_Condition"]].drop_duplicates(),
-            on="Run",
-        )
-
+    cond_report_data, unique_conditions = _merge_condition_data(report_data, sample_df, file_df)
+    if cond_report_data is not None:
         for sequence_protein, group in cond_report_data.groupby(
                 ["Stripped.Sequence", "Protein.Names"]
         ):
-
-            condition_data = dict()
-            for condition, sub_group in group.groupby("MSstats_Condition"):
-                condition_data[str(condition)] = np.log10(sub_group["Precursor.Normalised"].mean())
-
+            condition_data = {
+                str(cond): np.log10(sub_group["Precursor.Normalised"].mean())
+                for cond, sub_group in group.groupby("MSstats_Condition")
+            }
             table_dict[sequence_protein].update(condition_data)
 
-        for exp_condition in sample_df["MSstats_Condition"].drop_duplicates():
-            headers[str(exp_condition)] = {
-                "title": str(exp_condition),
-                "description": "MSstats Condition",
-                "format": "{:,.4f}",
-            }
+        _add_condition_headers(headers, unique_conditions)
 
     result_dict = {i: v for i, (_, v) in enumerate(table_dict.items(), start=1)}
-
     return result_dict, headers
 
 
 # DIA-NN: Protein Quantification Table
 def create_protein_table(report_df, sample_df, file_df):
-    # Validation: remove rows with 0 or NA Precursor.Normalised values
-    report_data = report_df[report_df["Precursor.Normalised"] > 0].copy()
-    report_data = drop_empty_row(report_data, ["Protein.Names", "Stripped.Sequence"])
+    """Create protein quantification table from DIA-NN report."""
+    report_data = _prepare_quant_table_data(report_df)
 
-    table_dict = dict()
+    table_dict = {}
     for protein_name, group in report_data.groupby("Protein.Names"):
         table_dict[protein_name] = {
             "ProteinName": protein_name,
@@ -807,40 +831,18 @@ def create_protein_table(report_df, sample_df, file_df):
         },
     }
 
-    if not sample_df.empty and not file_df.empty:
-
-        sample_cond_df = pd.merge(
-            sample_df[["Sample", "MSstats_Condition"]],
-            file_df[["Sample", "Spectra_Filepath"]],
-            on="Sample",
-        )
-        sample_cond_df["Run"] = sample_cond_df["Spectra_Filepath"].apply(
-            lambda x: os.path.splitext(x)[0]
-        )
-
-        cond_report_data = pd.merge(
-            report_data[["Stripped.Sequence", "Protein.Names", "Precursor.Normalised", "Run"]],
-            sample_cond_df[["Run", "MSstats_Condition"]].drop_duplicates(),
-            on="Run",
-        )
-
+    cond_report_data, unique_conditions = _merge_condition_data(report_data, sample_df, file_df)
+    if cond_report_data is not None:
         for protein_name, group in cond_report_data.groupby("Protein.Names"):
-
-            condition_data = dict()
-            for condition, sub_group in group.groupby("MSstats_Condition"):
-                condition_data[str(condition)] = np.log10(sub_group["Precursor.Normalised"].mean())
-
+            condition_data = {
+                str(cond): np.log10(sub_group["Precursor.Normalised"].mean())
+                for cond, sub_group in group.groupby("MSstats_Condition")
+            }
             table_dict[protein_name].update(condition_data)
 
-        for exp_condition in sample_df["MSstats_Condition"].drop_duplicates():
-            headers[str(exp_condition)] = {
-                "title": str(exp_condition),
-                "description": "MSstats Condition",
-                "format": "{:,.4f}",
-            }
+        _add_condition_headers(headers, unique_conditions)
 
     result_dict = {i: v for i, (_, v) in enumerate(table_dict.items(), start=1)}
-
     return result_dict, headers
 
 

pmultiqc/modules/common/ms/msinfo.py

@@ -139,11 +139,11 @@ def parse(self, **_kwargs) -> None:
                 self.enable_dia,
             )
 
-            for m in mzml_table.keys():
-                if mzml_table[m]["MS2_Num"] > 0:
-                    heatmap_charge[m] = mzml_table[m]["Charge_2"] / mzml_table[m]["MS2_Num"]
-                else:
-                    heatmap_charge[m] = 0
+            # Calculate heatmap_charge only for current file (not all files in each iteration)
+            if mzml_table[m_name]["MS2_Num"] > 0:
+                heatmap_charge[m_name] = mzml_table[m_name]["Charge_2"] / mzml_table[m_name]["MS2_Num"]
+            else:
+                heatmap_charge[m_name] = 0
 
             self.log.info(
                 "{}: Done aggregating ms_statistics dataframe {}...".format(

pmultiqc/modules/common/plots/dia.py

@@ -574,7 +574,21 @@ def draw_loess_rt_irt(sub_section, plot_data):
     )
 
 def calculate_dia_intensity_std(df, sdrf_file_df):
-
+    """
+    Calculate standard deviation of intensity for DIA data.
+
+    Parameters:
+    -----------
+    df : pd.DataFrame
+        DataFrame with Run, Modified.Sequence, Protein.Group, and log_intensity columns.
+    sdrf_file_df : pd.DataFrame
+        SDRF file DataFrame (can be empty).
+
+    Returns:
+    --------
+    dict or None: Dictionary mapping sample/condition to list of log intensity std values,
+                  or None if calculation cannot be performed.
+    """
     df_sub = df[["Run", "Modified.Sequence", "Protein.Group", "log_intensity"]].copy()
 
     if not sdrf_file_df.empty:
@@ -617,8 +631,9 @@ def calculate_dia_intensity_std(df, sdrf_file_df):
         }
 
         return plot_data
-    else:
-        log.warning("No SDRF available; failed to parse experimental groups; SD Intensity not generated.")
+
+    log.warning("No SDRF available; failed to parse experimental groups; SD Intensity not generated.")
+    return None
 
 def extract_condition_and_replicate(run_name):
 

pmultiqc/modules/common/plots/general.py

@@ -197,7 +197,9 @@ def draw_exp_design(sub_sections, exp_design):
             "description": "",
             "scale": False,
         }}
-        for k, _ in condition_split(sample_df_slice["MSstats_Condition"].iloc[0]).items():
+        # Use first row of sample_df for condition keys (safer than relying on loop variable)
+        first_condition = sample_df["MSstats_Condition"].iloc[0] if not sample_df.empty else ""
+        for k, _ in condition_split(first_condition).items():
             headers["MSstats_Condition_" + str(k)] = {
                 "title": "MSstats Condition: " + str(k),
                 "description": "",

pmultiqc/modules/common/stats.py

@@ -25,13 +25,25 @@ def nanmedian(values: np.ndarray, all_nan_fallback: np.float64) -> np.float64:
 
 def qual_uniform(group_df_rt):
     """
+    Calculate quality score based on uniformity of retention time distribution.
+
     Parameters:
     -----------
     group_df_rt: group["Retention time"] or group["retention_time"]
 
+    Returns:
+    --------
+    float: Quality score between 0 and 1, where 1 indicates perfect uniformity.
     """
-    x = group_df_rt / np.nansum(group_df_rt)
     n = group_df_rt.notna().sum()
+    if n == 0:
+        return 0.0
+
+    total_sum = np.nansum(group_df_rt)
+    if total_sum == 0:
+        return 0.0
+
+    x = group_df_rt / total_sum
     y = np.nansum(x) / n
     worst = ((1 - y) ** 0.5) * 1 / n + (y**0.5) * (n - 1) / n
     sc = np.sum(np.abs(x - y) ** 0.5) / n
@@ -40,25 +52,47 @@ def qual_uniform(group_df_rt):
     return result
 
 
-def cal_delta_mass_dict(df, col):
+def cal_delta_mass_dict(df, col, num_bins: int = 1000):
+    """
+    Calculate delta mass distribution as counts and frequencies.
+
+    Parameters:
+    -----------
+    df : pd.DataFrame
+        DataFrame containing the mass delta column.
+    col : str
+        Name of the column containing mass delta values.
+    num_bins : int, optional
+        Number of bins for histogram (default: 1000).
+
+    Returns:
+    --------
+    dict: Dictionary with 'count' and 'frequency' keys containing binned data.
+    """
+    # Compute value_counts once and derive frequency from counts
+    count_bin = df[col].value_counts(sort=False, bins=num_bins)
 
-    count_bin = df[col].value_counts(sort=False, bins=1000)
-    count_bin_data = dict()
-    for index in count_bin.index:
-        count_bin_data[float(index.mid)] = int(count_bin[index])
+    # Build count dictionary
+    count_bin_data = {
+        float(interval.mid): int(count)
+        for interval, count in count_bin.items()
+    }
 
-    frequency_bin = df[col].value_counts(sort=False, bins=1000, normalize=True)
-    frequency_bin_data = dict()
-    for index in frequency_bin.index:
-        frequency_bin_data[float(index.mid)] = float(frequency_bin[index])
+    # Derive frequency from counts (more efficient than calling value_counts twice)
+    total_count = count_bin.sum()
+    if total_count > 0:
+        frequency_bin_data = {
+            float(interval.mid): float(count / total_count)
+            for interval, count in count_bin.items()
+        }
+    else:
+        frequency_bin_data = {k: 0.0 for k in count_bin_data.keys()}
 
-    delta_mass = {
+    return {
         "count": count_bin_data,
         "frequency": frequency_bin_data,
     }
 
-    return delta_mass
-
 
 def cal_hm_charge(df: pd.DataFrame, run_col: str, charge_col: str):
 

pmultiqc/modules/maxquant/maxquant_utils.py

@@ -19,7 +19,7 @@
 from pmultiqc.modules.common.common_utils import (
     mods_statistics,
     evidence_rt_count,
-    recommpute_mass_error,
+    recompute_mass_error,
     evidence_calibrated_mass_error
 )
 
@@ -910,7 +910,7 @@ def evidence_uncalibrated_mass_error(evidence_data):
     if "potential contaminant" in evidence_data.columns:
         evidence_data = evidence_data[evidence_data["potential contaminant"] != "+"].copy()
 
-    evd_df = recommpute_mass_error(evidence_data)
+    evd_df = recompute_mass_error(evidence_data)
 
     if evd_df is None:
         if any(