diff --git a/coverage.xml b/coverage.xml new file mode 100644 index 0000000..b080aa7 --- /dev/null +++ b/coverage.xml @@ -0,0 +1,146 @@ + + + + + + /mnt/data/figure2a_fair_project/src + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/pixi.toml b/pixi.toml new file mode 100644 index 0000000..6fd1781 --- /dev/null +++ b/pixi.toml @@ -0,0 +1,38 @@ +[project] +name = "figure2a-reproduction" +version = "0.1.0" +description = "Reproduce Figure 2A from Kelliher et al. 2016 with a FAIR Python workflow." +channels = ["conda-forge"] +platforms = ["linux-64", "osx-64", "osx-arm64", "win-64"] + +[dependencies] +python = ">=3.10,<3.13" +pandas = ">=2.2" +openpyxl = ">=3.1" +numpy = ">=1.26" +matplotlib = ">=3.8" + +[feature.dev.dependencies] +ruff = ">=0.6" +pyright = ">=1.1" +pytest = ">=8.0" +pytest-cov = ">=5.0" +safety = ">=3.2" +trivy = ">=0.56" + +[environments] +default = { features = ["default"] } +dev = { features = ["default", "dev"] } + +[tasks] +run-figure = "python src/reproduce_figure_2a.py --input data/pgen.1006453.s002.xlsx --output results/figure_2A_reproduced.png" + +[feature.dev.tasks] +lint = "ruff check ." +format = "ruff format ." +typecheck = "pyright src/reproduce_figure_2a.py tests" +test = "pytest" +coverage = "pytest --cov=src --cov-report=term-missing" +safety-check = "safety check" +trivy-check = "trivy fs --scanners secret,vuln,misconfig ." +all-checks = { depends-on = ["lint", "typecheck", "test"] } diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..7dc58e7 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,29 @@ +[project] +name = "figure2a-reproduction" +version = "0.1.0" +description = "Reproducible Python project to recreate Figure 2A from Kelliher et al. 2016." +authors = [{ name = "Chaker Aloui", email = "aloui-choko@hotmail.fr" }] +readme = "README.md" +license = { text = "MIT" } +requires-python = ">=3.10" +keywords = ["bioinformatics", "transcriptomics", "reproducibility", "figure"] + +[tool.ruff] +line-length = 88 +target-version = "py310" + +[tool.ruff.lint] +select = ["E", "F", "W", "I", "D", "N", "UP", "B", "SIM"] +ignore = ["D203", "D212"] + +[tool.ruff.lint.pydocstyle] +convention = "numpy" + +[tool.pyright] +pythonVersion = "3.10" +typeCheckingMode = "basic" +include = ["src", "tests"] + +[tool.pytest.ini_options] +testpaths = ["tests"] +addopts = "--cov=src --cov-report=term-missing --cov-report=xml" diff --git a/results/figure_2A_reproduced.png b/results/figure_2A_reproduced.png new file mode 100644 index 0000000..b34d2cb Binary files /dev/null and b/results/figure_2A_reproduced.png differ diff --git a/tests/__pycache__/test_reproduce_figure_2a.cpython-313-pytest-9.0.2.pyc b/tests/__pycache__/test_reproduce_figure_2a.cpython-313-pytest-9.0.2.pyc new file mode 100644 index 0000000..3d75ef2 Binary files /dev/null and b/tests/__pycache__/test_reproduce_figure_2a.cpython-313-pytest-9.0.2.pyc differ diff --git a/tests/test_reproduce_figure_2a.py b/tests/test_reproduce_figure_2a.py new file mode 100644 index 0000000..9d9f1d3 --- /dev/null +++ b/tests/test_reproduce_figure_2a.py @@ -0,0 +1,158 @@ +"""Unit tests for the Figure 2A reproduction module.""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) + +from reproduce_figure_2a import ( # noqa: E402 + build_expression_matrix, + compute_row_z_scores, + get_time_columns, + plot_figure_2a, + select_and_order_periodic_genes, +) + + +def make_test_dataframe() -> pd.DataFrame: + """Create a small expression table with Figure 2A-like columns.""" + return pd.DataFrame( + { + "gene": ["gene_a", "gene_b", "gene_c"], + "Figure2A_order_peaktime": [2, np.nan, 1], + 0: [1.0, 2.0, 4.0], + 5: [2.0, 3.0, 4.0], + 10: [3.0, 4.0, 4.0], + }, + ) + + +def test_get_time_columns_returns_sorted_numeric_columns() -> None: + """Check that numeric time columns are detected and sorted.""" + dataframe = make_test_dataframe() + assert get_time_columns(dataframe) == [0, 5, 10] + + +def test_get_time_columns_raises_without_numeric_columns() -> None: + """Check that missing time columns raise a clear error.""" + dataframe = pd.DataFrame({"gene": ["gene_a"]}) + with pytest.raises(ValueError, match="No numeric time columns"): + get_time_columns(dataframe) + + +def test_select_and_order_periodic_genes_uses_descending_order() -> None: + """Check the descending sort needed to match the published orientation.""" + dataframe = make_test_dataframe() + ordered = select_and_order_periodic_genes(dataframe) + assert ordered["gene"].to_list() == ["gene_a", "gene_c"] + + +def test_select_and_order_periodic_genes_raises_for_missing_column() -> None: + """Check that a missing order column is rejected.""" + dataframe = pd.DataFrame({"gene": ["gene_a"], 0: [1.0]}) + with pytest.raises(ValueError, match="Required order column"): + select_and_order_periodic_genes(dataframe) + + +def test_compute_row_z_scores_nominal_case() -> None: + """Check row-wise z-score normalization on non-constant rows.""" + result = compute_row_z_scores([[1.0, 2.0, 3.0]]) + expected = np.array([[-1.22474487, 0.0, 1.22474487]]) + np.testing.assert_allclose(result, expected) + + +def test_compute_row_z_scores_constant_row_becomes_zero() -> None: + """Check that constant rows do not produce NaN values.""" + result = compute_row_z_scores([[4.0, 4.0, 4.0]]) + np.testing.assert_allclose(result, np.zeros((1, 3))) + + +def test_compute_row_z_scores_rejects_invalid_bounds() -> None: + """Check that invalid plotting bounds raise an error.""" + with pytest.raises(ValueError, match="vmin must be lower"): + compute_row_z_scores([[1.0, 2.0]], vmin=1.0, vmax=1.0) + + +def test_build_expression_matrix_shape() -> None: + """Check that expression extraction returns the expected shape.""" + dataframe = make_test_dataframe() + ordered = select_and_order_periodic_genes(dataframe) + matrix = build_expression_matrix(ordered, [0, 5, 10]) + assert matrix.shape == (2, 3) + + +def test_plot_figure_2a_writes_png(tmp_path: Path) -> None: + """Check that the plotting function creates an output file.""" + output_path = tmp_path / "figure.png" + z_scores = np.array([[0.0, 1.0, -1.0], [1.0, 0.0, -1.0]]) + figure = plot_figure_2a(z_scores, [0, 50, 100], output_path, dpi=80) + assert output_path.exists() + assert figure.axes[0].get_ylabel() == "Top Periodic Genes (2)" + + +def test_plot_figure_2a_rejects_dimension_mismatch(tmp_path: Path) -> None: + """Check that inconsistent matrix and time dimensions are rejected.""" + with pytest.raises(ValueError, match="Number of matrix columns"): + plot_figure_2a(np.array([[1.0, 2.0]]), [0], tmp_path / "figure.png") + +from reproduce_figure_2a import ( # noqa: E402 + build_argument_parser, + main, + read_expression_table, + reproduce_figure_2a, +) + + +def test_read_expression_table_rejects_negative_header(tmp_path: Path) -> None: + """Check that a negative header row is rejected before file reading.""" + with pytest.raises(ValueError, match="header_row"): + read_expression_table(tmp_path / "missing.xlsx", header_row=-1) + + +def test_read_expression_table_rejects_missing_file(tmp_path: Path) -> None: + """Check that a missing input file raises FileNotFoundError.""" + with pytest.raises(FileNotFoundError, match="Input Excel file not found"): + read_expression_table(tmp_path / "missing.xlsx") + + +def test_reproduce_figure_2a_end_to_end_with_small_excel(tmp_path: Path) -> None: + """Check the complete workflow on a minimal Excel fixture.""" + input_path = tmp_path / "input.xlsx" + output_path = tmp_path / "figure.png" + dataframe = make_test_dataframe() + with pd.ExcelWriter(input_path) as writer: + dataframe.to_excel(writer, sheet_name="Sheet1", index=False, startrow=2) + + result_path = reproduce_figure_2a(input_path, output_path) + + assert result_path == output_path + assert output_path.exists() + + +def test_build_argument_parser_defaults() -> None: + """Check that the CLI parser has usable default arguments.""" + parser = build_argument_parser() + args = parser.parse_args([]) + assert args.input.name == "pgen.1006453.s002.xlsx" + assert args.output.name == "figure_2A_reproduced.png" + assert not args.no_panel_label + + +def test_main_returns_zero_with_small_excel(tmp_path: Path) -> None: + """Check that the command-line entry point succeeds on valid inputs.""" + input_path = tmp_path / "input.xlsx" + output_path = tmp_path / "figure.png" + dataframe = make_test_dataframe() + with pd.ExcelWriter(input_path) as writer: + dataframe.to_excel(writer, sheet_name="Sheet1", index=False, startrow=2) + + status = main(["--input", str(input_path), "--output", str(output_path)]) + + assert status == 0 + assert output_path.exists() diff --git a/tp_python/reproduce_figure_2a.py b/tp_python/reproduce_figure_2a.py new file mode 100644 index 0000000..2bbd9c9 --- /dev/null +++ b/tp_python/reproduce_figure_2a.py @@ -0,0 +1,449 @@ +"""Reproduce Figure 2A from Kelliher et al. 2016. + +This module reads the supplementary Excel file ``pgen.1006453.s002.xlsx`` and +recreates the Saccharomyces cerevisiae periodic-gene heatmap shown in Figure 2A. +""" + +from __future__ import annotations + +import argparse +import warnings +from collections.abc import Sequence +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import numpy.typing as npt +import pandas as pd +from matplotlib.colors import LinearSegmentedColormap +from matplotlib.figure import Figure +from pandas.api.types import is_numeric_dtype + +DEFAULT_INPUT_PATH = Path("data/pgen.1006453.s002.xlsx") +DEFAULT_OUTPUT_PATH = Path("results/figure_2A_reproduced.png") +DEFAULT_SHEET_NAME = "Sheet1" +DEFAULT_HEADER_ROW = 2 +DEFAULT_ORDER_COLUMN = "Figure2A_order_peaktime" +DEFAULT_VMIN = -1.5 +DEFAULT_VMAX = 1.5 +DEFAULT_DPI = 300 +DEFAULT_XTICK_TIMES = (0, 50, 100, 150, 200) + + +def read_expression_table( + input_path: Path, + sheet_name: str = DEFAULT_SHEET_NAME, + header_row: int = DEFAULT_HEADER_ROW, +) -> pd.DataFrame: + """Read the supplementary Excel expression table. + + Parameters + ---------- + input_path: + Path to the Excel file. + sheet_name: + Name of the sheet to read. + header_row: + Zero-based row index containing column names. + + Returns + ------- + pandas.DataFrame + Expression and annotation table. + + Raises + ------ + FileNotFoundError + If the Excel file does not exist. + ValueError + If ``header_row`` is negative. + """ + if header_row < 0: + msg = "header_row must be a non-negative integer." + raise ValueError(msg) + if not input_path.exists(): + msg = f"Input Excel file not found: {input_path.resolve()}" + raise FileNotFoundError(msg) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return pd.read_excel(input_path, sheet_name=sheet_name, header=header_row) + + +def get_time_columns(dataframe: pd.DataFrame) -> list[int | float]: + """Identify numeric time-point columns and sort them by time. + + Parameters + ---------- + dataframe: + Input expression table. + + Returns + ------- + list[int | float] + Sorted time-point columns. + + Raises + ------ + ValueError + If no numeric time-point columns are found. + """ + time_columns: list[int | float] = [] + for column_name in dataframe.columns: + if isinstance(column_name, int | float | np.integer | np.floating): + time_columns.append(column_name) + elif is_numeric_dtype(dataframe[column_name]) and str(column_name).isdigit(): + time_columns.append(float(column_name)) + + time_columns = sorted(time_columns, key=float) + if not time_columns: + msg = "No numeric time columns were found in the expression table." + raise ValueError(msg) + return time_columns + + +def select_and_order_periodic_genes( + dataframe: pd.DataFrame, + order_column: str = DEFAULT_ORDER_COLUMN, +) -> pd.DataFrame: + """Select genes used in Figure 2A and order them as in the paper. + + The published visual orientation is reproduced by sorting + ``Figure2A_order_peaktime`` in decreasing order. + + Parameters + ---------- + dataframe: + Input expression table. + order_column: + Column containing the Figure 2A order index. + + Returns + ------- + pandas.DataFrame + Filtered and ordered table. + + Raises + ------ + ValueError + If ``order_column`` is missing or contains no valid value. + """ + if order_column not in dataframe.columns: + msg = f"Required order column is missing: {order_column}" + raise ValueError(msg) + + periodic_genes = dataframe.loc[dataframe[order_column].notna()].copy() + periodic_genes[order_column] = pd.to_numeric( + periodic_genes[order_column], + errors="coerce", + ) + periodic_genes = periodic_genes.loc[periodic_genes[order_column].notna()].copy() + if periodic_genes.empty: + msg = f"Column {order_column} does not contain valid Figure 2A order values." + raise ValueError(msg) + + return periodic_genes.sort_values(order_column, ascending=False) + + +def compute_row_z_scores( + expression_values: npt.ArrayLike, + vmin: float = DEFAULT_VMIN, + vmax: float = DEFAULT_VMAX, +) -> npt.NDArray[np.float64]: + """Compute row-wise z-scores and clip values to the plotting range. + + Parameters + ---------- + expression_values: + Two-dimensional matrix with genes in rows and time points in columns. + vmin: + Lower clipping bound. + vmax: + Upper clipping bound. + + Returns + ------- + numpy.ndarray + Row-wise normalized expression matrix. + + Raises + ------ + ValueError + If bounds are invalid or the expression matrix is not two-dimensional. + """ + if vmin >= vmax: + msg = "vmin must be lower than vmax." + raise ValueError(msg) + + matrix = np.asarray(expression_values, dtype=float) + if matrix.ndim != 2: + msg = "expression_values must be a two-dimensional matrix." + raise ValueError(msg) + + row_means = np.nanmean(matrix, axis=1, keepdims=True) + row_standard_deviations = np.nanstd(matrix, axis=1, ddof=0, keepdims=True) + row_standard_deviations[row_standard_deviations == 0] = np.nan + + z_scores = (matrix - row_means) / row_standard_deviations + z_scores = np.nan_to_num(z_scores, nan=0.0, posinf=vmax, neginf=vmin) + return np.clip(z_scores, vmin, vmax) + + +def build_expression_matrix( + periodic_genes: pd.DataFrame, + time_columns: Sequence[int | float], + vmin: float = DEFAULT_VMIN, + vmax: float = DEFAULT_VMAX, +) -> npt.NDArray[np.float64]: + """Build the normalized heatmap matrix from ordered genes. + + Parameters + ---------- + periodic_genes: + Ordered periodic-gene table. + time_columns: + Time-point columns to extract. + vmin: + Lower clipping bound for z-scores. + vmax: + Upper clipping bound for z-scores. + + Returns + ------- + numpy.ndarray + Normalized expression matrix. + """ + expression_dataframe = periodic_genes.loc[:, list(time_columns)].apply( + pd.to_numeric, + errors="coerce", + ) + return compute_row_z_scores(expression_dataframe.to_numpy(dtype=float), vmin, vmax) + + +def make_heatmap_colormap() -> LinearSegmentedColormap: + """Create the cyan-black-yellow colormap used for the heatmap. + + Returns + ------- + matplotlib.colors.LinearSegmentedColormap + Colormap where low values are cyan and high values are yellow. + """ + return LinearSegmentedColormap.from_list( + "cyan_black_yellow", + ["#00d5e8", "#111111", "#ffea00"], + N=256, + ) + + +def plot_figure_2a( + z_scores: npt.NDArray[np.float64], + time_columns: Sequence[int | float], + output_path: Path, + *, + add_panel_label: bool = True, + dpi: int = DEFAULT_DPI, + vmin: float = DEFAULT_VMIN, + vmax: float = DEFAULT_VMAX, +) -> Figure: + """Plot and save the Figure 2A heatmap. + + Parameters + ---------- + z_scores: + Normalized expression matrix. + time_columns: + Time-point columns matching the matrix columns. + output_path: + Destination PNG path. + add_panel_label: + Whether to draw the large ``A.`` panel label. + dpi: + Figure resolution. + vmin: + Lower color-scale bound. + vmax: + Upper color-scale bound. + + Returns + ------- + matplotlib.figure.Figure + The generated Matplotlib figure. + + Raises + ------ + ValueError + If matrix and time-column dimensions do not match. + """ + if z_scores.ndim != 2: + msg = "z_scores must be a two-dimensional matrix." + raise ValueError(msg) + if z_scores.shape[1] != len(time_columns): + msg = "Number of matrix columns must match number of time columns." + raise ValueError(msg) + + output_path.parent.mkdir(parents=True, exist_ok=True) + number_of_genes = z_scores.shape[0] + + figure, axis = plt.subplots(figsize=(3.05, 5.55), dpi=dpi) + axis.imshow( + z_scores, + aspect="auto", + interpolation="nearest", + cmap=make_heatmap_colormap(), + vmin=vmin, + vmax=vmax, + origin="upper", + ) + + tick_positions = [ + list(time_columns).index(time) for time in DEFAULT_XTICK_TIMES if time in time_columns + ] + tick_labels = [str(time) for time in DEFAULT_XTICK_TIMES if time in time_columns] + axis.set_xticks(tick_positions) + axis.set_xticklabels(tick_labels) + axis.set_yticks([]) + + axis.set_xlabel("time (minutes)", fontsize=12, labelpad=13) + axis.set_ylabel(f"Top Periodic Genes ({number_of_genes})", fontsize=12, labelpad=10) + axis.set_title(r"$\it{Saccharomyces\ cerevisiae}$", fontsize=12, pad=8) + + for spine in axis.spines.values(): + spine.set_visible(True) + spine.set_linewidth(1.2) + spine.set_color("black") + + axis.tick_params(axis="x", labelsize=11, length=5, width=1) + axis.tick_params(axis="y", length=0) + + if add_panel_label: + axis.text( + -0.23, + 1.04, + "A.", + transform=axis.transAxes, + fontsize=30, + fontweight="bold", + va="top", + ha="left", + ) + + figure.tight_layout() + figure.savefig(output_path, bbox_inches="tight") + return figure + + +def reproduce_figure_2a( + input_path: Path, + output_path: Path, + *, + sheet_name: str = DEFAULT_SHEET_NAME, + header_row: int = DEFAULT_HEADER_ROW, + order_column: str = DEFAULT_ORDER_COLUMN, + add_panel_label: bool = True, +) -> Path: + """Run the complete Figure 2A reproduction workflow. + + Parameters + ---------- + input_path: + Path to the supplementary Excel file. + output_path: + Destination path for the PNG figure. + sheet_name: + Excel sheet to read. + header_row: + Zero-based row index containing column names. + order_column: + Column containing the Figure 2A order index. + add_panel_label: + Whether to draw the large ``A.`` panel label. + + Returns + ------- + pathlib.Path + Path to the generated figure. + """ + expression_table = read_expression_table(input_path, sheet_name, header_row) + time_columns = get_time_columns(expression_table) + periodic_genes = select_and_order_periodic_genes(expression_table, order_column) + z_scores = build_expression_matrix(periodic_genes, time_columns) + figure = plot_figure_2a( + z_scores, + time_columns, + output_path, + add_panel_label=add_panel_label, + ) + plt.close(figure) + return output_path + + +def build_argument_parser() -> argparse.ArgumentParser: + """Build the command-line argument parser. + + Returns + ------- + argparse.ArgumentParser + Configured parser. + """ + parser = argparse.ArgumentParser( + description="Reproduce Figure 2A from Kelliher et al. 2016.", + ) + parser.add_argument( + "--input", + type=Path, + default=DEFAULT_INPUT_PATH, + help="Path to pgen.1006453.s002.xlsx.", + ) + parser.add_argument( + "--output", + type=Path, + default=DEFAULT_OUTPUT_PATH, + help="Output PNG path.", + ) + parser.add_argument( + "--sheet-name", + default=DEFAULT_SHEET_NAME, + help="Excel sheet name to read.", + ) + parser.add_argument( + "--header-row", + type=int, + default=DEFAULT_HEADER_ROW, + help="Zero-based Excel header row index.", + ) + parser.add_argument( + "--no-panel-label", + action="store_true", + help="Do not draw the large A. label.", + ) + return parser + + +def main(argv: Sequence[str] | None = None) -> int: + """Run the command-line interface. + + Parameters + ---------- + argv: + Optional command-line arguments. ``None`` uses ``sys.argv``. + + Returns + ------- + int + Process exit status. + """ + parser = build_argument_parser() + args = parser.parse_args(argv) + output_path = reproduce_figure_2a( + input_path=args.input, + output_path=args.output, + sheet_name=args.sheet_name, + header_row=args.header_row, + add_panel_label=not args.no_panel_label, + ) + print(f"Figure saved to: {output_path.resolve()}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main())