diff --git a/.dockerignore b/.dockerignore index 6f3017e..c6dd0b9 100644 --- a/.dockerignore +++ b/.dockerignore @@ -64,10 +64,26 @@ compose.yaml # Data & Output files data/ renders/ +*_renders/ *.pdf *.csv *.xlsx +benchmark/ +benchmark_data/ # Project files (not needed in container) MANIFEST.in -requirements.txt \ No newline at end of file +requirements.txt + +src/cython_wrapper/*.c +src/cython_wrapper/*.so +*.egg-info +tests/ +.old_venv/ +.bioprofilekit +iedb_test_report.html +iedb_sub_renders.zip +profile.out +test.html +tests_ +__pycache__ \ No newline at end of file diff --git a/.gitignore b/.gitignore index eff2f70..7f294bf 100644 --- a/.gitignore +++ b/.gitignore @@ -284,3 +284,6 @@ Temporary Items src/renders/ /data/ +*_renders/ +*.zip +renders/ \ No newline at end of file diff --git a/Makefile b/Makefile index d16f694..1ff2c92 100644 --- a/Makefile +++ b/Makefile @@ -13,13 +13,14 @@ help: @echo "$$help" install: + python -m pip install --upgrade pip pip install -e . python setup.py build_ext --inplace clean: rm -rf build/ dist/ */*.egg-info/ - find . -name ".so" -delete - find . -name ".c" -delete + find src/cython_wrapper -name "*.so" -delete + find src/cython_wrapper -name "*.c" -delete find . -name "__pycache__" -type d -exec rm -rf {} + reinstall: clean install diff --git a/README.md b/README.md index 33be5c3..e732c5b 100644 --- a/README.md +++ b/README.md @@ -1,27 +1,65 @@ # BioProfileKit -[![Python: 3.13](https://img.shields.io/badge/Python-3.13-green.svg)](https://www.python.org/downloads/release/python-3130/) [![License: MIT](https://img.shields.io/badge/License-MIT-darkred.svg)](https://opensource.org/licenses/MIT) +[![Python: 3.13](https://img.shields.io/badge/Python-3.13-green.svg)](https://www.python.org/downloads/release/python-3130/) [![License: MIT](https://img.shields.io/badge/License-MIT-darkred.svg)](https://opensource.org/licenses/MIT) [![pages-build-deployment](https://github.com/BioWarlock/BioProfileKit/actions/workflows/pages/pages-build-deployment/badge.svg?branch=main)](https://github.com/BioWarlock/BioProfileKit/actions/workflows/pages/pages-build-deployment) ## Overview -BioProfileKit is a specialized bioinformatics tool that enables scientists to analyze large and diverse datasets. -Unlike traditional profilers, it offers customized analyzes for genomics, proteomics, transcriptomics, and metabolomics, providing sequence analysis and reports on nucleotide or amino acid distribution and abundance. -It includes advanced visualizations for pattern and anomaly detection, along with interactive dashboards showing key metrics. -Designed to be user-friendly, BioProfileKit is accessible to scientists without extensive data science skills. -Example Results: [BioProfileKit Example](https://hansen-maria.github.io/BioProfileKit/) +BioProfileKit is a domain-specific bioinformatics profiling tool for tabular biological data. It generates self-contained, interactive HTML reports covering exploratory data analysis, biological sequence profiling, and metadata validation — designed to be accessible to scientists without extensive data science skills. + +Example report: [BioProfileKit Example](https://hansen-maria.github.io/BioProfileKit/) ### Features -- **EDA-Focused Analysis** - Offers a detailed overview of data structure, quality, and composition with automated detection of issues like missing values or correlations -- **Specialized Sequence Profiling** - Examines DNA, RNA, and protein sequences using relevant metrics (such as GC content, k-mer frequencies, and amino acid composition) -- **Biological Metadata Recognition** - Automatically identifies and verifies organism names, taxonomic identifiers, and biological annotations using controlled vocabularies from official databases -- **Rich Visualizations** - Creates histograms and interactive charts to help quickly identify patterns -- **Interactive HTML Reports** - Provides portable, user-friendly reports with dynamic filtering and cross-linked visualizations for seamless data exploration -## WebApplication +**Column-level EDA** +- Automated column type inference (numeric, categorical, sequence, metadata) +- Distributional metrics and quality flags per column type +- Numeric quality metrics: IQR, outlier detection (adjusted IQR, modified Z-score, standard Z-score), zero/negative/infinity counts +- Categorical imbalance metrics: class imbalance ratio, rare category detection, near-zero variance, top-1/top-5 coverage +- Collapsible filter bar with text search and quality flag toggles (AND/OR logic) + +**Dataset-level EDA** +- Correlation and pairwise structure analysis +- Duplicate detection +- Missing value pattern analysis: MCAR (Little's test), MAR (point-biserial correlation), MNAR (KS-test heuristic) +- Interactive Plotly visualizations + +**Biological Sequence Profiling** +- DNA/RNA: GC content, nucleotide distribution, k-mer frequencies, AT/GC skew, ambiguous base ratio, dinucleotide O/E ratios, codon completeness, low-complexity entropy, reverse complement duplicate detection +- Protein: amino acid composition and distribution +- Configurable k-mer size and top-N reporting + +**Biological Metadata Validation** +- Organism name and taxonomic identifier verification (Cython-accelerated) +- Functional annotation validation: GO terms and COG categories +- Controlled vocabulary checks against official databases + +**Interactive HTML Reports** +- Portable, single-file reports with no server dependency +- Dynamic column filtering with quality flag toggles +- Colorblind-safe theme (Okabe-Ito palette) with persistent toggle +- WCAG AA accessible + +## Project Structure + +``` +BioProfileKit/ +├── src/ +│ ├── analysis/ # Column-level and dataset-level EDA +│ ├── biological/ # Sequence profiling and metadata validation +│ ├── cli/ # Entry point and report writer +│ ├── cython_wrapper/ # Cython-accelerated components +│ ├── data_utils/ # File I/O and remote data access +│ ├── models/ # Data models per column type +│ └── templates/ # Jinja templates for reports +└── tests/ # Unit and integration tests +``` + +## Web Application https://bioprofilekit.computational.bio/ + ## Installation ```bash @@ -34,30 +72,40 @@ pip install -e . python setup.py build_ext --inplace ``` -## Parameters -Currently only supports .csv, .tsv, and .json as input files -```bash - bioprofilekit -i input.csv - Options: - -i, --input PATH Input file as .tsv, .csv or .json [required] - -t, --tax Enable taxonomy analysis - -f, --func [cog|go] Choose between COG or GO analysis, if validation is needed - -tc, --target_column TEXT Target column for further analysis - -k, --kmer INTEGER K-mer Size for sequence analysis. Default: 3 - -n, --top_n INTEGER Top N entries analysis. Default: 20 - -h, --help Show this message and exit. + +## Usage + +Supports `.csv`, `.tsv`, and `.json` as input formats. + +```bash +bioprofilekit -i input.csv + +Options: + -i, --input PATH Input file (.tsv, .csv, or .json) [required] + -t, --tax Enable taxonomy analysis + -f, --func [cog|go] Functional annotation validation (COG or GO) + -tc, --target_column TEXT Target column for further analysis + -k, --kmer INTEGER K-mer size for sequence analysis [default: 3] + -n, --top_n INTEGER Top N entries for reporting [default: 20] + -h, --help Show this message and exit. ``` +> [!WARNING] +> For high-dimensional datasets with many columns, report rendering times may increase significantly. In such cases, the hosted web application at https://bioprofilekit.computational.bio/ is recommended. + ## Contributing -Contributions to this project are welcome! Whether you find bugs, want to request features, or submit enhancements, please feel free to open an issue or submit a pull request. For major changes, it's recommended to discuss them first to ensure alignment with project goals. -Please read the [`CODE OF CONDUCT`](CODE_OF_CONDUCT.md) to learn more about our guidelines and the contribution process. +Contributions are welcome. For bug reports, feature requests, or enhancements, please open an issue or submit a pull request. For major changes, open an issue first to discuss alignment with project goals. + +Please read the [`CODE OF CONDUCT`](CODE_OF_CONDUCT.md) for contribution guidelines. + ## License -Licensed under MIT license ([LICENSE-MIT](LICENSE) or http://opensource.org/licenses/MIT). -Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in BioProfileKit by you, shall be licensed as above, without any additional terms or conditions. +Licensed under the MIT license ([LICENSE](LICENSE) or https://opensource.org/licenses/MIT). +Unless explicitly stated otherwise, any contribution intentionally submitted for inclusion in BioProfileKit shall be licensed under the same terms. + ## Contact -For inquiries or support regarding this project, you can reach out to the maintainers through GitHub issues or Discussion. \ No newline at end of file +For inquiries or support, reach out via GitHub Issues or Discussions. diff --git a/__main__.py b/__main__.py index 58d3e25..3eb087a 100644 --- a/__main__.py +++ b/__main__.py @@ -1,9 +1,11 @@ -from app import cli +from src.cli.app import cli if __name__ == '__main__': - cli(['--input', 'data/iedb_sub.tsv', '--tax', '-tc','bind_class']) - #cli(['-i', "data/Complete_nitrogenase.csv"]) - #cli(['--input', 'data/upstream_sd.complete.tsv']) - #cli(['-i', "data/iedb_test.tsv"]) + #cli(['--input', 'data/iedb.tsv', '--tax', '-tc','bind_class']) + # cli(['-i', "data/Complete_nitrogenase.csv"]) + #cli(['--input', 'data/sorfdb_sub.tsv']) + # cli(['-i', "data/iedb_sub.tsv", "--tax"]) + #cli(['-i', 'data/winequality-white.csv']) # cli(['-i', "data/proteinGroups.tsv"]) - #cli(['-i', 'data/sample_bakrep.tsv']) \ No newline at end of file + cli(['-i', 'data/bakrep_sagalactiae.tsv']) + #cli(['-i', 'data/iedb_test.tsv', '-tc', 'bind_class']) \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 663ba06..5519467 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -2,8 +2,8 @@ requires = [ "setuptools==82.0.0", "wheel>=0.44.0", - "cython>=3.2.1", - "numpy>=2.4.2" + "cython>=3.2.5", + "numpy>=2.4.6" ] build-backend = "setuptools.build_meta" @@ -16,10 +16,9 @@ authors = [ {name = "Sonja Diedrich", email = "sonja.diedrich@computational.bio.uni-giessen.de"}, {name = "Maria Hansen", email = "maria.hansen@computational.bio.uni-giessen.de"} ] -#license = "MIT" -#license-files = ["LICENSE"] +license = {file = "LICENSE"} requires-python = ">=3.13" -keywords = ["bioinformatics", "data-analysis", "proteomics"] +keywords = ["bioinformatics", "data-analysis", "proteomics", "EDA"] classifiers = [ "Development Status :: 3 - Alpha", "Intended Audience :: Science/Research", @@ -29,22 +28,23 @@ classifiers = [ "Topic :: Scientific/Engineering :: Bio-Informatics", ] dependencies = [ - "cython==3.2.4", - "click>=8.3.1", + "cython==3.2.5", + "click>=8.3.3", "Jinja2==3.1.6", - "numpy==2.4.2", - "pandas==3.0.1", - "plotly==6.5.2", + "numpy==2.4.6", + "pandas==3.0.3", + "plotly==6.8.0", "scipy==1.17.1", "termcolor==3.3.0", - "pyarrow==23.0.0", + "pyarrow==24.0.0", "peptides==0.5.0", "iteration_utilities==0.13.0", - "requests==2.32.5", - "goatools==1.5.2", - "importlib_resources==6.5.2", - "biopython==1.86", + "requests==2.34.2", + "goatools==1.6.5", + "importlib_resources==7.1.0", + "biopython==1.87", "weblogo==3.9.0", + "scikit-learn==1.9.0" ] [project.urls] Homepage="https://github.com/BioWarlock/BioProfileKit/" @@ -52,6 +52,10 @@ Repository="https://github.com/BioWarlock/BioProfileKit/" [project.optional-dependencies] cli = ["click"] +dev = [ + "pytest>=9.1.0", + "pytest-cov>=6.0.0", +] [tool.setuptools] package-dir = {"" = "src"} @@ -60,7 +64,8 @@ include-package-data = true [tool.setuptools.package-data] "static" = ["static/**/*"] -"templates"= ["*.jinja"] +"templates" = ["*.jinja"] +#"cython" = ["*.pyx"] [tool.setuptools.exclude-package-data] "*" = [ @@ -69,5 +74,26 @@ include-package-data = true "*.egg-info" ] +[tool.pytest.ini_options] +testpaths = ["tests"] +pythonpath = ["src"] +addopts = "--tb=short -q" + +[tool.coverage.run] +source = ["src"] +omit = ["src/cython_wrapper/*"] + +[tool.coverage.report] +exclude_lines = [ + "pragma: no cover", + "if __name__ == .__main__.:", + "raise NotImplementedError", +] +show_missing = true +skip_empty = true + [project.scripts] -bioprofilekit = "app:cli" +bioprofilekit = "cli.app:cli" + +[pytest] +strict = true \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 75856df..f8caee8 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,6 @@ peptides~=0.5.0 weblogo~=3.7.12 -requests~=2.32.5 +requests~=2.34.2 goatools~=1.5.1 pandas~=2.3.0 jinja2~=3.1.6 @@ -11,4 +11,5 @@ termcolor~=3.1.0 plotly~=6.3.0 scipy~=1.16.1 setuptools~=70.0.0 -matplotlib~=3.10.7 \ No newline at end of file +matplotlib~=3.10.7 +pytest~=9.1.0 \ No newline at end of file diff --git a/setup.py b/setup.py index c9c4176..1d8a4a9 100644 --- a/setup.py +++ b/setup.py @@ -1,5 +1,5 @@ import sysconfig - +import numpy as np from Cython.Build import cythonize from setuptools import setup, Extension, find_packages @@ -9,13 +9,16 @@ def get_python_include_dir(): extensions = [ Extension( - "qc_eda.basic.wrapper_utils", - ["src/qc_eda/basic/wrapper_utils.pyx"], + "cython_wrapper.wrapper_utils", + ["src/cython_wrapper/wrapper_utils.pyx"], ), Extension( - "qc_eda.basic.taxonomy_validator", - ["src/qc_eda/basic/taxonomy_validator.pyx"], + "cython_wrapper.taxonomy_validator", + ["src/cython_wrapper/taxonomy_validator.pyx"], ), + Extension("cython_wrapper.medcouple_fast", + ["src/cython_wrapper/medcouple_fast.pyx"], + include_dirs=[np.get_include()]) ] setup( @@ -23,6 +26,6 @@ def get_python_include_dir(): version="0.1", packages=find_packages(where="src"), package_dir={"": "src"}, - ext_modules=cythonize(extensions, compiler_directives={'language_level': "3"}), + ext_modules=cythonize(extensions, compiler_directives={'language_level': "3"}, force=True), zip_safe=False, ) \ No newline at end of file diff --git a/src/analysis/__init__.py b/src/analysis/__init__.py new file mode 100644 index 0000000..8a045f6 --- /dev/null +++ b/src/analysis/__init__.py @@ -0,0 +1,6 @@ +from .overview import overview, column_overview +from .numeric_analysis import numeric_columns +from .categorical_analysis import categorical_columns +from .multivariate_analysis import get_correlation, multivariate_analysis +from .data_quality import check_mixed_types, check_suspect_values +from .plot_utils import plot_overview \ No newline at end of file diff --git a/src/analysis/categorical_analysis.py b/src/analysis/categorical_analysis.py new file mode 100644 index 0000000..991e14c --- /dev/null +++ b/src/analysis/categorical_analysis.py @@ -0,0 +1,47 @@ +import numpy as np +import pandas as pd + +from models.categorical import CategoricalColumns + + +def categorical_columns(df: pd.DataFrame, col: str) -> CategoricalColumns: + value_counts = df[col].value_counts() + n = int(df[col].notna().sum()) + frequencies = value_counts / n if n > 0 else value_counts + rare = np.sum(value_counts <= len(df[col]) * 0.01) + cib_ratio = value_counts.max() / value_counts.min() + top_coverage = df[col].value_counts(normalize=True).head(5).sum() + + if n > 0: + entropy = -(frequencies * np.log2(frequencies)).sum() + gini = 1 - (frequencies ** 2).sum() + simpson = 1 / (frequencies ** 2).sum() + cardinality_ratio = df[col].nunique() / n + mode_value = df[col].mode().iloc[0] + else: + entropy = np.nan + gini = np.nan + simpson = np.nan + cardinality_ratio = np.nan + mode_value = "" + lengths = df[col].astype(str).str.len() + + return CategoricalColumns( + name=col, + unique_categories=df[col].nunique(), + mode=mode_value, + entropy=round(entropy, 2), + frequencies=df[col].value_counts(normalize=True).head(20).to_dict(), + gini=round(gini, 2), + simpson_diversity=round(simpson, 2), + value_counts=df[col].value_counts().head(20).to_dict(), + max_category_length=lengths.max(), + min_category_length=lengths.min(), + cardinality_ratio=round(cardinality_ratio, 3), + rare_categories=rare, + top_5_coverage=top_coverage, + cib_ratio=cib_ratio, + top_1_coverage=(df[col].value_counts(normalize=True).max() if n > 0 else np.nan), # Near Zero Variance (NZV) + effective_cardinality=((df[col].nunique()/len(df[col]))*100), + memory = df[col].memory_usage(deep=True), + ) diff --git a/src/analysis/data_quality.py b/src/analysis/data_quality.py new file mode 100644 index 0000000..4258b34 --- /dev/null +++ b/src/analysis/data_quality.py @@ -0,0 +1,62 @@ +import numpy as np +import pandas as pd +from pandas.api.types import infer_dtype + + +def _split_numeric_string(series): + inferred = infer_dtype(series) + + if 'mixed' in inferred: + is_numeric = series.apply(lambda x: isinstance(x, (int, float, complex))) + return series[is_numeric], series[~is_numeric] + + if str(series.dtype) in ('str', 'string', 'object'): + numeric_mask = pd.to_numeric(series, errors='coerce').notna() & series.notna() + return series[numeric_mask], series[~numeric_mask] + + return None, None + + +def check_mixed_types(df, col, suspect_threshold=0.125): + series = df[col].dropna() + if series.empty: + return None + + numeric_part, string_part = _split_numeric_string(series) + if numeric_part is None or len(numeric_part) == 0 or len(string_part) == 0: + return None + + minority_count = min(len(numeric_part), len(string_part)) + minority_ratio = minority_count / len(series) + if len(numeric_part) >= len(string_part): + majority_type = "Numeric" + minority = string_part + else: + majority_type = "String" + minority = numeric_part + + if minority_ratio < suspect_threshold: + return majority_type, minority + else: + return majority_type, None + + +def check_suspect_values(df, col): + series = df[col].dropna() + if series.empty: + return None + + if str(series.dtype) in ('str', 'string', 'object'): + numeric_part, _ = _split_numeric_string(series) + if numeric_part is None or len(numeric_part) == 0: + return None + if len(numeric_part) / len(series) < 0.15: + return numeric_part + + if pd.api.types.is_numeric_dtype(series): + sentinels = {np.inf, -np.inf} + found = [v for v in sentinels if v in series.values] + if found: + return pd.Series([str(v) for v in found]) + + return None \ No newline at end of file diff --git a/src/analysis/multivariate_analysis.py b/src/analysis/multivariate_analysis.py new file mode 100644 index 0000000..66ca693 --- /dev/null +++ b/src/analysis/multivariate_analysis.py @@ -0,0 +1,528 @@ +import numpy as np +import pandas as pd +import plotly.express as px +import plotly.graph_objects as go +from scipy.sparse import coo_matrix +from scipy.stats import chi2_contingency + +from models.multivariate import MultivariateAnalysis +from .plot_utils import apply_standard_axes + + +def multivariate_analysis(df: pd.DataFrame, target: str) -> MultivariateAnalysis: + types = _classify_columns(df) + values, methods = compute_correlation_matrix(df) + feat_target_corr = feature_target_correlation(df, target) if target else None + mi = mutual_information(df, target) if target else None + + return MultivariateAnalysis( + correlation_heatmap=correlation_heatmap(values, methods), + pearson_heatmap=pearson_correlation_heatmap(df), + cramers_heatmap=cramers_heatmap(values, types), + eta_heatmap=eta_heatmap(values, types), + missing_matrix=missing_matrix(df), + missing_values_barchart=missing_values_barchart(df), + balance_plot=balance_plot(df, target) if target else None, + boxplot=boxplot(df), + scatter_matrix=scatter_matrix(df), + correlation_matrix=values, + correlation_methods=methods, + top_associations=top_associations(values, methods), + feature_target_correlation=feat_target_corr, + feature_target_plot=feature_target_plot(feat_target_corr) if feat_target_corr else None, + mutual_information=mi, + mutual_information_plot=mutual_information_plot(mi) if mi else None, + mi_relationship_plots=mi_relationship_plots(df, target, mi) if mi else None, + mcar_result=littles_mcar_test(df), + target_name=target if target else None, + ) + +def _classify_columns(df: pd.DataFrame) -> dict: + num = df.select_dtypes(include='number') + num_std = num.std(ddof=0) + numeric = num_std[num_std > 0].index + + cat = df.select_dtypes(exclude='number') + cat_nunique = cat.nunique(dropna=True) + categorical = cat_nunique[cat_nunique > 1].index + + return {**{c: 'numeric' for c in numeric}, **{c: 'categorical' for c in categorical}} + +""" +Correlation +""" +# Num x Num +def _pearson(df: pd.DataFrame, a: str, b: str) -> float: + pair = df[[a, b]].dropna() + if len(pair) < 2: + return np.nan + return pair[a].corr(pair[b], method='pearson') + +# Cat x Cat +def _cramers_v(df: pd.DataFrame, a: str, b: str) -> float: + """Bias-corrected Cramér's V (Bergsma 2013).""" + pair = df[[a, b]].dropna() + if pair.empty: + return np.nan + + a_codes, _ = pair[a].factorize() + b_codes, _ = pair[b].factorize() + if len(a_codes) == 0 or len(b_codes) == 0: + return np.nan + r = int(a_codes.max()) + 1 + k = int(b_codes.max()) + 1 + if r < 2 or k < 2: + return np.nan + + n = len(pair) + obs = coo_matrix((np.ones(n), (a_codes, b_codes)), shape=(r, k)).tocsr() + row_sum = np.asarray(obs.sum(axis=1)).ravel() + col_sum = np.asarray(obs.sum(axis=0)).ravel() + + coo = obs.tocoo() + expected = row_sum[coo.row] * col_sum[coo.col] / n + chi2 = np.sum(coo.data ** 2 / expected) - n + + phi2 = chi2 / n + r, k = obs.shape + + # Bergsma-Bias Correction + phi2corr = max(0.0, phi2 - ((k - 1) * (r - 1)) / (n - 1)) + rcorr = r - ((r - 1) ** 2) / (n - 1) + kcorr = k - ((k - 1) ** 2) / (n - 1) + + denom = min((kcorr - 1), (rcorr - 1)) + if denom <= 0: + return np.nan + return np.sqrt(phi2corr / denom) + +# Num x Cat +def _eta_squared(df: pd.DataFrame, cat: str, num: str) -> float: + """Correlation ratio eta^2 (ANOVA): variance in num explained by groups of cat.""" + pair = df[[cat, num]].dropna() + if len(pair) < 2: + return np.nan + + grand_mean = pair[num].mean() + agg = pair.groupby(cat)[num].agg(['count', 'mean']) + counts = agg['count'].to_numpy() + means = agg['mean'].to_numpy() + + ss_between = np.sum(counts * (means - grand_mean) ** 2) + ss_total = ((pair[num] - grand_mean) ** 2).sum() + + if ss_total == 0: + return np.nan + return ss_between / ss_total + + +def _correlation_pair(df, a, b, types) -> tuple[float, str]: + ta, tb = types[a], types[b] + + if ta == 'numeric' and tb == 'numeric': + return abs(_pearson(df, a, b)), 'Pearson' + if ta == 'categorical' and tb == 'categorical': + return _cramers_v(df, a, b), "Cramér's V" + cat, num = (a, b) if ta == 'categorical' else (b, a) + return _eta_squared(df, cat, num), 'Eta²' + + +def compute_correlation_matrix(df: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]: + types = _classify_columns(df) + cols = list(types.keys()) + n = len(cols) + + values = pd.DataFrame(np.eye(n), index=cols, columns=cols) + methods = pd.DataFrame('', index=cols, columns=cols, dtype=object) + for c in cols: + methods.at[c, c] = 'self' + + for i in range(n): + for j in range(i + 1, n): + a, b = cols[i], cols[j] + value, method = _correlation_pair(df, a, b, types) + value = round(float(value), 3) if value == value else np.nan + + values.iat[i, j] = value + values.iat[j, i] = value + methods.iat[i, j] = method + methods.iat[j, i] = method + + return values, methods + + +def feature_target_correlation(df: pd.DataFrame, target: str) -> dict | None: + if target not in df.columns: + return None + + types = _classify_columns(df) + if target not in types: + return None + + result = {} + for col in types: + if col == target: + continue + value, method = _correlation_pair(df, col, target, types) + if value == value: + result[col] = {'value': round(float(value), 3), 'method': method} + return result or None + +def feature_target_plot(feature_target: dict | None): + if not feature_target: + return None + + ranked = sorted(feature_target.items(), key=lambda x: x[1]['value']) + features = [f for f, _ in ranked] + valuevec = [info['value'] for _, info in ranked] + methodvec = [info['method'] for _, info in ranked] + + method_colors = { + "Pearson": "#0F65A0", + "Cramér's V": "#65A1E1", + "Eta squared": "#994564", + } + bar_colors = [method_colors.get(m, "#A1ACBD") for m in methodvec] + + fig = go.Figure() + seen = set() + for f, v, m, c in zip(features, valuevec, methodvec, bar_colors): + fig.add_trace(go.Bar( + x=[v], y=[f], orientation='h', + marker_color=c, name=m, + legendgroup=m, showlegend=m not in seen, + text=[f"{v:.3f}"], textposition="outside", + hovertemplate=f"{f}
Association: {v:.3f}
Method: {m}", + )) + seen.add(m) + + fig.update_layout( + title="Feature–Target Association", + xaxis=dict(title="Association strength", range=[0, 1.05]), + yaxis=dict(title="Feature"), + template="plotly_white", + bargap=0.3, + height=max(400, 40 * len(features) + 150), + legend=dict(title="Method"), + ) + return fig.to_html(full_html=False, include_plotlyjs=False) + +#ToDo change for Column Overview +def get_correlation(df: pd.DataFrame, col) -> list | None: + ncols = df.select_dtypes(include='number').dropna(axis=1, how='all').columns + if col in ncols: + std = df[ncols].std(ddof=0) + ncols = std[std > 0].index + if col not in ncols: + return None + corr = df[ncols].corrwith(df[col], method='pearson') + corr.drop(labels=col, inplace=True) + corr = corr[corr.abs() >= 0.3] + if corr.empty: + return None + return list(zip(corr.index, corr)) + return None + +def pearson_correlation_heatmap(df: pd.DataFrame): + df_numeric = df.select_dtypes(include=['float64', 'int64']).dropna(axis=1, how='all') + if not df_numeric.empty: + std = df_numeric.std(ddof=0) + df_numeric = df_numeric.loc[:, std[std > 0].index] + corr_matrix = round(df_numeric.corr(), 3) + fig = px.imshow(corr_matrix, text_auto=True, + labels=dict(color="Correlation"), + color_continuous_scale="RdBu_r", aspect="auto", height=700) + fig.update_layout(title="Correlation Heatmap") + return fig.to_html(full_html=False, include_plotlyjs=False) + +def correlation_heatmap(df: pd.DataFrame, methods: pd.DataFrame): + fig = px.imshow( + df, text_auto=True, + labels=dict(color="Association"), + color_continuous_scale="Blues", aspect="auto", height=700, + zmin=0, zmax=1, + ) + fig.update_traces( + customdata=methods.values, + hovertemplate="%{x} ↔ %{y}
Association: %{z}
Method: %{customdata}", + ) + fig.update_layout( + title="Association Heatmap (mixed-type)", + plot_bgcolor="#FFFFFF", + margin=dict(b=120), + ) + return fig.to_html(full_html=False, include_plotlyjs=False) + + +def cramers_heatmap(df: pd.DataFrame, types: dict) -> str | None: + """Cramér's V matrix over categorical columns only (0..1).""" + cols = [c for c, t in types.items() if t == 'categorical' and c in df.index] + if len(cols) < 2: + return None + + mat = df.loc[cols, cols] + + colorscale = [ + [0.0, "#EAF4FB"], [0.25, "#A9E1FF"], [0.5, "#65A1E1"], + [0.75, "#4082C0"], [1.0, "#0F65A0"], + ] + + fig = px.imshow( + mat, text_auto=".3g", + labels=dict(color="Cramér's V"), + color_continuous_scale=colorscale, aspect="auto", + height=700, zmin=0, zmax=1, + ) + fig.update_traces( + hovertemplate="%{x} ↔ %{y}
Cramér's V: %{z}", + ) + fig.update_layout( + title="Cramér's V (categorical pairs)", + autosize=True, + xaxis=dict(tickangle=-45), + plot_bgcolor="#FFFFFF", + margin=dict(b=120), + ) + return fig.to_html(full_html=False, include_plotlyjs=False, config={"responsive": True}) + + +def eta_heatmap(df: pd.DataFrame, types: dict) -> str | None: + cats = [c for c, t in types.items() if t == 'categorical' and c in df.index] + nums = [c for c, t in types.items() if t == 'numeric' and c in df.index] + if not cats or not nums: + return None + + mat = df.loc[cats, nums] + + colorscale = [ + [0.0, "#F7E9F0"], [0.25, "#D27897"], [0.5, "#A83665"], + [0.75, "#994564"], [1.0, "#932263"], + ] + + fig = px.imshow( + mat, text_auto=".3g", + labels=dict(color="Eta²", x="Numeric", y="Categorical"), + color_continuous_scale=colorscale, aspect="auto", + height=700, zmin=0, zmax=1, + ) + fig.update_traces( + hovertemplate="%{y} (cat) → %{x} (num)
Eta²: %{z}", + ) + fig.update_layout( + title="Eta² — variance in numeric explained by categorical (directional)", + autosize=True, + xaxis=dict(tickangle=-45), + plot_bgcolor="#FFFFFF", + margin=dict(b=120), + ) + return fig.to_html(full_html=False, include_plotlyjs=False, config={"responsive": True}) + +def top_associations(values, methods, threshold=0.7): + pairs = [] + cols = list(values.columns) + for i in range(len(cols)): + for j in range(i + 1, len(cols)): + v = values.iat[i, j] + if v == v and v >= threshold: + pairs.append({ + "var1": cols[i], + "var2": cols[j], + "value": round(float(v), 3), + "method": methods.iat[i, j], + }) + return sorted(pairs, key=lambda p: p["value"], reverse=True) or None + + +def mutual_information(df: pd.DataFrame, target: str) -> dict | None: + from sklearn.feature_selection import mutual_info_classif, mutual_info_regression + + if target not in df.columns: + return None + types = _classify_columns(df) + if target not in types: + return None + + features = [c for c in types if c != target] + if not features: + return None + + work = df[features + [target]].dropna() + if len(work) < 3: + return None + + X = work[features].copy() + discrete_features = [] + for col in features: + if types[col] == 'categorical': + X[col], _ = X[col].factorize() + discrete_features.append(True) + else: + X[col] = X[col].astype(float) + discrete_features.append(False) + + if types[target] == 'categorical': + y, _ = work[target].factorize() + mi = mutual_info_classif(X, y, discrete_features=discrete_features, random_state=0) + else: + y = work[target].astype(float) + mi = mutual_info_regression(X, y, discrete_features=discrete_features, random_state=0) + + return {feat: {'value': round(float(score), 4)} for feat, score in zip(features, mi)} or None + + +def mutual_information_plot(mi: dict | None): + if not mi: + return None + + ranked = sorted(mi.items(), key=lambda x: x[1]['value']) + features = [f for f, _ in ranked] + values = [info['value'] for _, info in ranked] + + fig = go.Figure(go.Bar( + x=values, y=features, orientation='h', + marker_color="#0F65A0", + text=[f"{v:.3f}" for v in values], textposition="outside", + hovertemplate="%{y}
Mutual Information: %{x:.4f}", + )) + fig.update_layout( + title="Mutual Information (feature → target, non-linear)", + xaxis=dict(title="Mutual Information"), + yaxis=dict(title="Feature"), + template="plotly_white", + bargap=0.3, + height=max(400, 40 * len(features) + 150), + ) + return fig.to_html(full_html=False, include_plotlyjs=False) + + +def mi_relationship_plots(df: pd.DataFrame, target: str, mi: dict | None, top_n: int = 3): + if not mi or target not in df.columns: + return None + + types = _classify_columns(df) + if target not in types: + return None + + top = sorted(mi.items(), key=lambda x: x[1]['value'], reverse=True)[:top_n] + target_is_cat = types[target] == 'categorical' + plots = [] + for feature, info in top: + work = df[[feature, target]].dropna() + if work.empty: + continue + + if target_is_cat: + # distribution of the feature per target class + if types[feature] == 'numeric': + fig = px.violin( + work, x=target, y=feature, color=target, box=True, points=False, + color_discrete_sequence=["#0F65A0", "#994564", "#65A1E1", "#A83665"], + ) + else: + # categorical feature vs categorical target -> grouped counts + ct = work.groupby([feature, target]).size().reset_index(name="count") + fig = px.bar(ct, x=feature, y="count", color=target, barmode="group", + color_discrete_sequence=["#0F65A0", "#994564", "#65A1E1", "#A83665"]) + fig.update_layout(title=f"{feature} by {target} (MI={info['value']})") + else: + # numeric target: scatter + LOWESS trend + if types[feature] == 'numeric': + fig = px.scatter( + work, x=feature, y=target, trendline="lowess", + trendline_color_override="#994564", + color_discrete_sequence=["#0F65A0"], opacity=0.5, + ) + fig.update_traces(marker=dict(size=4), selector=dict(mode="markers")) + else: + # categorical feature vs numeric target -> box per category + fig = px.box(work, x=feature, y=target, + color_discrete_sequence=["#0F65A0"]) + fig.update_layout(title=f"{feature} vs {target} (MI={info['value']})") + + fig.update_layout(template="plotly_white", height=400) + plots.append({ + "feature": feature, + "value": info['value'], + "plot": fig.to_html(full_html=False, include_plotlyjs=False), + }) + return plots or None + +""" + +Missingness + +""" +def littles_mcar_test(df: pd.DataFrame): + return {"Hello":"World"} + +def missing_matrix(df: pd.DataFrame): + missing_values = df.isnull().astype(int) + fig = px.imshow(missing_values, labels=dict(color="Missing Values"), + aspect="auto", color_continuous_scale="blues_r", + title="Missing Values Matrix") + fig.update_yaxes(autorange='reversed') + fig.update_layout(coloraxis_showscale=False) + return fig.to_html(full_html=False, include_plotlyjs=False) + + +def missing_values_barchart(df: pd.DataFrame): + missing_counts = df.isna().sum() + fig = px.bar(x=missing_counts.index, y=missing_counts.values, + labels={'x': 'Columns', 'y': 'Missing Values'}, + color_discrete_sequence=['#0F65A0']) + fig.update_layout(title="Missing Values per Column", bargap=0.2, plot_bgcolor='white') + apply_standard_axes(fig, tick_angle=-45) + return fig.to_html(full_html=False, include_plotlyjs=False) + + +def balance_plot(df, target): + fig = px.histogram(df, x=target, color_discrete_sequence=["#0F65A0"], text_auto=True) + fig.update_layout(title="Class Balance (Target Distribution)", bargap=0.2, plot_bgcolor="white") + apply_standard_axes(fig) + return fig.to_html(full_html=False, include_plotlyjs=False) + + +def boxplot(df: pd.DataFrame): + df = df.select_dtypes(include=['float64', 'int64']).dropna(axis=1, how='all') + fig = go.Figure() + for col in df: + fig.add_trace(go.Box(y=df[col].values, name=df[col].name)) + fig.update_yaxes(type="log", title="Logarithmic", showticklabels=False) + fig.update_layout( + xaxis=dict(rangeslider=dict(visible=True), type="linear"), + title="Boxplot", xaxis_title="Columns", + yaxis_title="Values", legend_title="Columns", + ) + return fig.to_html(full_html=False, include_plotlyjs=False) + + +def scatter_matrix(df: pd.DataFrame): + df = df.select_dtypes(include=['float64', 'int64']).dropna(axis=1, how='all') + fig = px.scatter_matrix(df, color_discrete_sequence=["#0F65A0"], height=750) + fig.update_traces(diagonal_visible=False, + marker=dict(size=2, opacity=0.5, color="#0F65A0")) + fig.update_layout(title="Scatter Matrix", xaxis_title="Columns", + plot_bgcolor="white", bargap=0.2, dragmode="select") + + n_vars = len(df.columns) + for i in range(1, n_vars + 1): + for j in range(1, n_vars + 1): + if i == 1: + fig.update_layout({f'xaxis{j if j > 1 else ""}': dict( + mirror=True, ticks="outside", showline=True, + linecolor="black", linewidth=1, gridcolor="lightgrey", + title_standoff=25, + )}) + if j == 1: + fig.update_layout({f'yaxis{i if i > 1 else ""}': dict( + mirror=True, ticks="outside", showline=True, + linecolor="black", linewidth=1, gridcolor="lightgrey", + title_standoff=25, + )}) + + fig.for_each_annotation(lambda a: a.update( + textangle=0, + x=-0.08 if a.textangle == 90 or 'y' in str(a.yref) else a.x, + y=-0.08 if a.textangle == 0 and 'x' in str(a.xref) else a.y, + )) + return fig.to_html(full_html=False, include_plotlyjs=False) \ No newline at end of file diff --git a/src/analysis/numeric_analysis.py b/src/analysis/numeric_analysis.py new file mode 100644 index 0000000..2696322 --- /dev/null +++ b/src/analysis/numeric_analysis.py @@ -0,0 +1,73 @@ +import numpy as np +import pandas as pd +from scipy import stats + +from analysis.outlier_detection import detect_outliers +from models.numeric import NumericColumns + + +def _safe_round(val, decimals=2): + return round(val, decimals) if np.isfinite(val) else np.nan + + +def numeric_columns(df: pd.DataFrame, col) -> NumericColumns: + series = df[col].dropna() + + if series.empty: + coefficient_of_variation = np.nan + quantiles = np.array([np.nan, np.nan, np.nan]) + mode_value = np.nan + value_counts = {} + frequencies = {} + mad_value = np.nan + min_value = np.nan + max_value = np.nan + mean_value = np.nan + median_value = np.nan + std_value = np.nan + sum_value = np.nan + kurtosis_value = np.nan + skewness_value = np.nan + outliers = None + else: + mean_value = float(series.mean()) + std_value = float(series.std()) + coefficient_of_variation = np.nan if mean_value == 0 else float(std_value / abs(mean_value)) + quantiles = np.array(series.quantile([0.25, 0.5, 0.75]).to_list(), dtype=float) + mode_series = series.mode() + mode_value = float(mode_series.iloc[0]) if not mode_series.empty else np.nan + value_counts = series.value_counts().head(20).to_dict() + frequencies = series.value_counts(normalize=True).head(20).to_dict() + mad_value = float(stats.median_abs_deviation(series, nan_policy='omit')) + min_value = float(series.min()) + max_value = float(series.max()) + median_value = float(series.median()) + sum_value = float(series.sum()) + kurtosis_value = float(series.kurtosis()) + skewness_value = float(series.skew()) + outliers = detect_outliers(series.to_numpy(dtype=np.double)) + + return NumericColumns( + name=col, + min=round(min_value, 2) if np.isfinite(min_value) else np.nan, + max=round(max_value, 2) if np.isfinite(max_value) else np.nan, + mean=round(mean_value, 2) if np.isfinite(mean_value) else np.nan, + median=round(median_value, 2) if np.isfinite(median_value) else np.nan, + mode=round(mode_value, 2) if np.isfinite(mode_value) else np.nan, + std=round(std_value, 2) if np.isfinite(std_value) else np.nan, + sum=round(sum_value, 2) if np.isfinite(sum_value) else np.nan, + kurtosis=round(kurtosis_value, 2) if np.isfinite(kurtosis_value) else np.nan, + skewness=round(skewness_value, 2) if np.isfinite(skewness_value) else np.nan, + coefficient_of_variation=round(coefficient_of_variation, 2) if np.isfinite(coefficient_of_variation) else np.nan, + mad=mad_value if np.isfinite(mad_value) else np.nan, + quantiles=quantiles, + infinity= series.isin([np.inf, -np.inf]).sum(), + negative_count=np.sum((series < 0).values.ravel()), + zero_count= np.sum(series== 0), + memory=df[col].memory_usage(deep=True), + value_counts=value_counts, + frequencies=frequencies, + outliers=outliers, + ) + + diff --git a/src/analysis/outlier_detection.py b/src/analysis/outlier_detection.py new file mode 100644 index 0000000..d857b45 --- /dev/null +++ b/src/analysis/outlier_detection.py @@ -0,0 +1,64 @@ +from typing import Optional + +import numpy as np + +from cython_wrapper.medcouple_fast import fast_medcouple +from models.outliers import Outliers + +def detect_outliers(values: np.ndarray) -> Optional[Outliers]: + if len(values) < 20 or np.unique(values).size <= 10: + return None + + q1 = np.percentile(values, 25, axis=0) + q3 = np.percentile(values, 75, axis=0) + iqr = q3 - q1 + mean = np.mean(values) + std = np.std(values, ddof=1) + median = np.median(values) + mad = np.median(np.abs(values - median)) + + if iqr <= 0: + return None + # Adjusted IQR (Hubert & Vandervieren, 2008) + """ + Adjusted implementation adapted from robpy (v.0.0.6) + Source: https://robpy.readthedocs.io/en/latest/_modules/robpy/univariate/adjusted_boxplot.html + """ + mc = fast_medcouple(values) + + if mc >= 0: + lower_iqr = q1 - 1.5 * np.exp(-4 * mc) * iqr + upper_iqr = q3 + 1.5 * np.exp(3 * mc) * iqr + else: + lower_iqr = q1 - 1.5 * np.exp(-3 * mc) * iqr + upper_iqr = q3 + 1.5 * np.exp(4 * mc) * iqr + + lower_iqr = max(lower_iqr, values.min()) + upper_iqr = min(upper_iqr, values.max()) + + # Modified Z-Score (MAD-based, Iglewicz & Hoaglin, 1993) + mz_lower, mz_upper = 0, 0 + if mad > 0: + modified_z = 0.6745 * (values - median) / mad + mz_lower = int((modified_z < -3.5).sum()) + mz_upper = int((modified_z > 3.5).sum()) + + # Standard Z-Score + z_upper, z_lower = 0, 0 + if std > 0: + z = (values - mean) / std + z_lower = int((z < -3.0).sum()) + z_upper = int((z > 3.0).sum()) + + return Outliers( + lower_bound=np.round(lower_iqr, 4).astype(np.float64), + upper_bound=np.round(upper_iqr, 4).astype(np.float64), + n_lower_iqr=int((values < lower_iqr).sum()), + n_upper_iqr=int((values > upper_iqr).sum()), + medcouple=np.round(mc, 4), + n_lower_mzscore=mz_lower, + n_upper_mzscore=mz_upper, + n_lower_zscore=z_lower, + n_upper_zscore=z_upper, + ) + diff --git a/src/analysis/overview.py b/src/analysis/overview.py new file mode 100644 index 0000000..19586a9 --- /dev/null +++ b/src/analysis/overview.py @@ -0,0 +1,99 @@ +import numpy as np +import pandas as pd + +from biological.sequence_detection import check_sequence +from models.overview import DatasetSummary, ColumnOverview +from .data_quality import check_mixed_types, check_suspect_values +from .plot_utils import plot_overview +from .multivariate_analysis import get_correlation + +# Keywords that strongly suggest a taxonomy/organism column +_TAX_NAME_HINTS = frozenset({ + "organism", "species", "taxon", "taxonomy", "host", "strain", + "genus", "family", "phylum", "class", "order", "lineage", + "scientific_name", "scientificname", "tax", "ncbi", +}) + + +def _check_taxonomy_candidate(df: pd.DataFrame, col: str) -> bool: + """Lightweight heuristic: is this column likely to contain organism names? + + No network calls, no external data. Returns True when: + - The column is a string/object dtype (numeric columns can't be organism names), AND + - Either the column name contains a known taxonomy keyword, OR + the values look like binomial nomenclature (two capitalised words, e.g. "Homo sapiens"). + + Intentionally conservative — false negatives are better than flooding every + string column with a "Not Checked" badge. + """ + if not pd.api.types.is_object_dtype(df[col]): + return False + + col_lower = col.lower().replace(" ", "_") + if any(hint in col_lower for hint in _TAX_NAME_HINTS): + return True + + # Value-pattern heuristic: sample up to 50 non-null values and check + # whether a meaningful fraction looks like "Genus species [strain]" + sample = df[col].dropna().astype(str).head(50) + if sample.empty: + return False + + binomial_matches = sample.str.match( + r'^[A-Z][a-z]+ [a-z]' # "Homo s…" or "Escherichia c…" + ).sum() + return binomial_matches / len(sample) >= 0.5 + + +def overview(df: pd.DataFrame, file) -> DatasetSummary: + return DatasetSummary( + filename=file, + rows=df.shape[0], + cols=df.shape[1], + nulls=sum(df.isnull().sum()), + nulls_percentage=100 if df.isnull().all().all() else round(sum(df.isnull().sum()) * 100 / df.size, 2), + empty_rows=df.replace("", np.nan).isna().all(axis=1).sum(), + dup_row=int(df.duplicated().sum()), + dup_col=int(df.columns.duplicated().sum()), + ratio=round(df.shape[0] / df.shape[1], 3), + memory=int(df.memory_usage(deep=True).sum()), + alerts=0, + ) + + +def column_overview(df: pd.DataFrame, col) -> ColumnOverview: + seq_type, invalid = check_sequence(df, col) + mixed = check_mixed_types(df, col) if seq_type == "None" else None + if mixed is not None: + suspect = mixed[1] + else: + suspect = check_suspect_values(df, col) + missing_percentage = 100 if df[col].isnull().all() else round(sum(df[col].isnull()) * 100 / df[col].size, 2) + + # Only run the candidate check for non-sequence string columns — + # sequence columns are already classified; numeric/bool can't be taxonomy. + taxonomy_candidate = ( + _check_taxonomy_candidate(df, col) + if seq_type == "None" + else False + ) + + return ColumnOverview( + name=col, + number=int(df[col].notnull().sum()), + unique=df[col].nunique(), + missing=int(df[col].isnull().sum()), + missing_per=missing_percentage, + density=(100 - missing_percentage), + type=str(df[col].dtype), + sequence=seq_type, + invalid_seqs=invalid, + mixed_types=mixed, + suspect_values=suspect, + describe_plot=None if df[col].isnull().all() else plot_overview(df[col]), + constant=(df[col].nunique() == 1), + correlation=get_correlation(df, col), + cardinality_dimension_ratio=round(df[col].nunique() / len(df), 3), + monotonicity=True if df[col].is_monotonic_increasing or df[col].is_monotonic_decreasing else False, + taxonomy_candidate=taxonomy_candidate, + ) diff --git a/src/analysis/plot_utils.py b/src/analysis/plot_utils.py new file mode 100644 index 0000000..7c77209 --- /dev/null +++ b/src/analysis/plot_utils.py @@ -0,0 +1,31 @@ +import plotly.express as px + + +def apply_standard_axes(fig, tick_angle=None): + x_kwargs = dict(mirror=True, ticks='outside', showline=True, + linecolor='black', gridcolor='lightgrey') + if tick_angle is not None: + x_kwargs['tickangle'] = tick_angle + fig.update_xaxes(**x_kwargs) + fig.update_yaxes(mirror=True, ticks='outside', showline=True, + linecolor='black', gridcolor='lightgrey') + +#ToDo: Check if we want to keep nan in count plot. +def plot_overview(col): + col = col.dropna() + if col.nunique() == 1: + return " " #ToDo find better solution + if col.dtype in ('str', 'string'): + truncated = col.apply(lambda x: str(x)[:20] + '…' if len(str(x)) > 20 else str(x)) + if col.dtype != 'object': + bins = col.nunique() if col.nunique() < 20 else 20 + if col.dtype in ('str', 'string'): + fig = px.histogram(x=truncated, color_discrete_sequence=['#0F65A0']) + apply_standard_axes(fig, tick_angle=-45) + else: + fig = px.histogram(x=col, color_discrete_sequence=['#0F65A0'], nbins=bins) + apply_standard_axes(fig) + fig.update_layout(bargap=0.4, plot_bgcolor='white', xaxis_title=col.name,yaxis_title='Count', autosize=True) + fig.layout.xaxis.automargin = True + return fig.to_html(full_html=False, include_plotlyjs=False) + return None \ No newline at end of file diff --git a/src/app.py b/src/app.py deleted file mode 100644 index 6bf425f..0000000 --- a/src/app.py +++ /dev/null @@ -1,115 +0,0 @@ -#!/usr/bin/env python3 - -import shutil -from pathlib import Path - -import click -import pandas as pd -from jinja2 import Environment, FileSystemLoader -from termcolor import colored -from importlib_resources import files -from qc_eda.basic.numerical_data import overview, column_overview, numeric_columns, categorical_columns -from qc_eda.biological.biological_data import dna_rna_columns, protein_columns -from qc_eda.biological.functional_annotation import annotation_flags -from qc_eda.biological.measurement_data import measurement_columns -from qc_eda.biological.taxonomy import taxonomy_flags -from utils.download_metadata import get_tax_ids -from utils.file_reader import read_file, parse_parquet -from qc_eda.basic.general import general_plots - -CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) -TEMPLATE_DIR = files("templates").joinpath() -STATIC_DIR = files("static").joinpath() -env = Environment(loader=FileSystemLoader(str(TEMPLATE_DIR)), autoescape=True) - - -@click.group(context_settings=CONTEXT_SETTINGS, invoke_without_command=True) -@click.option("-i", "--input", type=click.Path(exists=True, resolve_path=True), required=True, - help="Input file as .tsv, .csv or .json", ) -@click.option('-t','--tax', is_flag=True,help='Enable taxonomy analysis') -@click.option('-f', '--func', type=click.Choice(['cog','go']), help='Enable functional annotation analysis. Choose between cog or go') -@click.option('-tc', '--target_column', type=str, help='Target column for Analysis') -@click.option('-k', '--kmer', type=int, default=3, help="K-mer Size for sequence analysis") -@click.option('-n', '--top_n', type=int, default=20, help="Top N entries analysis") -def cli(input: str, tax: bool = False, func: str = None, target_column: str = None, kmer: int = None, top_n: int = None): - input_path = Path(input) - print(colored(f'Reading file {input_path.name}', 'green')) - df = read_file(input_path) if input_path.suffix != '.parquet' else parse_parquet(input_path) - general = overview(df, input_path.name) - plots = general_plots(df, target_column) - - dups = df[df.duplicated(keep=False)] - dups = dups.reset_index() - duplicates_table = dups.to_html(classes="table table-hover table-responsive nowrap", border="0", - table_id="dup_table", - index=False) - print(colored(f'Analyse {len(df.columns)} columns', 'blue')) - - column_overviews = [column_overview(df, col) for col in df.columns] - tax_df = None - if tax: - tax_df = get_tax_ids() - for col_overview in column_overviews: - if tax and tax_df is not None: - col_overview.taxonomy = taxonomy_flags(df, col_overview.name, tax_df) - if func and tax_df is not None: - col_overview.annotation = [annotation_flags(df, col_overview.name, func)] - if hasattr(col_overview, "top_10") and isinstance(col_overview.top_10, pd.Series): - col_overview.top_10_items = list(col_overview.top_10.items()) - if col_overview.sequence == 'dna': - print(colored(f'Analyzing DNA/RNA sequences in column: {col_overview.name}', 'cyan')) - bio_data = dna_rna_columns(df[col_overview.name], k=kmer, top_n=top_n) - col_overview.dna_rna_data = bio_data - elif col_overview.sequence == 'protein': - print(colored(f'Analyzing protein sequences in column: {col_overview.name}', 'cyan')) - bio_data = protein_columns(df[col_overview.name], k=kmer, top_n=top_n) - col_overview.protein_data = bio_data - else: - col_overview.dna_rna_data = None - col_overview.protein_data = None - - if col_overview.sequence == 'None': - measurement_data = measurement_columns(col_overview, df) - if measurement_data: - print(colored(f'Analyzing lab measurements in column: {col_overview.name}', 'cyan')) - col_overview.measurement_data = measurement_data - else: - col_overview.measurement_data = None - else: - col_overview.measurement_data = None - - print(colored(f'Analyse {len(df.select_dtypes(include="number").columns)} numeric columns ', 'blue')) - - #ToDo: check cat_columns - exclude_cols = list() - #numeric_overviews = [if df[col].isnull().all() numeric_columns(df, col) for col in df.select_dtypes(include='number').columns] - numeric_overviews = [numeric_columns(df, col) if not df[col].isnull().all() else exclude_cols.append(col) for col in - df.select_dtypes(include="number").columns] - cat_columns = [col for col in df.select_dtypes(include=['str','object', 'bool', 'int64', 'float64']).columns if - any(i.sequence == 'None' for i in column_overviews if i.name == col)] - print(colored(f'Analyse {len(cat_columns)} object columns ', 'blue')) - print(exclude_cols) - categorical_overviews = [categorical_columns(df, col) if col not in exclude_cols else None for col in cat_columns] - - output_path = Path(input_path.stem+"_renders") - output_path.mkdir(parents=True, exist_ok=True) - - shutil.copytree(str(STATIC_DIR), output_path.as_posix()+"/static/", dirs_exist_ok=True) - - landing_template = env.get_template('LandingPage.jinja') - numeric_template = env.get_template('numeric_overview.jinja') - columns = env.get_template('columns.jinja') - stats = env.get_template('general_statistics.jinja') - - print(colored('Writing report …', 'green')) - with open(output_path.as_posix()+"/index.html", 'w', encoding="utf-8") as output: - print(landing_template.render(), file=output) - - with open(output_path.as_posix()+"/numeric_data.html", "w",encoding="utf-8") as output: - print(numeric_template.render(general=general, dups=duplicates_table), file=output) - - with open(output_path.as_posix()+"/columns.html", "w",encoding="utf-8") as output: - print(columns.render(columns=column_overviews, overview=numeric_overviews, categorical=categorical_overviews, top_n=top_n), file=output) - - with open(output_path.as_posix()+"/general_statistics.html", "w",encoding="utf-8") as output: - print(stats.render(plots=plots), file=output) diff --git a/src/biological/__init__.py b/src/biological/__init__.py new file mode 100644 index 0000000..cba3d97 --- /dev/null +++ b/src/biological/__init__.py @@ -0,0 +1,5 @@ +from .sequence_detection import check_sequence +from .sequence_data import dna_rna_columns, protein_columns +from .taxonomy import taxonomy_flags, build_lookups +from .functional_annotation import annotation_flags +from .measurement_data import measurement_columns \ No newline at end of file diff --git a/src/qc_eda/biological/functional_annotation.py b/src/biological/functional_annotation.py similarity index 95% rename from src/qc_eda/biological/functional_annotation.py rename to src/biological/functional_annotation.py index 3a726d4..07c5817 100644 --- a/src/qc_eda/biological/functional_annotation.py +++ b/src/biological/functional_annotation.py @@ -1,6 +1,6 @@ from dataclasses import dataclass import pandas as pd -from utils.download_metadata import get_clusters_of_orthologous_groups, get_gene_ontology +from data_utils.remote_data import get_clusters_of_orthologous_groups, get_gene_ontology @dataclass class AnnotationFlags: diff --git a/src/biological/measurement_data.py b/src/biological/measurement_data.py new file mode 100644 index 0000000..1cb1a20 --- /dev/null +++ b/src/biological/measurement_data.py @@ -0,0 +1,50 @@ +import re +from typing import List + +import pandas as pd + +from models.measurement import UNITColumns +from enums.measurement_enum import MEASUREMENTS + + +def measurement_columns(col: pd.Series, name: str, col_type: str) -> UNITColumns | bool: + name_match = MEASUREMENTS.UNIT_IN_COL_TITLE.value.match(name) + if name_match: + unit: str = name_match.group(0).lstrip('[').rstrip(']') + return UNITColumns( + units=[unit], + unit_counts={unit: 1}, + with_measurement=False + ) + if col_type == 'str': + values = col.dropna().unique().astype(str).tolist() + if all(len(x) > 1 for x in values): + measurement_and_unit = match_units(values, MEASUREMENTS.UNIT_COLUMN.value) + if len(measurement_and_unit) == len(values): + measurement_passed = True + else: + all_values = col.dropna().astype(str).tolist() + all_measurement_and_unit = match_units(values, MEASUREMENTS.UNIT_COLUMN.value) + measurement_passed = len(all_measurement_and_unit) / len(all_values) > 0.95 + if measurement_passed: + return UNITColumns( + units=[unit.split(' ')[1] if ' ' in unit else unit for unit in measurement_and_unit], + unit_counts=col.value_counts().to_dict(), + with_measurement=any([has_number_and_unit(unit) for unit in measurement_and_unit]) + ) + return False + +def match_units(entries: List[str], regex: re.Pattern) -> List[str]: + units: List[str] = [] + for entry in entries: + measurement_and_unit = regex.fullmatch(entry) + if measurement_and_unit: + units.append(measurement_and_unit.group(0)) + return units + +def has_number_and_unit(value: str) -> bool: + special_units = {'1/s', '1/m', '1/M', '1/h', '1/min'} # Add new ones, if needed + if value in special_units: + return False + pattern = re.compile(r'^-?\d+\.?\d*\s*[a-zA-Z°/%]+$') + return bool(pattern.match(value) and value not in special_units) diff --git a/src/biological/plot_utils.py b/src/biological/plot_utils.py new file mode 100644 index 0000000..dcde478 --- /dev/null +++ b/src/biological/plot_utils.py @@ -0,0 +1,199 @@ +from pandas import DataFrame +from plotly import express as px, graph_objs as go, graph_objects as go +from plotly.subplots import make_subplots + +config = { + 'toImageButtonOptions': { + 'format': 'png', # one of png, svg, jpeg, webp + 'filename': None, + 'height': 1200, + 'width': 1600, + 'scale': 4 + } +} + + +def length_distribution(all_overview: DataFrame, unit: str = "bp"): + bins = 300 if all_overview.shape[0] > 10000 else all_overview.shape[0] + fig = px.histogram( + all_overview, + x="lengths", + title="Sequence Length Distribution Across All Sequences", + labels={"lengths": f"Sequence Length ({unit})"}, + color_discrete_sequence=['#0F65A0'], + nbins=bins + ) + fig.update_traces( + marker_line_color="#616D78", + marker_line_width=1, + opacity=0.85 + ) + fig.update_layout( + template="plotly_white", + xaxis_title=f"Sequence Length ({unit})", + yaxis_title="Sequence Count", + bargap=0.5 + ) + config['toImageButtonOptions']['filename'] = "length_distribution" + return fig.to_html(full_html=False, include_plotlyjs=False, config=config) + +def gc_distribution(all_overview: DataFrame): + fig = px.violin( + all_overview, + y="gc_content", + box=True, + points='all', + title="Distribution of GC Content Across All Sequences", + labels={"gc_content": "GC Content (%)"}, + ) + fig.update_traces( + fillcolor="#AB7E8F", + line_color="#616D78", + box_fillcolor="#994564", + box_line_color="#616D78", + marker=dict( + color="#0F65A0", + size=3, + ), + opacity=0.85, + line_width=1, + ) + fig.update_layout( + template="plotly_white", + yaxis=dict( + range=[0, 100], + ticksuffix="%" + ) + ) + config['toImageButtonOptions']['filename'] = "gc_distribution" + return fig.to_html(full_html=False, include_plotlyjs=False, config=config) + +def ambiguous_distribution(all_overview: DataFrame, col: str = "N", label: str = "N-Count"): + fig = go.Figure() + + fig.add_trace(go.Scattergl( + x=all_overview.index, + y=all_overview[col], # .sort_values(ascending=False), + mode='lines', + line=dict(color='#994564', width=0.85), + name=label + )) + + fig.update_layout( + title=f"Complete Distribution Curve of {label}s Across All Sequences", + xaxis_title="Sequence Index (Chronological Order)", + yaxis_title=f"Number of '{label}' in Sequence", + template="plotly_white", + yaxis=dict(rangemode="tozero") + ) + + return fig.to_html(full_html=False, include_plotlyjs=False) + + +def at_gc_skewness(all_overview: DataFrame): + fig = make_subplots( + rows=2, cols=2, + row_heights=[0.12, 0.88], + column_widths=[0.88, 0.12], + shared_xaxes=True, + shared_yaxes=True, + horizontal_spacing=0.02, + vertical_spacing=0.02 + ) + + fig.add_trace( + go.Histogram( + x=all_overview['GC_skew'], marker=dict( + color='#0F65A0', + line=dict( + color='#616D78', + width=0.5 + ) + ), + name='GC Skew Dist.', + hovertemplate='GC Skew: %{x}
Count: %{y}', + showlegend=False, + ), + row=1, col=1 + ) + + hover_labels = [] + for idx, seq in zip(all_overview.index, all_overview['sequence']): + if isinstance(seq, str) and seq.strip(): + short_seq = f"{seq[:20]}..." if len(seq) > 20 else seq + hover_labels.append(f"#{idx}: {short_seq}") + else: + hover_labels.append(f"#{idx}: [No Sequence Data]") + + all_overview['hover_sequence_label'] = hover_labels + + fig.add_trace( + go.Scatter( + x=all_overview['GC_skew'], + y=all_overview['AT_skew'], + mode='markers', + marker=dict(color='#994564', opacity=0.3), + customdata=all_overview[['hover_sequence_label']], + name='%{customdata[0]}', + hovertemplate='Seq: %{customdata[0]}
GC Skew: %{x}
AT Skew: %{y}', + showlegend=False + ), + row=2, col=1 + ) + + fig.add_trace( + go.Histogram( + y=all_overview['AT_skew'], + marker=dict( + color='#0F65A0', + line=dict( + color='#616D78', + width=0.5 + ) + ), + name='AT Skew Dist.', + hovertemplate='AT Skew: %{y}
Count: %{x}', + showlegend=False, + + ), + row=2, col=2 + ) + + fig.update_layout( + title='Skewness Distribution', + template='plotly_white', + bargap=0.05 + ) + + fig.update_xaxes(title_text="GC Skew (G-C)/(G+C)", range=[-1.05, 1.05], row=2, col=1) + fig.update_yaxes(title_text="AT Skew (A-T)/(A+T)", range=[-1.05, 1.05], row=2, col=1) + + fig.add_hline(y=0, line_dash="dash", line_color="rgba(0,0,0,0.2)", row=2, col=1) + fig.add_vline(x=0, line_dash="dash", line_color="rgba(0,0,0,0.2)", row=2, col=1) + config['toImageButtonOptions']['filename'] = "at_gc_skewness" + return fig.to_html(full_html=False, include_plotlyjs=False, config=config) + +def aa_group_distribution(group_dist: dict): + groups = list(group_dist.keys()) + values = [round(v * 100, 2) for v in group_dist.values()] + + fig = px.bar( + x=groups, + y=values, + title="Amino Acid Group Distribution Across All Sequences", + labels={"x": "Amino Acid Group", "y": "Proportion (%)"}, + color_discrete_sequence=['#0F65A0'], + ) + fig.update_traces( + marker_line_color="#616D78", + marker_line_width=1, + opacity=0.85, + ) + fig.update_layout( + template="plotly_white", + xaxis_title="Amino Acid Group", + yaxis_title="Proportion (%)", + bargap=0.3, + ) + config['toImageButtonOptions']['filename'] = "aa_group_distribution" + return fig.to_html(full_html=False, include_plotlyjs=False, config=config) \ No newline at end of file diff --git a/src/biological/sequence_data.py b/src/biological/sequence_data.py new file mode 100644 index 0000000..5ec4324 --- /dev/null +++ b/src/biological/sequence_data.py @@ -0,0 +1,463 @@ +import re +import tempfile +from collections import defaultdict, Counter +from itertools import chain +from pathlib import Path +from typing import Tuple, List, Dict +from urllib.error import URLError, HTTPError + +from Bio import motifs +from Bio.SeqUtils.ProtParam import ProteinAnalysis +import numpy as np +import pandas as pd +import peptides +import plotly.express as px +import ssl + +from analysis.outlier_detection import detect_outliers +from analysis.plot_utils import apply_standard_axes +from biological.plot_utils import length_distribution, gc_distribution, ambiguous_distribution, at_gc_skewness, \ + aa_group_distribution +from models.sequence import DNARNAColumns, PROTEINColumns, SequenceMetricSummary + +ssl._create_default_https_context = ssl._create_stdlib_context + + +# ToDo: Add Composition over all + + +def count_nmer(sequence, n) -> defaultdict: + if n > len(sequence): + return defaultdict(int) + seqs = np.frombuffer(sequence.encode('utf-8'), dtype='S1') + nmere = np.lib.stride_tricks.sliding_window_view(seqs, n) + unique, counts = np.unique(nmere, return_counts=True, axis=0) + result = defaultdict(int) + for nmer, count in zip(unique, counts): + nmer_str = b"".join(nmer).decode('utf-8') + result[nmer_str] += count + return result + + +def top_mere(seq, n=3, top=5) -> List[Tuple[str, int]] | None: + if not seq or len(seq) < n: + return None + counts = count_nmer(seq, n) + return sorted(counts.items(), key=lambda x: x[1], reverse=True)[:top] + + +def biological_data_top_entries(seqs: pd.Series, top_k: int = 20) -> Tuple[ + np.ndarray, np.ndarray, int, int, np.ndarray]: + seqs = seqs.str.upper() + vc = seqs.value_counts() + uniq_tmp, counts_tmp = vc.index.to_numpy(), vc.values + + top_k = min(top_k, len(uniq_tmp)) + top_idx = np.argsort(counts_tmp)[::-1][:top_k] + + uniques = uniq_tmp[top_idx].astype(str) + counts = counts_tmp[top_idx] + + lengths = np.array([len(s) for s in uniques]) + min_len, max_len = lengths.min(), lengths.max() + + return uniques, counts, min_len, max_len, lengths + + +def dna_rna_columns(seqs: pd.Series, k: int = 3, top_n: int = 20, top: int = 5, invalid: list=None) -> DNARNAColumns: + if invalid: + seqs = seqs[~seqs.isin(invalid)] + seqs.dropna(inplace=True) + #Over Top N + uniques, counts, min_len, max_len, lengths = biological_data_top_entries(seqs, top_n) + nucleotide_count = [dict(Counter(s)) for s in uniques] + k_mers = _kmer_check(k, top, uniques) + + #Over all sequences + all_overview = pd.DataFrame(seqs, columns=['sequence']) + all_overview['sequence'] = seqs.str.upper() + all_overview["lengths"] = all_overview['sequence'].str.len() + + A = all_overview['sequence'].str.count('A') + T = all_overview['sequence'].str.count('T') + G = all_overview['sequence'].str.count('G') + C = all_overview['sequence'].str.count('C') + all_overview['GC_skew'] = (G - C) / (G + C) + all_overview['AT_skew'] = (A - T) / (A + T) + + all_overview["GC_count"] = all_overview['sequence'].str.count('[GC]') + all_overview["gc_content"] = np.round(np.where(all_overview["lengths"] > 0, all_overview["GC_count"] / all_overview["lengths"] * 100, 0.0), 2) + + all_overview["N"] = all_overview['sequence'].str.count('N') + all_overview["ambiguous_base_ratio"] = np.round(np.where(all_overview["lengths"] > 0, all_overview["N"] / all_overview["lengths"] * 100, 0.0), 2) + + all_overview["codon_complete"] = all_overview["lengths"] % 3 + all_overview["codon_pct"] = np.where(all_overview["lengths"] > 0, (all_overview["lengths"] - all_overview["codon_complete"]) / all_overview["lengths"] * 100,0.0) + codon_completeness = SequenceMetricSummary( + min=round(float(np.nanmin(all_overview["codon_pct"] )), 2), + max=round(float(np.nanmax(all_overview["codon_pct"] )), 2), + mean=round(float(np.nanmean(all_overview["codon_pct"] )), 2), + ) + + gc_content = SequenceMetricSummary( + min=round(float(all_overview["gc_content"].min()), 2), + max=round(float(all_overview["gc_content"].max()), 2), + mean=round(float(all_overview["gc_content"].mean()), 2), + ) + + ambiguous_base_ratio = SequenceMetricSummary( + min=round(float(all_overview["ambiguous_base_ratio"].min()), 2), + max=round(float(all_overview["ambiguous_base_ratio"].max()), 2), + mean=round(float(all_overview["ambiguous_base_ratio"].mean()), 2), + ) + + length_stats = SequenceMetricSummary( + min=round(float(all_overview["lengths"].min()), 2), + max=round(float(all_overview["lengths"].max()), 2), + mean=round(float(all_overview["lengths"].mean()), 2), + ) + + gc_skew = SequenceMetricSummary( + min=round(float(all_overview["GC_skew"].min()), 4), + max=round(float(all_overview["GC_skew"].max()), 4), + mean=round(float(all_overview["GC_skew"].mean()), 4), + ) + + at_skew = SequenceMetricSummary( + min=round(float(all_overview["AT_skew"].min()), 4), + max=round(float(all_overview["AT_skew"].max()), 4), + mean=round(float(all_overview["AT_skew"].mean()), 4), + ) + + cpg_oe, tpa_oe = _dinucleotide_oe(all_overview) + + all_overview["entropy"] = all_overview["sequence"].apply(_normalized_shanon_entropy) + + low_complexity = SequenceMetricSummary( + min=round(float(all_overview["entropy"].min()), 2), + max=round(float(all_overview["entropy"].max()), 2), + mean=round(float(all_overview["entropy"].mean()), 2), + ) + + reverse_complement_ratio, reverse_complement_list = _reverse_complement_duplicates(all_overview['sequence']) + at_gc_plot = at_gc_skewness(all_overview) + gc_dist_plot = gc_distribution(all_overview) + + affected_sequences = (all_overview["N"] > 0).sum() + ambiguous_dist_plot = None + if affected_sequences > 0: + ambiguous_dist_plot = ambiguous_distribution(all_overview) + + length_dist_plot = None + if all_overview["lengths"].nunique() > 1: + length_dist_plot = length_distribution(all_overview) + + length_outliers = detect_outliers(all_overview["lengths"].to_numpy(dtype=np.double)) + + + if min_len == max_len: + plot = make_logo(uniques, 'color_classic', seq_type="dna") + else: + flat_kmers = chain.from_iterable(k_mers) + df_kmers = pd.DataFrame(flat_kmers, columns=['kmer', 'count']) + aggregated = df_kmers.groupby('kmer', as_index=False)['count'].sum() + plot = plot_overview(aggregated['kmer'].tolist(), aggregated['count'].tolist()) + + return DNARNAColumns( + sequence=uniques.tolist(), + count=counts.tolist(), + length=lengths.tolist(), + gc_content=gc_content, + ambiguous_base_ratio=ambiguous_base_ratio, + length_stats=length_stats, + length_outliers=length_outliers, + codon_completeness=codon_completeness, + gc_skew=gc_skew, + at_skew=at_skew, + cpg_observed_expected=cpg_oe, + tpa_observed_expected=tpa_oe, + low_complexity=low_complexity, + reverse_complement_ratio=reverse_complement_ratio, + reverse_complement_list=reverse_complement_list, + nucleotide_count=nucleotide_count, + k_mers=k_mers, + plot=plot, + gc_distribution=gc_dist_plot, + length_distribution=length_dist_plot, + ambiguous_distribution=ambiguous_dist_plot, + at_gc_skewness=at_gc_plot, + ) + +def _kmer_check(k: int, top: int, uniques: np.ndarray) -> list: + check_length = any(len(i) <= k for i in uniques) + if not check_length: + k_mers = [top_mere(seq, n=k, top=top) for seq in uniques] + else: + print("Warning: Sequence length is smaller than choosen k-Mer Size. Setting k-Mer Size to 3.") + k_mers = [top_mere(seq, n=3, top=top) for seq in uniques] + return k_mers + +def _reverse_complement_duplicates(all_seqs: pd.Series) -> Tuple[float,set]: + complement = str.maketrans('ACGT', 'TGCA') + seq_set = set(all_seqs) + canonical = {min(seq, seq.translate(complement)[::-1]) for seq in seq_set} + redundant_seqs = seq_set - canonical + redundant = len(seq_set) - len(canonical) + return (round(redundant / len(seq_set) * 100, 2) if len(seq_set) > 0 else 0.0), redundant_seqs + +def _normalized_shanon_entropy(seq: str) -> np.float64: + n = len(seq) + if n == 0: + return np.float64(0.0) + counts = np.fromiter(Counter(seq).values(), dtype=np.float64) + max_entropy = np.log2(min(n, len(counts))) # ToDo or 4? + if max_entropy == 0: + return np.float64(0.0) + p = counts / n + h = -np.sum(p * np.log2(p)) + return np.float64(round(float(h / max_entropy * 100), 2)) + +def _dinucleotide_oe(df: pd.DataFrame) -> Tuple[SequenceMetricSummary, SequenceMetricSummary]: + seqs = df["sequence"] + lengths = df["lengths"] + + a_freq = seqs.str.count('A') / lengths + t_freq = seqs.str.count('T') / lengths + g_freq = seqs.str.count('G') / lengths + c_freq = seqs.str.count('C') / lengths + + tpa_obs = seqs.str.count('TA') / (lengths - 1) + tpa_exp = t_freq * a_freq + tpa_oe = np.where(tpa_exp > 0 , tpa_obs/tpa_exp, 0.0) + + cpg_obs = seqs.str.count('CG') / (lengths - 1) + cpg_exp = g_freq * c_freq + cpg_oe = np.where(cpg_exp > 0 , cpg_obs/cpg_exp, 0.0) + + #ToDo check if wie only use values above 0.0 + return ( + SequenceMetricSummary( + min=round(float(np.nanmin(cpg_oe)), 4), + max=round(float(np.nanmax(cpg_oe)), 4), + mean=round(float(np.nanmean(cpg_oe)), 4), + ), + SequenceMetricSummary( + min=round(float(np.nanmin(tpa_oe)), 4), + max=round(float(np.nanmax(tpa_oe)), 4), + mean=round(float(np.nanmean(tpa_oe)), 4), + ), + ) + +def protein_descriptors(peptide: str) -> Dict[str, str | float | dict[str, float]]: + descriptors: Dict[str, str | float | dict[str, float]] = {} + p: peptides.Peptide = peptides.Peptide(peptide) + descriptors["seq"] = peptide + descriptors["freq"] = p.frequencies() + + try: + descriptors["aidx"] = p.aliphatic_index() + except ZeroDivisionError: + descriptors["aidx"] = 0.0 + descriptors["boman"] = p.boman() + descriptors["charge"] = p.charge() + descriptors["hp"] = p.hydrophobicity() + descriptors["iep"] = p.isoelectric_point() + descriptors["iidx"] = p.instability_index() + descriptors["mol"] = p.molecular_weight() + descriptors["aroma"] = sum([peptide.count(aa) for aa in ('F', 'W', 'Y')]) / len(peptide) + return descriptors + +#ToDo Add Desclaimer; Removing invalid from list +def protein_columns(seqs: pd.Series, k: int = 3, top_n: int = 20, top: int = 5, invalid: list=None) -> PROTEINColumns: + # Top-N + if invalid: + seqs = seqs[~seqs.isin(invalid)] + uniques, counts, min_len, max_len, lengths = biological_data_top_entries(seqs, top_n) + aa_composition = [dict(Counter(seq)) for seq in uniques] + descriptors = [protein_descriptors(seq) for seq in uniques] + k_mers = _kmer_check(k, top, uniques) + + # Column-wide + all_overview = pd.DataFrame(seqs, columns=['sequence']) + all_overview['sequence'] = seqs.str.upper() + all_overview['lengths'] = all_overview['sequence'].str.len() + + all_overview['ambiguous'] = all_overview['sequence'].str.count('[XJU]') + all_overview['ambiguous_ratio'] = np.round(np.where(all_overview['lengths'] > 0, all_overview['ambiguous'] / all_overview['lengths'] * 100, 0.0), 2) + + all_overview['stop_codon'] = all_overview['sequence'].str.count(r'\*') + all_overview['entropy'] = all_overview['sequence'].apply(_normalized_shanon_entropy) + + ambiguous_residue_ratio = SequenceMetricSummary( + min=round(float(all_overview['ambiguous_ratio'].min()), 2), + max=round(float(all_overview['ambiguous_ratio'].max()), 2), + mean=round(float(all_overview['ambiguous_ratio'].mean()), 2), + ) + + length_stats = SequenceMetricSummary( + min=round(float(all_overview['lengths'].min()), 2), + max=round(float(all_overview['lengths'].max()), 2), + mean=round(float(all_overview['lengths'].mean()), 2), + ) + + length_outliers = detect_outliers(all_overview['lengths'].to_numpy(dtype=np.double)) + + stop_codon_ratio = round((all_overview['stop_codon'] > 0).sum() / len(all_overview) * 100, 2) + + low_complexity = SequenceMetricSummary( + min=round(float(all_overview['entropy'].min()), 2), + max=round(float(all_overview['entropy'].max()), 2), + mean=round(float(all_overview['entropy'].mean()), 2), + ) + + all_overview['gravy'] = all_overview['sequence'].apply(_gravy) + all_overview['cysteine'] = all_overview['sequence'].str.count('C') + all_overview['disorder'] = all_overview['sequence'].str.count(f"[{''.join(DISORDER_AA)}]") + all_overview['disorder_ratio'] = np.round( + np.where(all_overview['lengths'] > 0, all_overview['disorder'] / all_overview['lengths'] * 100, 0.0), 2 + ) + + gravy = SequenceMetricSummary( + min=round(float(all_overview['gravy'].min()), 4), + max=round(float(all_overview['gravy'].max()), 4), + mean=round(float(all_overview['gravy'].mean()), 4), + ) + + cysteine_count = SequenceMetricSummary( + min=round(float(all_overview['cysteine'].min()), 2), + max=round(float(all_overview['cysteine'].max()), 2), + mean=round(float(all_overview['cysteine'].mean()), 2), + ) + + disorder_propensity = SequenceMetricSummary( + min=round(float(all_overview['disorder_ratio'].min()), 2), + max=round(float(all_overview['disorder_ratio'].max()), 2), + mean=round(float(all_overview['disorder_ratio'].mean()), 2), + ) + + aa_group_dist = _aa_group_distribution(all_overview['sequence']) + aa_group_plot = aa_group_distribution(aa_group_dist) + + length_dist_plot = None + if all_overview['lengths'].nunique() > 1: + length_dist_plot = length_distribution(all_overview, unit="residues") + + affected = (all_overview['ambiguous'] > 0).sum() + ambiguous_dist_plot = None + if affected > 0: + ambiguous_dist_plot = ambiguous_distribution(all_overview, col='ambiguous', label='X/J/U') + + plot = None + if min_len == max_len: + plot = make_logo(uniques, "chemistry", seq_type="protein") + if plot is None: + flat_kmers = chain.from_iterable(k_mers) + df_kmers = pd.DataFrame(flat_kmers, columns=['kmer', 'count']) + aggregated = df_kmers.groupby('kmer', as_index=False)['count'].sum() + plot = plot_overview(aggregated['kmer'].tolist(), aggregated['count'].tolist()) + + return PROTEINColumns( + sequence=uniques.tolist(), + count=counts.tolist(), + length=lengths.tolist(), + composition=aa_composition, + frequency=[descriptor['freq'] for descriptor in descriptors], + hydrophobicity=[descriptor['hp'] for descriptor in descriptors], + charge=[descriptor['charge'] for descriptor in descriptors], + molecular_weight=[descriptor['mol'] for descriptor in descriptors], + isoelectric_point=[descriptor['iep'] for descriptor in descriptors], + aliphatic_index=[descriptor['aidx'] for descriptor in descriptors], + boman=[descriptor['boman'] for descriptor in descriptors], + aromaticity=[descriptor['aroma'] for descriptor in descriptors], + instability=[descriptor['iidx'] for descriptor in descriptors], + ambiguous_residue_ratio=ambiguous_residue_ratio, + length_stats=length_stats, + length_outliers=length_outliers, + stop_codon_ratio=stop_codon_ratio, + low_complexity=low_complexity, + k_mers=k_mers, + plot=plot, + length_distribution=length_dist_plot, + ambiguous_distribution=ambiguous_dist_plot, + gravy=gravy, + cysteine_count=cysteine_count, + disorder_propensity=disorder_propensity, + aa_group_distribution=aa_group_dist, + aa_group_plot=aa_group_plot + ) + +AA_GROUPS = { + "Unpolar": set("GAVLIMP"), + "Aromatic": set("FWY"), + "Polar": set("STCNQ"), + "Positive": set("KRH"), + "Negative": set("DE"), +} + +DISORDER_AA = set("PESQK") +KYTE_DOOLITTLE = {'A': 1.8, 'R': -4.5, 'N': -3.5, 'D': -3.5, 'C': 2.5,'Q': -3.5, 'E': -3.5, 'G': -0.4, 'H': -3.2, + 'I': 4.5, 'L': 3.8, 'K': -3.9, 'M': 1.9, 'F': 2.8, 'P': -1.6, 'S': -0.8, 'T': -0.7, 'W': -0.9, 'Y': -1.3, 'V': 4.2,} + +def _gravy(seq: str) -> float: + if len(seq) == 0: + return 0.0 + try: + return ProteinAnalysis(seq).gravy() + except (KeyError, ValueError): + # Fallback for sequences with ambiguous residues (X, J, U, *) + return sum(KYTE_DOOLITTLE.get(aa, 0.0) for aa in seq) / len(seq) + +def _aa_group_distribution(all_seqs: pd.Series) -> dict: + concat = all_seqs.str.cat() + total = len(concat) + if total == 0: + return {g: 0.0 for g in AA_GROUPS} + counts = Counter(concat) + return { + group: round(sum(counts.get(aa, 0) for aa in members) / total, 4) + for group, members in AA_GROUPS.items() + } + +def make_logo(seqs, color, seq_type): + if seq_type == "protein": + m = motifs.create(seqs, alphabet="ACDEFGHIKLMNPQRSTVWY") + else: + m = motifs.create(seqs, alphabet="ACGT") + with tempfile.NamedTemporaryFile(suffix='.svg', delete=False) as tmp_file: + tmp_path = tmp_file.name + + try: + m.weblogo(tmp_path, format="svg", sequence_type=seq_type, color_scheme=color, logo_font="Calibri", logo_margin=3, + fontsize=12, scale_width=True) + + with open(tmp_path, 'r', encoding='utf-8') as svg_file: + raw = svg_file.read() + + # WebLogo 3 wraps the SVG in a full HTML document — extract only the SVG fragment + # to avoid injecting Bootstrap 3, jQuery 1.x and other legacy dependencies into the report + svg_match = re.search(r'()', raw, re.IGNORECASE) + if svg_match: + svg_content = svg_match.group(1) + svg_content = svg_content.replace( + ' 0.3: + return "None", [] + + if all(len(x) > 2 for x in values): + match, invalid = fast_check_sequence(values, Sequence.DNA.value, threshold) + if match: + return "dna", _get_invalid(values, invalid) + match, invalid = fast_check_sequence(values, Sequence.RNA.value, threshold) + if match: + return "rna", _get_invalid(values, invalid) + match, invalid = fast_check_sequence(values, Sequence.PROTEIN.value, threshold) + if match: + if not invalid or char_entropy(values, PROTEIN_ALPHABET) >= ENTROPY_THRESHOLDS["protein"]: + return "protein", _get_invalid(values, invalid) + return "None", [] + + +def _get_invalid(values, invalid_indices): + if not invalid_indices: + return [] + return [values[i] for i in invalid_indices] \ No newline at end of file diff --git a/src/biological/taxonomy.py b/src/biological/taxonomy.py new file mode 100644 index 0000000..7388bd6 --- /dev/null +++ b/src/biological/taxonomy.py @@ -0,0 +1,160 @@ +from dataclasses import dataclass +from typing import Tuple + +import pandas as pd + +@dataclass +class TaxonomyFlags: + name: str + is_taxonomy: bool + taxid: set | str | None + taxonomy: list | str | None + rank_distribution: dict | None = None + is_mixed: bool = False + invalid_names: list | None = None + outdated_names : dict | None = None + + +def taxonomy_flags(df, col, valid_names, valid_tax_ids, name_to_rank, taxid_to_rank, + name_to_scientific) -> TaxonomyFlags: + if df[col].dtype in ['int64', 'float64'] or pd.api.types.is_numeric_dtype(df[col]): + taxid_result = is_taxid(df[col], valid_tax_ids) + if taxid_result is not None: + distribution, is_mixed = taxid_rank_distribution(df[col], taxid_to_rank) + return TaxonomyFlags( + name=col, + is_taxonomy=True, + taxid=taxid_result, + taxonomy=None, + rank_distribution=distribution, + is_mixed=is_mixed, + invalid_names=None, + ) + else: + taxonomy_result = is_taxonomy(df[col], valid_names, name_to_rank, name_to_scientific) + if taxonomy_result is not None: + return TaxonomyFlags( + name=col, + is_taxonomy=True, + taxid=None, + taxonomy=taxonomy_result["invalid_names"] or "Valid", + rank_distribution=taxonomy_result["rank_distribution"], + is_mixed=taxonomy_result["is_mixed"], + invalid_names=taxonomy_result["invalid_names"], + outdated_names=taxonomy_result["outdated"], + ) + return TaxonomyFlags( + name=col, + is_taxonomy=False, + taxid=None, + taxonomy=None, + rank_distribution=None, + is_mixed=False, + invalid_names=None + ) + + +def build_lookups(vocab: pd.DataFrame) -> Tuple[set, set, dict, dict, dict]: + valid_names = set(vocab['name_txt']) + valid_tax_ids = set(vocab['tax_id']) + name_to_rank = dict(zip(vocab['name_txt'], vocab['rank'])) + name_to_scientific = dict(zip(vocab['name_txt'], vocab['scientific_name'])) + + sci = vocab[vocab['name_class'] == 'scientific name'].drop_duplicates('tax_id') + taxid_to_rank = dict(zip(sci['tax_id'], sci['rank'])) + + return valid_names, valid_tax_ids, name_to_rank, taxid_to_rank, name_to_scientific + + +def is_taxid(col: pd.Series, valid_tax_ids: set, threshold: float = 0.9) -> set | str | None: + excluded_cols = ["length", "start", "end"] + + if col.name and str(col.name).lower() in excluded_cols: + return None + + tmp_series = pd.to_numeric(col, errors='coerce') + is_numeric_candidate = tmp_series.notna().sum() / len(col) > threshold + + if is_numeric_candidate: + is_valid = tmp_series.isin(valid_tax_ids) + validity_rate = is_valid.sum() / len(col) + + if validity_rate > threshold: + invalid_mask = ~is_valid & tmp_series.notna() + invalid_ids = set(col.loc[invalid_mask].tolist()) + + if invalid_ids: + return invalid_ids + else: + return "all tax IDs valid" + + return None + + +def is_taxonomy(col: pd.Series, valid_names: set, name_to_rank: dict, name_to_scientific: dict, threshold: float = 0.8) -> dict | None: + col_obj = col.astype(object) + uniques = pd.unique(col_obj.dropna()) + valid_unique = {u for u in uniques if u in valid_names} + + is_valid = col_obj.isin(valid_unique) + validity_rate = is_valid.sum() / len(col) + + cleaned_names = col_obj + if validity_rate < threshold: + cleaned_names = col_obj.astype(str).str.extract(r'^([^(]+)')[0].str.strip() + cleaned_uniques = pd.unique(cleaned_names.dropna()) + valid_cleaned_unique = {u for u in cleaned_uniques if u in valid_names} + is_valid_cleaned = cleaned_names.isin(valid_cleaned_unique) + validity_rate_cleaned = is_valid_cleaned.sum() / len(col) + + if validity_rate_cleaned > validity_rate: + is_valid = is_valid_cleaned + validity_rate = validity_rate_cleaned + + if validity_rate > threshold: + distribution, is_mixed, invalid_names = rank_distribution(cleaned_names, name_to_rank) + outdated = find_outdated_names(cleaned_names, valid_names, name_to_scientific) + return{ + "valid":True, + "rank_distribution":distribution, + "is_mixed":is_mixed, + "invalid_names":invalid_names if invalid_names else None, + "outdated":outdated if outdated else None, + } + return None + + +def taxid_rank_distribution(col: pd.Series, taxid_to_rank: dict, threshold: float = 0.05) -> Tuple[dict, bool]: + numeric = pd.to_numeric(col, errors='coerce') + ranks = numeric.map(taxid_to_rank) + + rank_counts = ranks.value_counts(normalize=True) + distribution = {rank: round(float(freq), 4) for rank, freq in rank_counts.items()} + + is_mixed = len([f for f in distribution.values() if f >= threshold]) > 1 + return distribution, is_mixed + +def rank_distribution(col: pd.Series, name_to_rank: dict, threshold: float = 0.05) -> Tuple[dict, bool, list]: + cleaned = col.astype(str).str.strip() + ranks = cleaned.map(name_to_rank) + + invalid_names = cleaned[ranks.isna() & col.notnull()].unique().tolist() + + rank_counts = ranks.value_counts(normalize=True) + distribution = {rank: round(float(freq), 4) for rank, freq in rank_counts.items()} + + is_mixed = len([f for f in distribution.values() if f >= threshold]) > 1 + return distribution, is_mixed, sorted(invalid_names) + + +def find_outdated_names(col: pd.Series, valid_names: set, name_to_scientific: dict) -> dict: + cleaned = col.astype(str).str.strip() + uniques = pd.unique(cleaned.dropna()) + valid_used = [u for u in uniques if u in valid_names] + + outdated = {} + for name in valid_used: + current = name_to_scientific.get(name) + if current is not None and name != current: + outdated[name] = current + return outdated \ No newline at end of file diff --git a/src/cli/__init__.py b/src/cli/__init__.py new file mode 100644 index 0000000..247b1cb --- /dev/null +++ b/src/cli/__init__.py @@ -0,0 +1 @@ +from .app import cli diff --git a/src/cli/app.py b/src/cli/app.py new file mode 100644 index 0000000..aa22eb8 --- /dev/null +++ b/src/cli/app.py @@ -0,0 +1,160 @@ +#!/usr/bin/env python3 + +import time +from pathlib import Path + +import click +import pandas as pd +from termcolor import colored + +from analysis.categorical_analysis import categorical_columns +from analysis.multivariate_analysis import multivariate_analysis +from analysis.numeric_analysis import numeric_columns +from analysis.overview import overview, column_overview +from biological.functional_annotation import annotation_flags +from biological.measurement_data import measurement_columns +from biological.sequence_data import dna_rna_columns, protein_columns +from biological.taxonomy import taxonomy_flags, build_lookups +from cli.report_json import write_result_json +from cli.report_writer import write_report +from cli.utils import _fmt_duration, print_step, info +from data_utils.file_reader import read_file, parse_parquet +from data_utils.remote_data import get_tax_ids +from quality_assessment.quality_assessment import quality_assessment, print_quality_report + +CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) + + +@click.group(context_settings=CONTEXT_SETTINGS, invoke_without_command=True) +@click.option("-i", "--input", type=click.Path(exists=True, resolve_path=True), required=True, + help="Input file as .tsv, .csv or .json") +@click.option('-t', '--tax', is_flag=True, help='Enable taxonomy analysis') +@click.option('-f', '--func', type=click.Choice(['cog', 'go']), + help='Enable functional annotation analysis. Choose between cog or go') +@click.option('-tc', '--target_column', type=str, help='Target column for Analysis') +@click.option('-k', '--kmer', type=int, default=3, help="K-mer Size for sequence analysis") +@click.option('-n', '--top_n', type=int, default=20, help="Top N entries analysis") +def cli(input: str, tax: bool = False, func: str = None, + target_column: str = None, kmer: int = None, top_n: int = None): + run_start = time.perf_counter() + input_path = Path(input) + + print(f"\n{'=' * 64}") + print(colored(f" BioProfileKit — {input_path.name}", "magenta", attrs=["bold"])) + print(f"{'=' * 64}") + + done = print_step(f"Reading {input_path.name}") + df = read_file(input_path) if input_path.suffix != '.parquet' else parse_parquet(input_path) + general = overview(df, input_path.name) + done(f"{df.shape[0]:,} rows × {df.shape[1]} columns") + + dups = df[df.duplicated(keep=False)].reset_index() + duplicates_table = dups.to_html( + classes="table table-hover table-responsive nowrap", + border="0", table_id="dup_table", index=False, + ) + + done = print_step(f"Column overview ({len(df.columns)} columns)") + column_overviews = [column_overview(df, col) for col in df.columns] + done() + + tax_df = get_tax_ids() if tax else None + valid_names, valid_tax_ids, name_to_rank, taxid_to_rank, name_to_scientific = None, None, None, None, None + if tax and tax_df is not None: + done = print_step("Building taxonomy lookups") + valid_names, valid_tax_ids, name_to_rank, taxid_to_rank, name_to_scientific = build_lookups(tax_df) + done(f"{len(valid_names):,} scientific names") + + done = print_step("Per-column analysis") + seq_count = 0 + for col_ov in column_overviews: + if tax and tax_df is not None: + col_ov.taxonomy = taxonomy_flags(df, col_ov.name, valid_names, valid_tax_ids, name_to_rank, taxid_to_rank, name_to_scientific) + + if func and col_ov.taxonomy is not None: + col_ov.annotation = [annotation_flags(df, col_ov.name, func)] + + if hasattr(col_ov, "top_10") and isinstance(col_ov.top_10, pd.Series): + col_ov.top_10_items = list(col_ov.top_10.items()) + + if col_ov.sequence == 'dna': + info(f"DNA sequences: {col_ov.name}") + col_ov.dna_rna_data = dna_rna_columns(df[col_ov.name], k=kmer, top_n=top_n, invalid=col_ov.invalid_seqs) + seq_count += 1 + elif col_ov.sequence == 'rna': + info(f"RNA sequences: {col_ov.name}") + col_ov.dna_rna_data = dna_rna_columns(df[col_ov.name], k=kmer, top_n=top_n, invalid=col_ov.invalid_seqs) + seq_count += 1 + elif col_ov.sequence == 'protein': + info(f"protein sequences: {col_ov.name}") + col_ov.protein_data = protein_columns(df[col_ov.name], k=kmer, top_n=top_n, invalid=col_ov.invalid_seqs) + seq_count += 1 + else: + col_ov.dna_rna_data = None + col_ov.protein_data = None + + if col_ov.sequence == 'None': + measurement_data = measurement_columns(df[col_ov.name], col_ov.name, col_ov.type) + if measurement_data: + info(f"lab measurements: {col_ov.name}") + col_ov.measurement_data = measurement_data if measurement_data else None + else: + col_ov.measurement_data = None + done(f"{seq_count} sequence column(s)") + + empty_cols = [col for col in df.columns if df[col].isnull().all()] + sequence_cols = {col.name for col in column_overviews if col.sequence in ('dna', 'rna', 'protein')} + + numeric_cols = [col for col in df.select_dtypes(include="number").columns + if col not in empty_cols and col not in sequence_cols] + done = print_step(f"Numeric analysis ({len(numeric_cols)} columns)") + numeric_overviews = [numeric_columns(df, col) for col in numeric_cols] + done() + + cat_columns = [col for col in df.select_dtypes(include=['str', 'object', 'bool']).columns + if col not in empty_cols and col not in sequence_cols] + done = print_step(f"Categorical analysis ({len(cat_columns)} columns)") + categorical_overviews = [categorical_columns(df, col) for col in cat_columns] + done() + + done = print_step("Multivariate analysis") + plots = multivariate_analysis(df, target_column) + done() + + done = print_step("Quality assessment") + quality = quality_assessment(general, column_overviews, numeric_overviews, + categorical_overviews, plots) + done() + print_quality_report(quality) + + general.n_number = len(numeric_cols) + general.n_categorical = len(cat_columns) + general.n_empty = len(empty_cols) + general.n_dna = sum(1 for col in column_overviews if col.sequence == 'dna') + general.n_rna = sum(1 for col in column_overviews if col.sequence == 'rna') + general.n_protein = sum(1 for col in column_overviews if col.sequence == 'protein') + general.n_taxonomy = sum(1 for col in column_overviews + if getattr(col, 'taxonomy', None) and col.taxonomy.is_taxonomy) + general.n_unit = sum(1 for col in column_overviews if col.measurement_data is not None) + general.n_functional = sum(1 for col in column_overviews if getattr(col, 'annotation', None)) + general.n_taxonomy_candidates = sum( + 1 for col in column_overviews + if getattr(col, 'taxonomy_candidate', False) + ) + + output_path = Path(input_path.stem + "_renders") + done = print_step("Writing report") + ctx = click.get_current_context() + parameters = {k: v for k, v in ctx.params.items() if k != "input"} + parameters["input_file"] = Path(input).name + + write_report(output_path, general, plots, duplicates_table, + column_overviews, numeric_overviews, categorical_overviews, top_n, quality) + write_result_json(output_path / "results.json", general, column_overviews, + numeric_overviews, categorical_overviews, plots, quality, empty_cols, parameters) + done(f"→ {output_path}/") + + total = _fmt_duration(time.perf_counter() - run_start) + print(f"\n{'=' * 64}") + print(colored(f" Done in {total}", "green", attrs=["bold"])) + print(f"{'=' * 64}\n") diff --git a/src/cli/report_json.py b/src/cli/report_json.py new file mode 100644 index 0000000..84a4756 --- /dev/null +++ b/src/cli/report_json.py @@ -0,0 +1,342 @@ +import json +import math +from datetime import datetime, timezone + +import numpy as np + +SCHEMA_VERSION = "1.0" +TOOL_VERSION = "0.1.0" +CORR_THRESHOLD = 0.9 #ToDo check consistencies + +def _num(v): + if v is None: + return None + if isinstance(v, np.integer): + return int(v) + if isinstance(v, (np.floating, float)): + f = float(v) + return f if math.isfinite(f) else None + if isinstance(v, (int, bool, str)): + return v + if isinstance(v, np.ndarray): + return [_num(x) for x in v.tolist()] + return v + + +def _present(v): + if v is None: + return False + if isinstance(v, np.ndarray): + return v.size > 0 + if hasattr(v, "empty"): # pandas Series / DataFrame + return not v.empty + if isinstance(v, (list, tuple, set, dict, str)): + return len(v) > 0 + return bool(v) + + +def _listify(v): + if not _present(v): + return None + if isinstance(v, np.ndarray): + return [_num(x) for x in v.tolist()] + if hasattr(v, "tolist"): # pandas Series + return [_num(x) for x in v.tolist()] + if isinstance(v, (list, tuple, set)): + return [_num(x) for x in v] + return _num(v) + + +def _summary(s): + if s is None: + return None + return {"min": _num(s.min), "max": _num(s.max), "mean": _num(s.mean)} + + +def _outliers(o): + if o is None: + return None + return { + "lower_bound": _num(o.lower_bound), + "upper_bound": _num(o.upper_bound), + "medcouple": _num(o.medcouple), + "n_lower_iqr": _num(o.n_lower_iqr), + "n_upper_iqr": _num(o.n_upper_iqr), + "n_lower_mzscore": _num(o.n_lower_mzscore), + "n_upper_mzscore": _num(o.n_upper_mzscore), + "n_lower_zscore": _num(o.n_lower_zscore), + "n_upper_zscore": _num(o.n_upper_zscore), + } + +def _numeric_block(n): + if n is None: + return None + return { + "min": _num(n.min), "max": _num(n.max), "mean": _num(n.mean), + "median": _num(n.median), "mode": _num(n.mode), "std": _num(n.std), + "sum": _num(n.sum), "kurtosis": _num(n.kurtosis), "skewness": _num(n.skewness), + "coefficient_of_variation": _num(n.coefficient_of_variation), "mad": _num(n.mad), + "quantiles": _num(n.quantiles), + "infinity": _num(n.infinity), "negative_count": _num(n.negative_count), + "zero_count": _num(n.zero_count), + "outliers": _outliers(n.outliers), + } + + +def _categorical_block(c): + if c is None: + return None + return { + "unique_categories": _num(c.unique_categories), "mode": c.mode, + "entropy": _num(c.entropy), "gini": _num(c.gini), + "simpson_diversity": _num(c.simpson_diversity), + "max_category_length": _num(c.max_category_length), + "min_category_length": _num(c.min_category_length), + "cardinality_ratio": _num(c.cardinality_ratio), + "rare_categories": _num(c.rare_categories), + "top_5_coverage": _num(c.top_5_coverage), + "top_1_coverage": _num(c.top_1_coverage), + "cib_ratio": _num(c.cib_ratio), + "effective_cardinality": _num(c.effective_cardinality), + } + + +def _dna_block(d): + return { + "gc_content": _summary(d.gc_content), + "ambiguous_base_ratio": _summary(d.ambiguous_base_ratio), + "length_stats": _summary(d.length_stats), + "length_outliers": _outliers(d.length_outliers), + "codon_completeness": _summary(d.codon_completeness), + "gc_skew": _summary(d.gc_skew), "at_skew": _summary(d.at_skew), + "cpg_observed_expected": _summary(d.cpg_observed_expected), + "tpa_observed_expected": _summary(d.tpa_observed_expected), + "low_complexity": _summary(d.low_complexity), + "reverse_complement_ratio": _num(d.reverse_complement_ratio), + } + + +def _protein_block(p): + return { + "ambiguous_residue_ratio": _summary(p.ambiguous_residue_ratio), + "length_stats": _summary(p.length_stats), + "length_outliers": _outliers(p.length_outliers), + "stop_codon_ratio": _num(p.stop_codon_ratio), + "low_complexity": _summary(p.low_complexity), + "gravy": _summary(p.gravy), + "cysteine_count": _summary(p.cysteine_count), + "disorder_propensity": _summary(p.disorder_propensity), + "aa_group_distribution": {k: _num(v) for k, v in p.aa_group_distribution.items()} + if p.aa_group_distribution else None, + } + + +def _taxonomy_block(tax): + if tax is None or not getattr(tax, "is_taxonomy", False): + return None + block = {"is_taxonomy": True, "is_mixed": bool(getattr(tax, "is_mixed", False))} + if getattr(tax, "rank_distribution", None): + block["rank_distribution"] = {k: _num(v) for k, v in tax.rank_distribution.items()} + if getattr(tax, "outdated_names", None): + block["name_corrections"] = tax.outdated_names + if getattr(tax, "invalid_names", None): + block["invalid_values"] = tax.invalid_names + return block + + +def _measurement_block(m): + if m is None: + return None + return { + "units": m.units, + "unit_counts": {str(k): _num(v) for k, v in m.unit_counts.items()}, + "with_measurement": bool(m.with_measurement), + } + + +def _correlated_partners(col_name, top_associations, threshold): + if not top_associations: + return [] + best = {} + for p in top_associations: + if p["value"] < threshold: + continue + if p["var1"] == col_name: + partner = p["var2"] + elif p["var2"] == col_name: + partner = p["var1"] + else: + continue + if partner not in best or p["value"] > best[partner]["value"]: + best[partner] = {"column": partner, "value": _num(p["value"]), "method": p["method"]} + return sorted(best.values(), key=lambda d: d["value"], reverse=True) + + +def _role(col_ov, is_empty): + if is_empty: + return "empty" + if col_ov.sequence == "dna": + return "dna_rna_sequence" + if col_ov.sequence == "protein": + return "protein_sequence" + tax = getattr(col_ov, "taxonomy", None) + if tax is not None and getattr(tax, "is_taxonomy", False): + return "taxonomy" + if getattr(col_ov, "measurement_data", None) is not None: + return "measurement" + return None + + +def _column_entry(col_ov, num, cat, multivariate, is_empty): + name = col_ov.name + entry = {"name": name, "dtype": col_ov.type} + + role = _role(col_ov, is_empty) + if role is not None: + entry["role"] = role + + if is_empty: + entry["empty"] = True + return entry + + entry["general"] = { + "count": _num(col_ov.number), + "n_unique": _num(col_ov.unique), + "missing": _num(col_ov.missing), + "missing_rate": _num(round((col_ov.missing_per or 0.0) / 100.0, 4)), + "density": _num(col_ov.density), + "constant": bool(col_ov.constant) if col_ov.constant is not None else None, + "mixed_types": _listify(col_ov.mixed_types), + "suspect_values": _listify(col_ov.suspect_values), + "monotonicity": bool(col_ov.monotonicity) if col_ov.monotonicity is not None else None, + "cardinality_dimension_ratio": _num(col_ov.cardinality_dimension_ratio), + "invalid_seqs": _listify(col_ov.invalid_seqs), + } + + nb = _numeric_block(num) + if nb is not None: + entry["numeric"] = nb + cb = _categorical_block(cat) + if cb is not None: + entry["categorical"] = cb + + # sequence block + pdata = getattr(col_ov, "protein_data", None) + ddata = getattr(col_ov, "dna_rna_data", None) + if pdata is not None: + entry["sequence"] = _protein_block(pdata) + elif ddata is not None: + entry["sequence"] = _dna_block(ddata) + + tb = _taxonomy_block(getattr(col_ov, "taxonomy", None)) + if tb is not None: + entry["taxonomy"] = tb + + mb = _measurement_block(getattr(col_ov, "measurement_data", None)) + if mb is not None: + entry["measurement"] = mb + + # multivariate relations attached per column + target = {} + ftc = getattr(multivariate, "feature_target_correlation", None) + if ftc and name in ftc: + target["association"] = _num(ftc[name]["value"]) + target["method"] = ftc[name]["method"] + mi = getattr(multivariate, "mutual_information", None) + if mi and name in mi: + target["mutual_information"] = _num(mi[name]["value"]) + if target: + entry["target_relation"] = target + + partners = _correlated_partners(name, getattr(multivariate, "top_associations", None), CORR_THRESHOLD) + if partners: + entry["correlated_with"] = partners + + return entry + +def _quality_block(quality): + if quality is None: + return None + return { + "overall": quality.overall, + "passed": quality.passed, "warnings": quality.warnings, + "failed": quality.failed, "total": quality.total, + "categories": [ + { + "name": cat.name, "status": cat.status, + "checks": [ + {"name": c.name, "status": c.status, "message": c.message} + for c in cat.checks + ], + } + for cat in quality.categories + ], + } + + +def build_result(general, column_overviews, numeric_overviews, + categorical_overviews, multivariate, quality, empty_cols, + parameters=None) -> dict: + num_by = {n.name: n for n in numeric_overviews if n is not None} + cat_by = {c.name: c for c in categorical_overviews if c is not None} + empty = set(empty_cols) + + columns = [ + _column_entry(col, num_by.get(col.name), cat_by.get(col.name), + multivariate, col.name in empty) + for col in column_overviews + ] + + return { + "meta": { + "schema_version": SCHEMA_VERSION, + "tool_version": TOOL_VERSION, + "generated_at": datetime.now(timezone.utc).isoformat(), + "parameters": parameters or {}, + "correlation_threshold": CORR_THRESHOLD, + }, + "summary": { + "filename": general.filename, + "rows": _num(general.rows), "cols": _num(general.cols), + "nulls": _num(general.nulls), "nulls_percentage": _num(general.nulls_percentage), + "empty_rows": _num(general.empty_rows), + "duplicate_rows": _num(general.dup_row), "duplicate_cols": _num(general.dup_col), + "rows_per_feature": _num(general.ratio), + "memory_bytes": _num(general.memory), + "column_types": { + "number": _num(general.n_number), "dna": _num(general.n_dna), + "rna": _num(general.n_rna), "protein": _num(general.n_protein), + "taxonomy": _num(general.n_taxonomy), "unit": _num(general.n_unit), + "functional": _num(general.n_functional), + "categorical": _num(general.n_categorical), "empty": _num(general.n_empty), + }, + }, + "target": getattr(multivariate, "target_name", None), + "columns": columns, + } + +def _json_default(o): + """Last-resort converter for anything that slipped through un-cleaned, so + json.dump never crashes on a Series/ndarray/numpy scalar.""" + if isinstance(o, np.integer): + return int(o) + if isinstance(o, np.floating): + f = float(o) + return f if math.isfinite(f) else None + if isinstance(o, np.ndarray): + return o.tolist() + if hasattr(o, "tolist"): + return o.tolist() + if hasattr(o, "item"): + return o.item() + return str(o) + +def write_result_json(path, general, column_overviews, numeric_overviews, + categorical_overviews, multivariate, quality, empty_cols, + parameters=None): + result = build_result(general, column_overviews, numeric_overviews, + categorical_overviews, multivariate, quality, empty_cols, + parameters) + with open(path, "w", encoding="utf-8") as f: + json.dump(result, f, indent=2, ensure_ascii=False, default=_json_default) + return path \ No newline at end of file diff --git a/src/cli/report_writer.py b/src/cli/report_writer.py new file mode 100644 index 0000000..820d89f --- /dev/null +++ b/src/cli/report_writer.py @@ -0,0 +1,48 @@ +import shutil +from pathlib import Path + +from jinja2 import Environment, FileSystemLoader +from importlib_resources import files + +TEMPLATE_DIR = files("templates").joinpath() +STATIC_DIR = files("static").joinpath() +env = Environment(loader=FileSystemLoader(str(TEMPLATE_DIR)), autoescape=True) + +PAGE_MAP = { + "overview": "numeric_data.html", + "columns": "columns.html", + "multivariate": "general_statistics.html", +} + +def _resolve_detail_links(quality): + if quality is None: + return + for category in quality.categories: + for check in category.checks: + if not check.detail_link: + continue + anchor = check.detail_link.lstrip("#") + check.detail_link = PAGE_MAP.get(anchor) + +def write_report(output_path: Path, general, plots, duplicates_table, + column_overviews, numeric_overviews, categorical_overviews, + top_n, quality=None): + output_path.mkdir(parents=True, exist_ok=True) + shutil.copytree(str(STATIC_DIR), str(output_path / "static"), dirs_exist_ok=True) + + _resolve_detail_links(quality) + + _render_to_file(output_path / "index.html", 'LandingPage.jinja') + _render_to_file(output_path / "numeric_data.html", 'numeric_overview.jinja', + general=general, dups=duplicates_table, quality=quality) + _render_to_file(output_path / "columns.html", 'columns.jinja', + columns=column_overviews, overview=numeric_overviews, + categorical=categorical_overviews, top_n=top_n) + _render_to_file(output_path / "general_statistics.html", 'general_statistics.jinja', + plots=plots) + + +def _render_to_file(filepath: Path, template_name: str, **context): + template = env.get_template(template_name) + with open(filepath, 'w', encoding="utf-8") as f: + f.write(template.render(**context)) \ No newline at end of file diff --git a/src/cli/utils.py b/src/cli/utils.py new file mode 100644 index 0000000..e498f79 --- /dev/null +++ b/src/cli/utils.py @@ -0,0 +1,27 @@ +import time + +from termcolor import colored + + +def _fmt_duration(seconds: float) -> str: + if seconds < 60: + return f"{seconds:.1f}s" + m, s = divmod(seconds, 60) + return f"{int(m)}m {s:.0f}s" + + +def print_step(label: str): + """Print a phase header and return a callable that closes it with timing.""" + print(f"▶ {label} …") + start = time.perf_counter() + + def done(detail: str = ""): + elapsed = _fmt_duration(time.perf_counter() - start) + suffix = f" — {detail}" if detail else "" + print(colored(f" ✓", "green"), f"{label} ({elapsed}){suffix}") + + return done + + +def info(message: str): + print(colored(f" {message}", "blue")) diff --git a/src/qc_eda/basic/__init__.py b/src/cython_wrapper/__init__.py similarity index 100% rename from src/qc_eda/basic/__init__.py rename to src/cython_wrapper/__init__.py diff --git a/src/cython_wrapper/medcouple_fast.pyx b/src/cython_wrapper/medcouple_fast.pyx new file mode 100644 index 0000000..98dfb4e --- /dev/null +++ b/src/cython_wrapper/medcouple_fast.pyx @@ -0,0 +1,176 @@ +# cython: language_level=3 + +cimport cython +import numpy as np +cimport numpy as cnp +from libc.math cimport fabs + +cnp.import_array() + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline double _sign(double x) noexcept nogil: + if x > 0: + return 1.0 + elif x < 0: + return -1.0 + return 0.0 + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline double _h_kern(double *zp, double *zm, int i, int j, + int np_, int nm, double eps) noexcept nogil: + cdef double a = zp[i] + cdef double b = zm[j] + cdef double denom = a - b + if fabs(denom) <= 2 * eps: + return _sign((np_ - 1 - i - j)) + return (a + b) / denom + + +def fast_medcouple(cnp.ndarray[double, ndim=1] x_in): + cdef int n = x_in.shape[0] + if n < 3: + return 0.0 + + cdef double eps1 = np.finfo(np.float64).eps + cdef double eps2 = np.finfo(np.float64).tiny + + # Step 1: Sort descending + cdef cnp.ndarray[double, ndim=1] x = np.sort(x_in)[::-1].copy() + + # Step 1.2: Median + cdef double median + if n % 2 == 0: + median = (x[n // 2 - 1] + x[n // 2]) / 2.0 + else: + median = x[n // 2] + + cdef double x_eps = eps1 * (eps1 + fabs(median)) + + if fabs(x[0] - median) < x_eps: + return -1.0 + if fabs(x[n - 1] - median) < x_eps: + return 1.0 + + # Scale for numerical stability + cdef double scale = 2.0 * max(x[0] - median, median - x[n - 1]) + if scale == 0: + return 0.0 + + # Step 2: Build Z+ and Z- by value comparison + cdef list zp_list = [] + cdef list zm_list = [] + cdef int i, j + + for i in range(n): + if x[i] >= median - x_eps: + zp_list.append((x[i] - median) / scale) + if x[i] <= median + x_eps: + zm_list.append((x[i] - median) / scale) + + cdef int np_ = len(zp_list) + cdef int nm = len(zm_list) + + if np_ == 0 or nm == 0: + return 0.0 + + cdef cnp.ndarray[double, ndim=1] zp = np.array(zp_list, dtype=np.float64) + cdef cnp.ndarray[double, ndim=1] zm = np.array(zm_list, dtype=np.float64) + + cdef double *zp_ptr = &zp[0] + cdef double *zm_ptr = &zm[0] + + # Steps 3-4: Johnson & Mizoguchi iteration + cdef cnp.ndarray[long, ndim=1] left = np.zeros(np_, dtype=np.int64) + cdef cnp.ndarray[long, ndim=1] right = np.full(np_, nm - 1, dtype=np.int64) + + cdef long l_total = 0 + cdef long r_total = np_ * nm + cdef long mc_idx = r_total // 2 + + cdef long mid, r_tent_total, l_tent_total + cdef double wm, wm_eps, h + cdef int n_mid + + cdef cnp.ndarray[double, ndim=1] row_med = np.empty(np_, dtype=np.float64) + cdef cnp.ndarray[double, ndim=1] wts = np.empty(np_, dtype=np.float64) + cdef cnp.ndarray[long, ndim=1] r_tent = np.empty(np_, dtype=np.int64) + cdef cnp.ndarray[long, ndim=1] l_tent = np.empty(np_, dtype=np.int64) + + while r_total - l_total > np_: + + # Weighted median of row medians + n_mid = 0 + for i in range(np_): + if left[i] <= right[i]: + mid = (left[i] + right[i]) // 2 + row_med[n_mid] = _h_kern(zp_ptr, zm_ptr, i, mid, np_, nm, eps2) + wts[n_mid] = right[i] - left[i] + 1 + n_mid += 1 + + if n_mid == 0: + break + + order = np.argsort(row_med[:n_mid]) + cum = np.cumsum(wts[:n_mid][order]) + half = cum[n_mid - 1] / 2.0 + wm = row_med[order[0]] + for i in range(n_mid): + if cum[i] >= half: + wm = row_med[order[i]] + break + + wm_eps = eps1 * (eps1 + fabs(wm)) + + # P[i]: right tentative border + j = 0 + for i in range(np_ - 1, -1, -1): + while j < nm and _h_kern(zp_ptr, zm_ptr, i, j, np_, nm, eps2) - wm > wm_eps: + j += 1 + r_tent[i] = j - 1 + + # Q[i]: left tentative border + j = nm - 1 + for i in range(np_): + while j >= 0 and _h_kern(zp_ptr, zm_ptr, i, j, np_, nm, eps2) - wm < -wm_eps: + j -= 1 + l_tent[i] = j + 1 + + r_tent_total = 0 + l_tent_total = 0 + for i in range(np_): + r_tent_total += r_tent[i] + 1 + l_tent_total += l_tent[i] + + if mc_idx <= r_tent_total - 1: + for i in range(np_): + right[i] = r_tent[i] + r_total = r_tent_total + elif mc_idx > l_tent_total - 1: + for i in range(np_): + left[i] = l_tent[i] + l_total = l_tent_total + else: + return wm + + # Step 5: Collect the remaining and select + remaining = [] + for i in range(np_): + for j in range(left[i], right[i] + 1): + remaining.append(_h_kern(zp_ptr, zm_ptr, i, j, np_, nm, eps2)) + + if len(remaining) == 0: + return 0.0 + + cdef cnp.ndarray[double, ndim=1] rem = np.array(remaining, dtype=np.float64) + rem.sort() + cdef long sel = mc_idx - l_total + if sel < 0: + sel = 0 + if sel >= len(rem): + sel = len(rem) - 1 + + return rem[sel] \ No newline at end of file diff --git a/src/qc_eda/basic/taxonomy_validator.pyx b/src/cython_wrapper/taxonomy_validator.pyx similarity index 100% rename from src/qc_eda/basic/taxonomy_validator.pyx rename to src/cython_wrapper/taxonomy_validator.pyx diff --git a/src/qc_eda/basic/wrapper_utils.pyx b/src/cython_wrapper/wrapper_utils.pyx similarity index 100% rename from src/qc_eda/basic/wrapper_utils.pyx rename to src/cython_wrapper/wrapper_utils.pyx diff --git a/src/data_utils/__init__.py b/src/data_utils/__init__.py new file mode 100644 index 0000000..9d3f6e8 --- /dev/null +++ b/src/data_utils/__init__.py @@ -0,0 +1,2 @@ +from .file_reader import read_file, parse_parquet +from .remote_data import get_tax_ids, get_clusters_of_orthologous_groups, get_gene_ontology \ No newline at end of file diff --git a/src/utils/file_reader.py b/src/data_utils/file_reader.py similarity index 59% rename from src/utils/file_reader.py rename to src/data_utils/file_reader.py index 8ee2b9e..43aacba 100644 --- a/src/utils/file_reader.py +++ b/src/data_utils/file_reader.py @@ -7,17 +7,25 @@ import click def read_file(file: click.Path) -> pd.DataFrame | None: - file = pathlib.Path(file).absolute() + file = pathlib.Path(file.__str__()).absolute() ext = pathlib.Path(file.__str__()).suffix seperator = _get_sep(file) with open(file, encoding="utf-8") as csv_file: csv_bytes = "".join(csv_file.readline() for _ in range(1000)) - dialect = csv.Sniffer().sniff(csv_bytes) - header = csv.Sniffer().has_header(csv_bytes) + try: + dialect = csv.Sniffer().sniff(csv_bytes) + header = csv.Sniffer().has_header(csv_bytes) + except csv.Error: + # Fallback for single-column or un-sniffable files + class DummyDialect: + delimiter = seperator + dialect = DummyDialect() + header = True csv_file.seek(0) - if dialect.delimiter is not seperator: - dialect.delimiter=seperator + + if dialect.delimiter != seperator: + dialect.delimiter = seperator head_col = 0 idx_col = None else: @@ -36,19 +44,24 @@ def read_file(file: click.Path) -> pd.DataFrame | None: return df elif ext == ".json": - df = pd.read_json(file.__str__(), orient='values') - cols = [i for i in df.columns if isinstance(df[i][0], dict)] - if not cols: - df = df.T - cols = [i for i in df.columns if isinstance(df[i][0], dict)] + # Removed orient='values' to allow pandas to natively infer dictionaries vs arrays + df = pd.read_json(file.__str__()) + + if not df.empty: + # Use .iloc[0] to look up positionally and avoid KeyErrors + cols = [i for i in df.columns if isinstance(df[i].iloc[0], dict)] if not cols: - return df - else: - data_frames = list() - for i in cols: - tmp = pd.json_normalize(df[i]) - data_frames.append(tmp) - df = reduce(lambda left, right: pd.merge(left, right, left_index=True, right_index=True), data_frames) + df = df.T + cols = [i for i in df.columns if isinstance(df[i].iloc[0], dict)] + if not cols: + return df + else: + data_frames = list() + for i in cols: + tmp = pd.json_normalize(df[i]) + data_frames.append(tmp) + df = reduce(lambda left, right: pd.merge(left, right, left_index=True, right_index=True), + data_frames) return df return None @@ -74,5 +87,5 @@ def parse_parquet(file: click.Path) -> pd.DataFrame: if not all([isinstance(col, str) for col in df.columns]): df.columns = df.columns.map(str) if '' in df.columns or None in df.columns: - df = df.rename({"": "Unknown", None: "Unknown_None"}) + df = df.rename(columns={"": "Unknown", None: "Unknown_None"}) return df diff --git a/src/data_utils/remote_data.py b/src/data_utils/remote_data.py new file mode 100644 index 0000000..69add93 --- /dev/null +++ b/src/data_utils/remote_data.py @@ -0,0 +1,119 @@ +import io +import os +import pathlib +import zipfile +from pathlib import Path + +import pandas as pd +import requests +from goatools.base import download_go_basic_obo +from goatools.obo_parser import GODag +import time +from pathlib import Path + +BPK_CACHE_ROOT = Path(os.environ.get("BPK_CACHE_DIR", Path.cwd() / ".bioprofilekit")) + +CACHE_DIR = BPK_CACHE_ROOT / "taxonomy" +#CACHE_DIR = Path.cwd() / ".bioprofilekit" / "taxonomy" #ToDo change to Home --> Backend +CACHE_TIL_DAYS = 30 +TAXONOMY_FILE = "taxonomy_raw.parquet" +TAXONOMY_VOCAB = "taxonomy_vocab.parquet" + + +def get_gene_ontology(): + obo_path = download_go_basic_obo() + go_dag = GODag(obo_path) + print(go_dag.version) + data = list() + for go_id in go_dag.keys(): + term = go_dag[go_id] + namespace = getattr(term, "namespace", "") + data.append([go_id, term.name, namespace]) + + df = pd.DataFrame(data, columns=["GO_ID", "Name", "Namespace"]) + if Path(obo_path).is_file(): + pathlib.Path(obo_path).unlink(missing_ok=True) + print(f"Removed {obo_path}") + return df + + +def get_clusters_of_orthologous_groups(): + url: str = "https://ftp.ncbi.nlm.nih.gov/pub/COG/COG2024/data/cog-24.def.tab" + fields = ["COG_ID", "Functional Category", "COG name"] + response = requests.get(url) + + if response.status_code == 200: + df = pd.read_csv(io.StringIO(response.text), sep="\t", skipinitialspace=True, usecols=[0, 1, 2], names=fields) + else: + print(f"Error: {response.status_code}") + return df + + +def get_tax_ids(force_refresh: bool = False): + CACHE_DIR.mkdir(parents=True, exist_ok=True) + raw_path = CACHE_DIR / TAXONOMY_FILE + vocab_path = CACHE_DIR / TAXONOMY_VOCAB + + if not force_refresh and vocab_path.is_file(): + age_days = (time.time() - vocab_path.stat().st_mtime) / 86400 + if age_days < CACHE_TIL_DAYS: + return pd.read_parquet(vocab_path) + + """if not force_refresh and raw_path.is_file(): + age_days = (time.time() - vocab_path.stat().st_mtime) / 86400 + if age_days < CACHE_TIL_DAYS: + raw = pd.read_parquet(raw_path)""" + + raw = _download_taxonomy() + raw.to_parquet(raw_path, index=False) + + vocab = build_taxonomy(raw) + vocab.to_parquet(vocab_path, index=False) + return vocab + +def _download_taxonomy(): + url = "https://ftp.ncbi.nih.gov/pub/taxonomy/taxdmp.zip" + + print(f"Downloading {url} ...") + resp = requests.get(url, stream=True) + resp.raise_for_status() + + with zipfile.ZipFile(io.BytesIO(resp.content)) as zf: + print("Files inside ZIP:", zf.namelist()) + + with zf.open("names.dmp") as fh: + print(fh) + names = pd.read_csv( + fh, + sep="|", + header=None, + index_col=False, + names=["tax_id", "name_txt", "unique_name", "name_class"], + engine="c" + ) + + with zf.open("nodes.dmp") as fh: + nodes = pd.read_csv(fh, sep="|", header=None, index_col=False, usecols=[0,2], names=["tax_id", "rank"], engine="c") + + names = names.map(lambda x: x.strip() if isinstance(x, str) else x) + nodes = nodes.map(lambda x: x.strip() if isinstance(x, str) else x) + df = names.merge(nodes, on="tax_id", how="left") + + return df + +def build_taxonomy(tax_df: pd.DataFrame) -> pd.DataFrame: + name_classes = ['scientific name', 'synonym', 'equivalent name','genbank common name', 'common name'] + sci = (tax_df[tax_df['name_class'] == 'scientific name'].drop_duplicates('tax_id').set_index('tax_id')['name_txt']) + names = tax_df[tax_df['name_class'].isin(name_classes)][['tax_id', 'name_txt', 'name_class', 'rank']].copy() + strains = tax_df[tax_df['name_class'] == 'type material'][['tax_id', 'name_txt']] + strains = strains.rename(columns={'name_txt': 'strain'}) + + combos = names.merge(strains, on="tax_id", how="inner") + combos['name_txt'] = combos['name_txt'] + ' ' + combos['strain'] + combos['name_class']= 'type_strain' + combos = combos[['tax_id', 'name_txt', 'name_class', 'rank']] + + vocab = pd.concat([names, combos], ignore_index=True) + vocab['scientific_name'] = vocab['tax_id'].map(sci) + vocab = vocab.drop_duplicates().reset_index(drop=True) + return vocab.sort_values(by=['tax_id'], ascending=True) diff --git a/src/enums/__init__.py b/src/enums/__init__.py new file mode 100644 index 0000000..03d1441 --- /dev/null +++ b/src/enums/__init__.py @@ -0,0 +1,2 @@ +from .sequence_enum import Sequence +from .measurement_enum import MEASUREMENTS diff --git a/src/qc_eda/basic/measurement_enum.py b/src/enums/measurement_enum.py similarity index 57% rename from src/qc_eda/basic/measurement_enum.py rename to src/enums/measurement_enum.py index 6d719f4..7278d8c 100644 --- a/src/qc_eda/basic/measurement_enum.py +++ b/src/enums/measurement_enum.py @@ -2,7 +2,7 @@ from enum import Enum class MEASUREMENTS(Enum): - UNIT_IN_COL_TITLE = re.compile(r'\[[a-zA-Zµμ°/]+]', re.I) + UNIT_IN_COL_TITLE = re.compile(r'\[[a-zA-Zµμ°/%]+]', re.I) UNIT_PREFIXES = r'(Q|R|Y|Z|E|P|T|G|M|k|h|da|d|c|m|µ|μ|micro|n|p|f|a|z|y|r|q)' - UNITS = r'(g|mol|M|m|L|l|s|min|unit|units|°C|C|K|A|cd|h|24 h|d|g creatinine|Eq|%|U|kU|IU|L (%)|% of total Hb|Fraction of total Hb|Osm|AU|angstroms|\d)' - UNIT_COLUMN = re.compile(rf'^((\d\.,)+\s*)?{UNIT_PREFIXES}?{UNITS}(/{UNIT_PREFIXES}?{UNITS})*$', re.I) + UNITS = r'(g|mol|M|m|L|l|s|min|unit|units|°C|C|K|A|cd|h|24 h|d|g creatinine|Eq|%|U|kU|IU|L (%)|% of total Hb|Fraction of total Hb|Osm|AU|angstroms)' + UNIT_COLUMN = re.compile(rf'^((\d\.,)+\s+)?{UNIT_PREFIXES}?({UNITS}|\d)(/{UNIT_PREFIXES}?{UNITS})*$', re.I) diff --git a/src/qc_eda/basic/sequence_enum.py b/src/enums/sequence_enum.py similarity index 58% rename from src/qc_eda/basic/sequence_enum.py rename to src/enums/sequence_enum.py index 86e98cb..9bf68a8 100644 --- a/src/qc_eda/basic/sequence_enum.py +++ b/src/enums/sequence_enum.py @@ -2,6 +2,6 @@ from enum import Enum class Sequence(Enum): - DNA = re.compile('^[acgtn]*$', re.I) - RNA = re.compile('^[acgun]*$', re.I) + DNA = re.compile('^[-acgtn]+$', re.I) + RNA = re.compile('^[-acgun]*$', re.I) PROTEIN = re.compile('^[acdefghiklmnpqrstvwyxju]*$', re.I) \ No newline at end of file diff --git a/src/models/__init__.py b/src/models/__init__.py new file mode 100644 index 0000000..68d9847 --- /dev/null +++ b/src/models/__init__.py @@ -0,0 +1,8 @@ +from .overview import DatasetSummary, ColumnOverview +from .numeric import NumericColumns +from .categorical import CategoricalColumns +from .sequence import DNARNAColumns, PROTEINColumns +from .measurement import UNITColumns +from .outliers import Outliers +from .multivariate import MultivariateAnalysis +from .quality import QualityCategory, QualityCheck, QualityAssessment \ No newline at end of file diff --git a/src/models/categorical.py b/src/models/categorical.py new file mode 100644 index 0000000..1fbc600 --- /dev/null +++ b/src/models/categorical.py @@ -0,0 +1,22 @@ +from dataclasses import dataclass + + +@dataclass +class CategoricalColumns: + name: str + unique_categories: int + mode: str + entropy: float + frequencies: dict + gini: float + simpson_diversity: float + value_counts: dict + max_category_length: int + min_category_length: int + cardinality_ratio: float + rare_categories: float + top_5_coverage: float + cib_ratio: float + top_1_coverage: float + effective_cardinality: float + memory: int diff --git a/src/models/measurement.py b/src/models/measurement.py new file mode 100644 index 0000000..6a338c1 --- /dev/null +++ b/src/models/measurement.py @@ -0,0 +1,9 @@ +from dataclasses import dataclass +from typing import List, Dict + + +@dataclass +class UNITColumns: + units: List[str] + unit_counts: Dict[str | None, int] + with_measurement: bool diff --git a/src/models/multivariate.py b/src/models/multivariate.py new file mode 100644 index 0000000..1f1ed58 --- /dev/null +++ b/src/models/multivariate.py @@ -0,0 +1,30 @@ +from dataclasses import dataclass +from typing import Optional + +import pandas as pd + +@dataclass +class MultivariateAnalysis: + # Plots (HTML strings) + correlation_heatmap: str + pearson_heatmap: Optional[str] # nur numerisch, -1..1 + cramers_heatmap: Optional[str] # nur kategorisch, 0..1 + eta_heatmap: Optional[str] + missing_matrix: str + missing_values_barchart: str + balance_plot: str | None + boxplot: str + scatter_matrix: str + mutual_information_plot: str | None + mi_relationship_plots: list | None + + # Computed metrics + correlation_matrix: pd.DataFrame # combined, absolute values + correlation_methods: pd.DataFrame # method used per pair + top_associations: list | None + + feature_target_correlation: dict | None # needs target + feature_target_plot: str | None # needs target + mutual_information: dict | None # needs target + mcar_result: dict | None # ToDo replace with finale Littles MCAR + target_name: str | None = None diff --git a/src/models/numeric.py b/src/models/numeric.py new file mode 100644 index 0000000..fe9b943 --- /dev/null +++ b/src/models/numeric.py @@ -0,0 +1,30 @@ +from dataclasses import dataclass +from typing import Optional + +from numpy import ndarray + +from .outliers import Outliers + + +@dataclass +class NumericColumns: + name: str + min: float + max: float + mean: float + median: float + mode: float + std: float + sum: float + kurtosis: float + skewness: float + coefficient_of_variation: float + mad: float + quantiles: ndarray + infinity: float + negative_count: float + zero_count: float + memory: int + value_counts: dict + frequencies: dict + outliers: Optional[Outliers] diff --git a/src/models/outliers.py b/src/models/outliers.py new file mode 100644 index 0000000..23e6d72 --- /dev/null +++ b/src/models/outliers.py @@ -0,0 +1,16 @@ +from dataclasses import dataclass + +import numpy as np + + +@dataclass +class Outliers: + lower_bound: np.float64 + upper_bound: np.float64 + n_lower_iqr: int + n_upper_iqr: int + medcouple: np.float64 + n_lower_mzscore: int + n_upper_mzscore: int + n_lower_zscore: int + n_upper_zscore: int \ No newline at end of file diff --git a/src/models/overview.py b/src/models/overview.py new file mode 100644 index 0000000..d6cd350 --- /dev/null +++ b/src/models/overview.py @@ -0,0 +1,49 @@ +from dataclasses import dataclass +from typing import Optional + + +@dataclass +class DatasetSummary: + filename: str + rows: int + cols: int + nulls: int + nulls_percentage: float + empty_rows: int + dup_row: int + dup_col: int + ratio: float + memory: float + alerts: int + n_number: int = 0 + n_dna: int = 0 + n_rna: int = 0 + n_protein: int = 0 + n_taxonomy: int = 0 + n_unit: int = 0 + n_functional: int = 0 + n_categorical: int = 0 + n_empty: int = 0 + + +@dataclass +class ColumnOverview: + name: str + number: int | None + unique: int | None + missing: int | None + missing_per: float | None + density: float | None + type: str + sequence: str | None + invalid_seqs: list[str] | None + mixed_types: list[str] | None + suspect_values: list[str] | None + describe_plot: str | None + constant: bool | None + correlation: list[str] | None + cardinality_dimension_ratio: float | None + monotonicity: bool | None + taxonomy: Optional = None + measurement_data: Optional = None + taxonomy_candidate: bool = False diff --git a/src/models/quality.py b/src/models/quality.py new file mode 100644 index 0000000..c0cb03d --- /dev/null +++ b/src/models/quality.py @@ -0,0 +1,26 @@ +from dataclasses import dataclass + + +@dataclass +class QualityCheck: + name: str + status: str # pass | warn | fail + message: str + detail_link: str | None + + +@dataclass +class QualityCategory: + name: str + status: str + checks: list # List[QualityCheck] + + +@dataclass +class QualityAssessment: + categories: list # List[QualityCategory] + passed: int + warnings: int + failed: int + total: int + overall: str diff --git a/src/models/sequence.py b/src/models/sequence.py new file mode 100644 index 0000000..86a5cd3 --- /dev/null +++ b/src/models/sequence.py @@ -0,0 +1,77 @@ +from dataclasses import dataclass +from typing import List, Dict, Tuple, Optional + +from .outliers import Outliers + + +@dataclass +class SequenceMetricSummary: + min: float + max: float + mean: float + +@dataclass +class DNARNAColumns: + # Per-sequence (Top-20) + sequence: List[str] + count: List[int] + length: List[int] + + # Column-wide summaries + gc_content: SequenceMetricSummary + ambiguous_base_ratio: SequenceMetricSummary + length_stats: SequenceMetricSummary + length_outliers:Optional[Outliers] + codon_completeness: SequenceMetricSummary + gc_skew: SequenceMetricSummary + at_skew: SequenceMetricSummary + cpg_observed_expected: SequenceMetricSummary + tpa_observed_expected: SequenceMetricSummary + low_complexity: SequenceMetricSummary + reverse_complement_ratio: float + reverse_complement_list: set[str] + # ToDo: Check if needed from Top N or overall + nucleotide_count: List[Dict[str, int]] + k_mers: List[List[Tuple[str, int]]] + + + # Plots + plot: str + gc_distribution: Optional[str] + length_distribution: Optional[str] + ambiguous_distribution: Optional[str] + at_gc_skewness: Optional[str] + +# ToDo: Add Composition over all +@dataclass +class PROTEINColumns: + sequence: List[str] + length: List[int] + count: List[int] + composition: List[Dict[str, int]] + frequency: List[float] + hydrophobicity: List[float] + charge: List[float] + molecular_weight: List[float] + isoelectric_point: List[float] + aliphatic_index: List[float] + boman: List[float] + aromaticity: List[float] + instability: List[float] + + # Column-wide metrics + ambiguous_residue_ratio: SequenceMetricSummary + length_stats: SequenceMetricSummary + length_outliers: Optional[Outliers] + stop_codon_ratio: float + low_complexity: SequenceMetricSummary + gravy: SequenceMetricSummary + cysteine_count: SequenceMetricSummary + disorder_propensity: SequenceMetricSummary + aa_group_distribution: Dict[str, float] + + k_mers: List[List[Tuple[str, int]]] + plot: str + length_distribution: Optional[str] + ambiguous_distribution: Optional[str] + aa_group_plot: Optional[str] diff --git a/src/qc_eda/basic/general.py b/src/qc_eda/basic/general.py deleted file mode 100644 index 1bd13d0..0000000 --- a/src/qc_eda/basic/general.py +++ /dev/null @@ -1,161 +0,0 @@ -import plotly.express as px -import pandas as pd -import plotly.graph_objects as go -from dataclasses import dataclass - -@dataclass -class GeneralPlots: - correlation_heatmap: str - missing_matrix: str - missing_values_barchart: str - balance_plot: str | None - boxplot: str - scatter_matrix: str - -def general_plots(df: pd.DataFrame, target: str) -> GeneralPlots: - return GeneralPlots( - correlation_heatmap=correlation_heatmap(df), - missing_matrix=missing_matrix(df), - missing_values_barchart=missing_values_barchart(df), - balance_plot=balance_plot(df, target) if target else None, - boxplot=boxplot(df), - scatter_matrix=scatter_matrix(df) - ) - -def correlation_heatmap(df: pd.DataFrame): - df_numeric = df.select_dtypes(include=['float64', 'int64']).dropna(axis=1, how='all') - if not df_numeric.empty: - std = df_numeric.std(ddof=0) - df_numeric = df_numeric.loc[:, std[std > 0].index] - corr_matrix = df_numeric.corr() - corr_matrix = round(corr_matrix, 3) - fig = px.imshow(corr_matrix, text_auto=True, labels=dict(color="Correlation"), color_continuous_scale="RdBu_r", aspect="auto", height=700) - fig.update_layout(title="Correlation Heatmap") - return fig.to_html(full_html=False, include_plotlyjs=False) - - -def missing_matrix(df: pd.DataFrame): - missing_values = df.isnull().astype(int) - - fig = px.imshow( - missing_values, - labels=dict(color="Missing Values"), - aspect="auto", - color_continuous_scale="blues_r", - title="Missing Values Matrix" - ) - fig.update_yaxes(autorange='reversed') - fig.update_layout(coloraxis_showscale=False) - - - return fig.to_html(full_html=False, include_plotlyjs=False) - -def missing_values_barchart(df: pd.DataFrame): - missing_counts = df.isna().sum() - fig = px.bar(x=missing_counts.index, y=missing_counts.values, labels={'x': 'Columns', 'y': 'Missing Values'},color_discrete_sequence=['#0F65A0']) - fig.update_layout(title="Missing Values per Column", bargap=0.2, plot_bgcolor='white') - fig.update_xaxes( - tickangle=-45, - mirror=True, - ticks='outside', - showline=True, - linecolor='black', - gridcolor='lightgrey' - ) - fig.update_yaxes( - mirror=True, - ticks='outside', - showline=True, - linecolor='black', - gridcolor='lightgrey' - ) - - return fig.to_html(full_html=False, include_plotlyjs=False) - -def balance_plot(df, target): - fig = px.histogram(df, x=target, color_discrete_sequence=["#0F65A0"], text_auto=True) - fig.update_layout( - title="Class Balance (Target Distribution)", - bargap=0.2, - plot_bgcolor="white" - ) - fig.update_xaxes( - mirror=True, - ticks="outside", - showline=True, - linecolor="black", - gridcolor="lightgrey" - ) - fig.update_yaxes( - mirror=True, - ticks="outside", - showline=True, - linecolor="black", - gridcolor="lightgrey" - ) - return fig.to_html(full_html=False, include_plotlyjs=False) - -def boxplot(df: pd.DataFrame): - df = df.select_dtypes(include=['float64', 'int64']).dropna(axis=1, how='all') - fig = go.Figure() - - for col in df: - fig.add_trace(go.Box(y=df[col].values, name=df[col].name)) - fig.update_yaxes(type="log", title="Logarithmic",showticklabels=False) - fig.update_layout( - xaxis=dict(rangeslider=dict(visible=True), type="linear"), - title="Boxplot", - xaxis_title="Columns", - yaxis_title="Values", - legend_title="Columns" - ) - - return fig.to_html(full_html=False, include_plotlyjs=False) - -def scatter_matrix(df: pd.DataFrame): - df = df.select_dtypes(include=['float64', 'int64']).dropna(axis=1, how='all') - fig = px.scatter_matrix(df, color_discrete_sequence=["#0F65A0"], height=750) - - fig.update_traces( - diagonal_visible=False, - marker=dict(size=2, opacity=0.5, color="#0F65A0") - ) - - # Layout anpassen - fig.update_layout( - title="Scatter Matrix", - xaxis_title="Columns", - plot_bgcolor="white", - bargap=0.2, - dragmode="select" - ) - n_vars = len(df.columns) - for i in range(1, n_vars + 1): - for j in range(1, n_vars + 1): - if i == 1: - fig.update_layout({f'xaxis{j if j > 1 else ""}': dict( - mirror=True, - ticks="outside", - showline=True, - linecolor="black", - linewidth=1, - gridcolor="lightgrey", - title_standoff=25 - )}) - if j == 1: - fig.update_layout({f'yaxis{i if i > 1 else ""}': dict( - mirror=True, - ticks="outside", - showline=True, - linecolor="black", - linewidth=1, - gridcolor="lightgrey", - title_standoff=25 - )}) - fig.for_each_annotation(lambda a: a.update( - textangle=0, - x=-0.08 if a.textangle == 90 or 'y' in str(a.yref) else a.x, - y=-0.08 if a.textangle == 0 and 'x' in str(a.xref) else a.y - )) - - return fig.to_html(full_html=False, include_plotlyjs=False) diff --git a/src/qc_eda/basic/numerical_data.py b/src/qc_eda/basic/numerical_data.py deleted file mode 100644 index 3f43512..0000000 --- a/src/qc_eda/basic/numerical_data.py +++ /dev/null @@ -1,302 +0,0 @@ -import re -from dataclasses import dataclass - -import math -import numpy as np -import pandas as pd -import plotly.express as px -from numpy import ndarray -from pandas.api.types import infer_dtype -from scipy import stats - -from .sequence_enum import Sequence -from .wrapper_utils import fast_check_sequence, char_entropy - -def _alphabet_from_pattern(pattern): - match = re.search(r'\[([^]]+)]', pattern.pattern) - return set(match.group(1).upper()) if match else set() - -# ToDo test if DNA and RNA Entropy is needed -DNA_ALPHABET = _alphabet_from_pattern(Sequence.DNA.value) -RNA_ALPHABET = _alphabet_from_pattern(Sequence.RNA.value) -PROTEIN_ALPHABET = _alphabet_from_pattern(Sequence.PROTEIN.value) - -ENTROPY_THRESHOLDS = { - "dna": 0.5 * math.log2(len(DNA_ALPHABET)), - "rna": 0.5 * math.log2(len(RNA_ALPHABET)), - "protein": 0.5 * math.log2(len(PROTEIN_ALPHABET)), -} - -@dataclass -class NumericalData: - filename: str - rows: int - cols: int - nulls: int - nulls_percentage: float - dup_row: int - dup_col: int - memory: float - alerts: int - - -@dataclass -class ColumnOverview: - name: str - number: int | None - unique: int | None - missing: int | None - missing_per: float | None - type: str - sequence: str | None - invalid_seqs: list[str] | None - describe_plot: str | None - constant: bool | None - correlation: list[str] | None - # taxonomy: bool - - -@dataclass -class NumericColumns: - name: str - min: float - max: float - mean: float - median: float - mode: float - std: float - sum: float - kurtosis: float - skewness: float - coefficient_of_variation: float - mad: float - quantiles: ndarray - memory: int - value_counts: dict - frequencies: dict - - # cardinalities: list[int] - - -@dataclass -class CategoricalColumns: - name: str - unique_categories: int - mode: str - entropy: float - frequencies: dict - gini: float - simpson_diversity: float - value_counts: dict - max_category_length: int - min_category_length: int - memory: int - cardinality_ratio: float - - -def overview(df: pd.DataFrame, file) -> NumericalData: - return NumericalData( - filename=file, - rows=df.shape[0], - cols=df.shape[1], - nulls=sum(df.isnull().sum()), - nulls_percentage=100 if df.isnull().all().all() else round(sum(df.isnull().sum()) * 100 / df.size, 2), - dup_row=int(df.duplicated().sum()), - dup_col=int(df.columns.duplicated().sum()), - memory=int(df.memory_usage(deep=True).sum()), - alerts=0 - ) - -#ToDo check for empty Column and return -def column_overview(df: pd.DataFrame, col) -> ColumnOverview: - seq_type, invalid = check_sequence(df, col) - if seq_type != "None": - print(f"Column: {col:10s} is of type: {seq_type:10s} with invalid sequences: {invalid}.") - return ColumnOverview( - name=col, - number=int(df[col].notnull().sum()), - unique=df[col].nunique(), - missing=int(df[col].isnull().sum()), - missing_per=100 if df[col].isnull().all() else round(sum(df[col].isnull()) * 100 / df[col].size, 2), - type=str(df[col].dtype), - sequence=seq_type, - invalid_seqs=invalid, - describe_plot=None if df[col].isnull().all() else plot_overview(df[col]), - constant=True if (df[col].nunique() == 1) else False, - correlation=get_correlation(df, col), - ) - -# -def numeric_columns(df: pd.DataFrame, col) -> NumericColumns: - series = df[col].dropna() - - if series.empty: - coefficient_of_variation = np.nan - quantiles = np.array([np.nan, np.nan, np.nan]) - mode_value = np.nan - value_counts = {} - frequencies = {} - mad_value = np.nan - min_value = np.nan - max_value = np.nan - mean_value = np.nan - median_value = np.nan - std_value = np.nan - sum_value = np.nan - kurtosis_value = np.nan - skewness_value = np.nan - else: - mean_value = float(series.mean()) - std_value = float(series.std()) - coefficient_of_variation = np.nan if mean_value == 0 else float(std_value / abs(mean_value)) - quantiles = np.array(series.quantile([0.25, 0.5, 0.75]).to_list(), dtype=float) - mode_series = series.mode() - mode_value = float(mode_series.iloc[0]) if not mode_series.empty else np.nan - value_counts = series.value_counts().head(20).to_dict() - frequencies = series.value_counts(normalize=True).head(20).to_dict() - mad_value = float(stats.median_abs_deviation(series, nan_policy='omit')) - min_value = float(series.min()) - max_value = float(series.max()) - median_value = float(series.median()) - sum_value = float(series.sum()) - kurtosis_value = float(series.kurtosis()) - skewness_value = float(series.skew()) - - return NumericColumns( - name=col, - min=round(min_value, 2) if np.isfinite(min_value) else np.nan, - max=round(max_value, 2) if np.isfinite(max_value) else np.nan, - mean=round(mean_value, 2) if np.isfinite(mean_value) else np.nan, - median=round(median_value, 2) if np.isfinite(median_value) else np.nan, - mode=round(mode_value, 2) if np.isfinite(mode_value) else np.nan, - std=round(std_value, 2) if np.isfinite(std_value) else np.nan, - sum=round(sum_value, 2) if np.isfinite(sum_value) else np.nan, - kurtosis=round(kurtosis_value, 2) if np.isfinite(kurtosis_value) else np.nan, - skewness=round(skewness_value, 2) if np.isfinite(skewness_value) else np.nan, - coefficient_of_variation=round(coefficient_of_variation, 2) if np.isfinite(coefficient_of_variation) else np.nan, - mad=mad_value if np.isfinite(mad_value) else np.nan, - quantiles=quantiles, - memory=df[col].memory_usage(deep=True), - value_counts=value_counts, - frequencies=frequencies - ) - - -def categorical_columns(df: pd.DataFrame, col: str) -> CategoricalColumns: - value_counts = df[col].value_counts() - n = int(df[col].notna().sum()) - frequencies = value_counts / n if n > 0 else value_counts - if n > 0: - entropy = -(frequencies * np.log2(frequencies)).sum() - gini = 1 - (frequencies ** 2).sum() - simpson = 1 / (frequencies ** 2).sum() - cardinality_ratio = df[col].nunique() / n - mode_value = df[col].mode().iloc[0] - else: - entropy = np.nan - gini = np.nan - simpson = np.nan - cardinality_ratio = np.nan - mode_value = "" - lengths = df[col].astype(str).str.len() - - return CategoricalColumns( - name=col, - unique_categories=df[col].nunique(), - mode=mode_value, - entropy=round(entropy, 2), - frequencies=df[col].value_counts(normalize=True).head(20).to_dict(), - gini=round(gini, 2), - simpson_diversity=round(simpson, 2), - value_counts=df[col].value_counts().head(20).to_dict(), - max_category_length=lengths.max(), - min_category_length=lengths.min(), - memory=df[col].memory_usage(deep=True), - cardinality_ratio=round(cardinality_ratio, 3) - ) - - -def get_correlation(df: pd.DataFrame, col) -> list | None: - ncols = df.select_dtypes(include='number').dropna(axis=1, how='all').columns - if col in ncols: - std = df[ncols].std(ddof=0) - ncols = std[std > 0].index - if col not in ncols: - return None - corr = df[ncols].corrwith(df[col], method='pearson') - corr.drop(labels=col, inplace=True) - corr = corr.drop(corr[corr < .3].index) - if corr.empty: - return None - return list(zip(corr.index, corr)) - return None - - -# ToDo: move to plot_utils -def plot_overview(col): - if col.dtype != 'object': - bins = None if col.nunique() < 10 else 10 - fig = px.histogram(col, nbins=bins, color_discrete_sequence=['#0F65A0']) - fig.update_layout(bargap=0.2, plot_bgcolor='white') - fig.update_xaxes( - mirror=True, - ticks='outside', - showline=True, - linecolor='black', - gridcolor='lightgrey' - ) - fig.update_yaxes( - mirror=True, - ticks='outside', - showline=True, - linecolor='black', - gridcolor='lightgrey' - ) - return fig.to_html(full_html=False, include_plotlyjs=False) - return None - - -# ToDo: move to sequence_utils -def check_sequence(df, col, threshold=0.95): - if df[col].name in df.select_dtypes(include=['number', 'bool']).columns or infer_dtype(df[col]).__contains__('mixed'): - return "None", [] - if df[col].astype(str).str.len().eq(1).all(): - return "None", [] - values = df[col].dropna().astype(str).tolist() - - #ToDo real cardinality - unique_count = len(set(values)) - if unique_count < 10: - return "None", [] - - if all(len(x) > 2 for x in values): - match, invalid = fast_check_sequence(values, Sequence.DNA.value, threshold) - if match: - return "dna", _get_invalid(values, invalid) - match, invalid = fast_check_sequence(values, Sequence.RNA.value, threshold) - if match: - return "rna", _get_invalid(values, invalid) - match, invalid = fast_check_sequence(values, Sequence.PROTEIN.value, threshold) - if match: - if not invalid or char_entropy(values, PROTEIN_ALPHABET) >= ENTROPY_THRESHOLDS["protein"]: - return "protein", _get_invalid(values, invalid) - return "None", [] - -def _get_invalid(values, invalid_indices): - if not invalid_indices: - return [] - return [values[i] for i in invalid_indices] - - -"""def rank_taxonomy(df, col): - if df[col].dtype != 'object': - return False - - results = df[col].astype(str).apply(lambda x: validate_taxonomy(x)) - results = results[~results.str.len().eq(0)] - - if not results.empty: - results = set(deepflatten(results.value_counts().index.tolist(), depth=1)) - print(results) - - return False""" diff --git a/src/qc_eda/biological/biological_data.py b/src/qc_eda/biological/biological_data.py deleted file mode 100644 index 2fecb06..0000000 --- a/src/qc_eda/biological/biological_data.py +++ /dev/null @@ -1,222 +0,0 @@ -import tempfile -import warnings -from collections import defaultdict, Counter -from dataclasses import dataclass -from itertools import chain -from pathlib import Path -from typing import Tuple -from Bio import motifs -from weblogo import * -import numpy as np -import pandas as pd -import peptides -import plotly.express as px -import ssl -ssl._create_default_https_context = ssl._create_stdlib_context - -@dataclass -class DNARNAColumns: - sequence: List[str] - gc_content: List[float] - length: List[int] - count: List[int] - nucleotide_count: List[Dict[str, int]] - k_mers: List[List[Tuple[str, int]]] - plot: str - -#ToDo: Add Composition over all -@dataclass -class PROTEINColumns: - sequence: List[str] - length: List[int] - count: List[int] - composition: List[Dict[str, int]] - frequency: List[float] - hydrophobicity: List[float] - charge: List[float] - molecular_weight: List[float] - isoelectric_point: List[float] - aliphatic_index: List[float] - boman: List[float] - aromaticity: List[float] - instability: List[float] - k_mers: List[List[Tuple[str, int]]] - plot: str - - -def count_nmer(sequence, n) -> defaultdict: - seqs = np.frombuffer(sequence.encode('utf-8'), dtype='S1') - nmere = np.lib.stride_tricks.sliding_window_view(seqs, n) - unique, counts = np.unique(nmere, return_counts=True, axis=0) - result = defaultdict(int) - for nmer, count in zip(unique, counts): - nmer_str = b"".join(nmer).decode('utf-8') - result[nmer_str] += count - return result - - -def top_mere(seq, n=3, top=5) -> List[Tuple[str, int]] | None: - if not seq or len(seq) < n: - return None - counts = count_nmer(seq, n) - return sorted(counts.items(), key=lambda x: x[1], reverse=True)[:top] - - -def biological_data_top_entries(seqs: pd.Series, top_k: int = 20) -> Tuple[np.ndarray, np.ndarray, int, int, np.ndarray]: - seqs =seqs.str.upper() - vc = seqs.value_counts() - uniq_tmp, counts_tmp =vc.index.to_numpy(), vc.values - - top_k = min(top_k, len(uniq_tmp)) - top_idx = np.argsort(counts_tmp)[::-1][:top_k] - - uniques = uniq_tmp[top_idx].astype(str) - counts = counts_tmp[top_idx] - - lengths = np.array([len(s) for s in uniques]) - min_len, max_len = lengths.min(), lengths.max() - - return uniques, counts, min_len, max_len, lengths - - -def dna_rna_columns(seqs: pd.Series, k: int = 3, top_n: int = 20, top: int = 5) -> DNARNAColumns: - uniques, counts, min_len, max_len, lengths = biological_data_top_entries(seqs, top_n) - - gc_count = np.char.count(uniques, 'G') + np.char.count(uniques, 'C') - gc_content = np.round(np.where(lengths > 0, gc_count / lengths * 100, 0.0), 2).tolist() - nucleotide_count = [dict(Counter(seq)) for seq in uniques] - - k_mers = _kmer_check(k, top, uniques) - - if min_len == max_len: - plot = make_logo(uniques,'color_classic', seq_type="dna") - else: - flat_kmers = chain.from_iterable(k_mers) - df_kmers = pd.DataFrame(flat_kmers, columns=['kmer', 'count']) - aggregated = df_kmers.groupby('kmer', as_index=False)['count'].sum() - plot = plot_overview(aggregated['kmer'].tolist(), aggregated['count'].tolist()) - - return DNARNAColumns( - sequence=uniques.tolist(), - gc_content=gc_content, - length=lengths.tolist(), - count=counts.tolist(), - nucleotide_count=nucleotide_count, - k_mers=k_mers, - plot=plot - ) - -def _kmer_check(k: int, top: int, uniques: np.ndarray) -> list: - check_length = any(len(i) <= k for i in uniques) - if not check_length: - k_mers = [top_mere(seq, n=k, top=top) for seq in uniques] - else: - warnings.warn("Warning: Sequence length is smaller than choosen k-Mer Size. Setting k-Mer Size to 3.") - k_mers = [top_mere(seq, n=3, top=top) for seq in uniques] - return k_mers - - -def protein_descriptors(peptide: str) -> Dict[str, str | float | dict[str, float]]: - descriptors: Dict[str, str | float | dict[str, float]] = {} - p: peptides.Peptide = peptides.Peptide(peptide) - descriptors["seq"] = peptide - descriptors["freq"] = p.frequencies() - - try: - descriptors["aidx"] = p.aliphatic_index() - except ZeroDivisionError: - descriptors["aidx"] = 0.0 - descriptors["boman"] = p.boman() - descriptors["charge"] = p.charge() - descriptors["hp"] = p.hydrophobicity() - descriptors["iep"] = p.isoelectric_point() - descriptors["iidx"] = p.instability_index() - descriptors["mol"] = p.molecular_weight() - descriptors["aroma"] = sum([peptide.count(aa) for aa in ('F', 'W', 'Y')]) / len(peptide) - return descriptors - - -def protein_columns(seqs: pd.Series, k: int = 3, top_n: int = 20, top: int = 5) -> PROTEINColumns: - uniques, counts, min_len, max_len, lengths = biological_data_top_entries(seqs, top_n) - - aa_composition = [dict(Counter(seq)) for seq in uniques] - descriptors = [protein_descriptors(seq) for seq in uniques] - - k_mers = _kmer_check(k, top, uniques) - #ToDo Add Desclaimer - if min_len == max_len: - plot = make_logo(uniques, "chemistry", seq_type="protein") - else: - flat_kmers = chain.from_iterable(k_mers) - df_kmers = pd.DataFrame(flat_kmers, columns=['kmer', 'count']) - aggregated = df_kmers.groupby('kmer', as_index=False)['count'].sum() - plot = plot_overview(aggregated['kmer'].tolist(), aggregated['count'].tolist()) - - - return PROTEINColumns( - sequence=uniques.tolist(), - length=lengths.tolist(), - count=counts.tolist(), - composition=aa_composition, - frequency=[descriptor['freq'] for descriptor in descriptors], - hydrophobicity=[descriptor['hp'] for descriptor in descriptors], - charge=[descriptor['charge'] for descriptor in descriptors], - molecular_weight=[descriptor['mol'] for descriptor in descriptors], - isoelectric_point=[descriptor['iep'] for descriptor in descriptors], - aliphatic_index=[descriptor['aidx'] for descriptor in descriptors], - boman=[descriptor['boman'] for descriptor in descriptors], - aromaticity=[descriptor['aroma'] for descriptor in descriptors], - instability=[descriptor['iidx'] for descriptor in descriptors], - k_mers=k_mers, - plot=plot - ) - - -def make_logo(seqs, color, seq_type): - if seq_type == "protein": - m = motifs.create(seqs, alphabet="ACDEFGHIKLMNPQRSTVWY") - else: - m = motifs.create(seqs, alphabet="ACGT") - with tempfile.NamedTemporaryFile(suffix='.svg', delete=False) as tmp_file: - tmp_path = tmp_file.name - - try: - m.weblogo(tmp_path, format="svg",sequence_type=seq_type, color=color,logo_font="Calibri",logo_margin=3, fontsize=12) - - with open(tmp_path, 'r', encoding='utf-8') as svg_file: - svg_content = svg_file.read() - if ' UNITColumns | bool: - if column_overview.sequence == 'None': - name_match = MEASUREMENTS.UNIT_IN_COL_TITLE.value.match(column_overview.name) - if name_match: - unit: str = name_match.group(0).lstrip('[').rstrip(']') - return UNITColumns( - units=[unit], - unit_counts=[{unit: 1}], - with_measurement=False - ) - if column_overview.type == 'object': - values = df[column_overview.name].dropna().unique().astype(str).tolist() - if all(len(x) > 1 for x in values): - measurement_and_unit = match_units(values, MEASUREMENTS.UNIT_COLUMN.value) - if len(measurement_and_unit) > 0: - return UNITColumns( - units=[unit.split(' ')[1] if ' ' in unit else unit for unit in measurement_and_unit], - unit_counts=df[column_overview.name].value_counts(dropna=False).to_dict(), - with_measurement=any([has_number_and_unit(unit) for unit in measurement_and_unit]) - ) - return False - -def match_units(entries: List[str], regex: re.Pattern) -> List[str]: - units: List[str] = [] - for entrie in entries: - measurement_and_unit = regex.fullmatch(entrie) - if measurement_and_unit: - units.append(measurement_and_unit.group(0)) - return units - -def has_number_and_unit(value: str) -> bool: - special_units = {'1/s', '1/m', '1/M', '1/h', '1/min'} # Add new ones, if needed - if value in special_units: - return False - pattern = re.compile(r'^-?\d+\.?\d*\s*[a-zA-Z°/%]+$') - return bool(pattern.match(value) and value not in special_units) diff --git a/src/qc_eda/biological/taxonomy.py b/src/qc_eda/biological/taxonomy.py deleted file mode 100644 index 0621339..0000000 --- a/src/qc_eda/biological/taxonomy.py +++ /dev/null @@ -1,86 +0,0 @@ -from dataclasses import dataclass -import pandas as pd - -@dataclass -class TaxonomyFlags: - name: str - is_taxonomy: bool - taxid: set | str | None - taxonomy: list | str | None - - -def taxonomy_flags(df, col, tax_df) -> TaxonomyFlags: - - if df[col].dtype in ['int64', 'float64'] or pd.api.types.is_numeric_dtype(df[col]): - taxid_result = is_taxid(df[col], tax_df) - if taxid_result is not None: - return TaxonomyFlags( - name=col, - is_taxonomy=True, - taxid=taxid_result, - taxonomy=None - ) - else: - taxonomy_result = is_taxonomy(df[col], tax_df) - if taxonomy_result is not None: - return TaxonomyFlags( - name=col, - is_taxonomy=True, - taxid=None, - taxonomy=taxonomy_result - ) - - return TaxonomyFlags( - name=col, - is_taxonomy=False, - taxid=None, - taxonomy=None - ) - - -def is_taxid(col: pd.Series, tax_df: pd.DataFrame, threshold: float = 0.9) -> set | str | None: - valid_tax_ids = set(tax_df['tax_id']) - excluded_cols = ["length", "start", "end"] - - if col.name and str(col.name).lower() in excluded_cols: - return None - - tmp_series = pd.to_numeric(col, errors='coerce') - is_numeric_candidate = tmp_series.notna().sum() / len(col) > threshold - - if is_numeric_candidate: - is_valid = tmp_series.isin(valid_tax_ids) - validity_rate = is_valid.sum() / len(col) - - if validity_rate > threshold: - invalid_mask = ~is_valid & tmp_series.notna() - invalid_ids = set(col.loc[invalid_mask].tolist()) - - if invalid_ids: - return invalid_ids - else: - return "all tax IDs valid" - - return None - - -def is_taxonomy(col: pd.Series, tax_df: pd.DataFrame, threshold: float = 0.8) -> list | str | None: - valid_names = set(tax_df['name_txt']) - is_valid = col.isin(valid_names) - validity_rate = is_valid.sum() / len(col) - if validity_rate < threshold: - cleaned_names = col.astype(str).str.extract(r'^([^(]+)')[0].str.strip() - is_valid_cleaned = cleaned_names.isin(valid_names) - validity_rate_cleaned = is_valid_cleaned.sum() / len(col) - - if validity_rate_cleaned > validity_rate: - is_valid = is_valid_cleaned - validity_rate = validity_rate_cleaned - - if validity_rate > threshold: - invalid_names_list = col.loc[~is_valid & col.notnull()].unique().tolist() - if invalid_names_list: - return sorted(invalid_names_list) - else: - return "Valid" - return None diff --git a/src/quality_assessment/__init__.py b/src/quality_assessment/__init__.py new file mode 100644 index 0000000..350b1bd --- /dev/null +++ b/src/quality_assessment/__init__.py @@ -0,0 +1 @@ +from .quality_assessment import quality_assessment, print_quality_report \ No newline at end of file diff --git a/src/quality_assessment/biological_quality.py b/src/quality_assessment/biological_quality.py new file mode 100644 index 0000000..67532db --- /dev/null +++ b/src/quality_assessment/biological_quality.py @@ -0,0 +1,104 @@ +from models import QualityCheck +from .utils import _worst + +AMBIGUOUS_SEQ_WARN = 5.0 +AMBIGUOUS_SEQ_FAIL = 15.0 +STOP_CODON_WARN = 1.0 +RC_REDUNDANCY_WARN = 10.0 +RC_REDUNDANCY_FAIL = 30.0 + + +def _check_sequence_validity(column_overviews) -> QualityCheck: + seq_cols = [c for c in column_overviews if c.sequence in ("dna", "protein")] + if not seq_cols: + return QualityCheck(name="Sequence Validity", status="pass", + message="No sequence columns present", detail_link=None) + statuses, notes = [], [] + for c in seq_cols: + n_invalid = len(c.invalid_seqs) if c.invalid_seqs else 0 + if n_invalid > 0: + statuses.append("warn") + notes.append(f"{c.name}: {n_invalid} invalid") + data = getattr(c, "dna_rna_data", None) or getattr(c, "protein_data", None) + if data is None: + continue + amb = getattr(data, "ambiguous_base_ratio", None) or getattr(data, "ambiguous_residue_ratio", None) + if amb is not None: + if amb.mean > AMBIGUOUS_SEQ_FAIL: + statuses.append("fail") + notes.append(f"{c.name}: {amb.mean:.1f}% mean ambiguous") + elif amb.mean > AMBIGUOUS_SEQ_WARN: + statuses.append("warn") + notes.append(f"{c.name}: {amb.mean:.1f}% mean ambiguous") + stop = getattr(data, "stop_codon_ratio", None) + if stop is not None and stop > STOP_CODON_WARN: + statuses.append("warn") + notes.append(f"{c.name}: {stop:.1f}% with stop codons") + status = _worst(statuses) if statuses else "pass" + return QualityCheck(name="Sequence Validity", status=status, + message="; ".join(notes) if notes else "All sequences valid", + detail_link="#columns") + + +def _check_sequence_redundancy(column_overviews) -> QualityCheck: + dna_cols = [c for c in column_overviews if c.sequence == "dna"] + if not dna_cols: + return QualityCheck(name="Sequence Redundancy (Reverse Complement)", status="pass", + message="No DNA/RNA columns present", detail_link=None) + statuses, notes = [], [] + for c in dna_cols: + data = getattr(c, "dna_rna_data", None) + if data is None: + continue + rc = getattr(data, "reverse_complement_ratio", None) + if rc is None: + continue + if rc > RC_REDUNDANCY_FAIL: + statuses.append("fail") + notes.append(f"{c.name}: {rc:.1f}% RC duplicates") + elif rc > RC_REDUNDANCY_WARN: + statuses.append("warn") + notes.append(f"{c.name}: {rc:.1f}% RC duplicates") + status = _worst(statuses) if statuses else "pass" + return QualityCheck(name="Sequence Redundancy (Reverse Complement)", status=status, + message="; ".join(notes) if notes else "No significant reverse-complement redundancy", + detail_link="#columns") + + +def _check_taxonomy_validity(column_overviews) -> QualityCheck: + tax_cols = [c for c in column_overviews if c.taxonomy is not None] + if not tax_cols: + return QualityCheck(name="Taxonomy Validity", status="pass", + message="No taxonomy columns present", detail_link=None) + statuses, notes = [], [] + for c in tax_cols: + n_invalid = len(c.taxonomy.invalid_names) if c.taxonomy.invalid_names else 0 + if n_invalid > 0: + statuses.append("warn") + notes.append(f"{c.name}: {n_invalid} invalid") + + status = _worst(statuses) if statuses else "pass" + return QualityCheck(name="Taxonomy Validity", status=status, + message="; ".join(notes) if notes else "Taxonomy valid", + detail_link="#columns") + + +def _check_unit_validity(column_overviews) -> QualityCheck: + unit_cols = [c for c in column_overviews if c.measurement_data is not None] + if not unit_cols: + return QualityCheck(name="Unit Validity", status="pass", + message="No unit columns present", detail_link=None) + statuses, notes = [], [] + for c in unit_cols: + ucount = len(c.measurement_data.unit_counts) if c.measurement_data.unit_counts else 0 + if ucount > 1: + statuses.append("warn") + notes.append(f"{c.name}: {ucount} units found.") + if c.measurement_data.with_measurement: + statuses.append("warn") + notes.append(f"{c.name}: measurement with unit found.") + + status = _worst(statuses) if statuses else "pass" + return QualityCheck(name="Unit Validity", status=status, + message="; ".join(notes) if notes else "Unit valid", + detail_link="#columns") diff --git a/src/quality_assessment/column_quality.py b/src/quality_assessment/column_quality.py new file mode 100644 index 0000000..1300a73 --- /dev/null +++ b/src/quality_assessment/column_quality.py @@ -0,0 +1,131 @@ +from quality_assessment.utils import _worst, _rate +from models import QualityCheck + +MISSING_COL_WARN = 0.30 +MISSING_COL_FAIL = 0.50 + +QUASI_CONSTANT_WARN = 0.95 +HIGH_CARDINALITY_WARN = 0.90 + +SKEW_WARN = 1.0 +SKEW_FAIL = 2.0 + +OUTLIER_WARN = 0.05 +OUTLIER_FAIL = 0.15 + +def _missing_rate(col): + mp = col.missing_per or 0.0 + return mp / 100.0 if mp > 1 else mp + + +def _check_missing(column_overviews) -> list[QualityCheck]: + flagged, empty, worst = [], [], "pass" + for c in column_overviews: + rate = _missing_rate(c) + if rate >= 1.0: + empty.append(c.name) + elif rate > MISSING_COL_FAIL: + flagged.append((c.name, rate)); worst = "fail" + elif rate > MISSING_COL_WARN: + flagged.append((c.name, rate)) + worst = "fail" if worst == "fail" else "warn" + + checks = [] + if flagged: + flagged.sort(key=lambda x: x[1], reverse=True) + listed = ", ".join(f"{n} ({r*100:.0f}%)" for n, r in flagged) + checks.append(QualityCheck("Missing Values", worst, + f"Columns above 30% missing: {listed}", "#multivariate")) + if empty: + checks.append(QualityCheck("Empty Column", "fail", + f"Empty column(s): {', '.join(empty)}", "#multivariate")) + if not checks: + checks.append(QualityCheck("Missing Values", "pass", + "No column above 30% missing", "#multivariate")) + return checks + +def _check_variance(column_overviews, categorical_overviews) -> QualityCheck: + constant_cols = [c.name for c in column_overviews if c.constant] + quasi = [c.name for c in categorical_overviews + if c is not None and c.top_1_coverage is not None and c.top_1_coverage >= QUASI_CONSTANT_WARN] + if constant_cols: + return QualityCheck(name="Feature Variance", status="fail", + message=f"Constant column(s): {', '.join(constant_cols)}", detail_link="#columns") + if quasi: + return QualityCheck(name="Feature Variance", status="warn", + message=f"Quasi-constant column(s) (top category ≥{QUASI_CONSTANT_WARN:.0%}): {', '.join(quasi)}", + detail_link="#columns") + return QualityCheck(name="Feature Variance", status="pass", + message="No constant or quasi-constant columns", detail_link="#columns") + + +def _check_high_cardinality(column_overviews) -> QualityCheck: + flagged = [c.name for c in column_overviews + if c.cardinality_dimension_ratio is not None and c.cardinality_dimension_ratio >= HIGH_CARDINALITY_WARN] + status = "warn" if flagged else "pass" + return QualityCheck( + name="High Cardinality", status=status, + message=(f"High-cardinality column(s) (≥{HIGH_CARDINALITY_WARN:.0%} unique): {', '.join(flagged)}" + if flagged else "No high-cardinality columns"), + detail_link="#columns", + ) + + +def _check_mixed_types(column_overviews) -> QualityCheck: + flagged = [c.name for c in column_overviews if c.mixed_types] + status = "fail" if flagged else "pass" + return QualityCheck( + name="Mixed Types", status=status, + message=(f"Mixed-type column(s): {', '.join(flagged)}" if flagged else "No mixed-type columns"), + detail_link="#columns", + ) + + +def _check_infinity(numeric_overviews) -> QualityCheck: + total_inf = sum(int(n.infinity) for n in numeric_overviews if n.infinity) + status = "fail" if total_inf > 0 else "pass" + return QualityCheck(name="Infinity Check", status=status, + message=(f"{total_inf} infinite value(s)" if total_inf else "No infinite values"), + detail_link="#columns") + + +def _check_skewness(numeric_overviews) -> QualityCheck: + statuses, notes = [], [] + for n in numeric_overviews: + s = abs(n.skewness) if n.skewness == n.skewness else 0.0 # guard NaN + if s > SKEW_FAIL: + statuses.append("fail") + notes.append(f"{n.name} ({n.skewness:.1f})") + elif s > SKEW_WARN: + statuses.append("warn") + notes.append(f"{n.name} ({n.skewness:.1f})") + status = _worst(statuses) if statuses else "pass" + return QualityCheck( + name="Skewness", status=status, + message=(f"Strongly skewed (transform recommended): {', '.join(notes)}" + if notes else "No strongly skewed numeric columns"), + detail_link="#columns", + ) + + +def _check_outliers(numeric_overviews, n_rows) -> QualityCheck: + statuses, notes = [], [] + for n in numeric_overviews: + o = n.outliers + if o is None: + continue + flagged = (o.n_lower_iqr or 0) + (o.n_upper_iqr or 0) # adjusted IQR + rate = _rate(flagged, n_rows) + if rate > OUTLIER_FAIL: + statuses.append("fail") + notes.append(f"{n.name} ({rate * 100:.1f}%)") + elif rate > OUTLIER_WARN: + statuses.append("warn") + notes.append(f"{n.name} ({rate * 100:.1f}%)") + status = _worst(statuses) if statuses else "pass" + return QualityCheck( + name="Outliers (Adjusted IQR)", status=status, + message=(f"High outlier share: {', '.join(notes)}" if notes + else "No columns with excessive outliers"), + detail_link="#columns", + ) diff --git a/src/quality_assessment/dataset_quality.py b/src/quality_assessment/dataset_quality.py new file mode 100644 index 0000000..a2d826e --- /dev/null +++ b/src/quality_assessment/dataset_quality.py @@ -0,0 +1,66 @@ +from models import QualityCheck +from quality_assessment.utils import _rate + +MIN_ROWS_PER_FEATURE_WARN = 10 +MIN_ROWS_PER_FEATURE_FAIL = 3 + +DUP_ROW_WARN = 0.05 +DUP_ROW_FAIL = 0.20 +EMPTY_ROW_WARN = 0.01 + +MISSING_COL_WARN = 0.30 +MISSING_COL_FAIL = 0.50 + + +def _check_sample_size(general) -> QualityCheck: + ratio = general.ratio + if ratio < MIN_ROWS_PER_FEATURE_FAIL: + status = "fail" + elif ratio < MIN_ROWS_PER_FEATURE_WARN: + status = "warn" + else: + status = "pass" + return QualityCheck( + name="Sample Size vs Feature Count", status=status, + message=f"{ratio:.1f} rows per feature " + f"({general.rows} rows / {general.cols} columns)", + detail_link="#overview", + ) + + +def _check_duplicate_rows(general) -> QualityCheck: + dup = general.dup_row + rate = _rate(dup, general.rows) + if rate > DUP_ROW_FAIL: + status = "fail" + elif rate > DUP_ROW_WARN: + status = "warn" + else: + status = "pass" + unique_rows = general.rows - dup + return QualityCheck( + name="Duplicate Rows", status=status, + message=(f"{dup} duplicate rows ({rate * 100:.1f}%); " + f"{unique_rows} unique of {general.rows}" if dup else "No duplicate rows"), + detail_link="#overview", + ) + + +def _check_duplicate_columns(general) -> QualityCheck: + dup = general.dup_col + status = "warn" if dup > 0 else "pass" + return QualityCheck( + name="Duplicate Columns", status=status, + message=(f"{dup} duplicate column(s)" if dup else "No duplicate columns"), + detail_link="#overview", + ) + + +def _check_empty_rows(general) -> QualityCheck: + empty = general.empty_rows + status = "warn" if empty > 0 else "pass" + return QualityCheck( + name="Empty Rows", status=status, + message=(f"{empty} completely empty row(s)" if empty else "No empty rows"), + detail_link="#overview", + ) diff --git a/src/quality_assessment/quality_assessment.py b/src/quality_assessment/quality_assessment.py new file mode 100644 index 0000000..26180ec --- /dev/null +++ b/src/quality_assessment/quality_assessment.py @@ -0,0 +1,81 @@ +from termcolor import colored + +from models.quality import QualityCategory, QualityAssessment +from quality_assessment.biological_quality import _check_sequence_validity, _check_sequence_redundancy, \ + _check_taxonomy_validity, _check_unit_validity +from quality_assessment.column_quality import _check_missing, _check_variance, _check_high_cardinality, \ + _check_mixed_types, _check_infinity, _check_skewness, _check_outliers +from quality_assessment.dataset_quality import _check_sample_size, _check_duplicate_rows, _check_duplicate_columns, \ + _check_empty_rows +from quality_assessment.relationships import _check_leakage, _check_multicollinearity +from quality_assessment.utils import _worst + + +def quality_assessment(general, column_overviews, numeric_overviews, categorical_overviews, multivariate) -> QualityAssessment: + numeric_overviews = [n for n in numeric_overviews if n is not None] + categorical_overviews = [c for c in categorical_overviews if c is not None] + + structure = [ + _check_sample_size(general), + _check_duplicate_rows(general), + _check_duplicate_columns(general), + _check_empty_rows(general), + ] + column_quality = [ + *_check_missing(column_overviews), + _check_variance(column_overviews, categorical_overviews), + _check_high_cardinality(column_overviews), + _check_mixed_types(column_overviews), + _check_infinity(numeric_overviews), + _check_skewness(numeric_overviews), + _check_outliers(numeric_overviews, general.rows), + ] + biological = [ + _check_sequence_validity(column_overviews), + _check_sequence_redundancy(column_overviews), + _check_taxonomy_validity(column_overviews), + _check_unit_validity(column_overviews), + ] + relationships = [ + _check_leakage(multivariate), + _check_multicollinearity(multivariate), + ] + + categories = [ + QualityCategory("Dataset Structure", _worst([c.status for c in structure]), structure), + QualityCategory("Column Quality", _worst([c.status for c in column_quality]), column_quality), + QualityCategory("Biological", _worst([c.status for c in biological]), biological), + QualityCategory("Relationships", _worst([c.status for c in relationships]), relationships), + ] + + all_checks = [c for cat in categories for c in cat.checks] + passed = sum(1 for c in all_checks if c.status == "pass") + warnings = sum(1 for c in all_checks if c.status == "warn") + failed = sum(1 for c in all_checks if c.status == "fail") + + overall = "not_ready" if failed else ("caution" if warnings else "ready") + + return QualityAssessment( + categories=categories, passed=passed, warnings=warnings, + failed=failed, total=len(all_checks), overall=overall, + ) + + +def print_quality_report(quality: QualityAssessment) -> None: + symbols = {"pass": colored("[PASS]", 'green'), "warn": colored("[WARN]", 'yellow'), "fail": colored("[FAIL]", 'red')} + overall_label = {"ready": colored("NO QUALITY ISSUES", 'green'), "caution": colored("MINOR QUALITY ISSUES", 'yellow'), "not_ready": colored("MAJOR QUALITY ISSUES", 'red')} + + print("\n" + "=" * 64) + print("QUALITY ASSESSMENT".center(64)) + print("=" * 64) + print(f" Overall: {overall_label.get(quality.overall, quality.overall)}") + print(f" {quality.passed} passed | {quality.warnings} warnings | " + f"{quality.failed} failed (of {quality.total})") + + for cat in quality.categories: + print("-" * 64) + print(f" {symbols.get(cat.status)} {cat.name.upper()}") + for c in cat.checks: + print(f" {symbols.get(c.status)} {c.name}") + print(f" {c.message}") + print("=" * 64) diff --git a/src/quality_assessment/relationships.py b/src/quality_assessment/relationships.py new file mode 100644 index 0000000..0a6fae2 --- /dev/null +++ b/src/quality_assessment/relationships.py @@ -0,0 +1,37 @@ +from models import QualityCheck + +LEAKAGE_WARN = 0.90 +LEAKAGE_FAIL = 0.98 +MULTICOLLINEAR_WARN = 0.90 + + +def _check_leakage(multivariate) -> QualityCheck: + ftc = getattr(multivariate, "feature_target_correlation", None) + if not ftc: + return QualityCheck(name="Target Leakage", status="pass", + message="No target set; leakage check skipped", detail_link=None) + suspects = [(f, info["value"]) for f, info in ftc.items() if info["value"] >= LEAKAGE_WARN] + if not suspects: + return QualityCheck(name="Target Leakage", status="pass", + message="No feature strongly associated with the target", + detail_link="#multivariate") + suspects.sort(key=lambda x: x[1], reverse=True) + status = "fail" if suspects[0][1] >= LEAKAGE_FAIL else "warn" + listed = ", ".join(f"{f} ({v:.2f})" for f, v in suspects) + return QualityCheck(name="Target Leakage", status=status, + message=f"Suspiciously high feature-target association: {listed}", + detail_link="#multivariate") + + +def _check_multicollinearity(multivariate) -> QualityCheck: + pairs = getattr(multivariate, "top_associations", None) + if not pairs: + return QualityCheck(name="Multicollinearity", status="pass", + message="No strongly associated feature pairs", detail_link="#multivariate") + high = [p for p in pairs if p["value"] >= MULTICOLLINEAR_WARN] + if not high: + return QualityCheck(name="Multicollinearity", status="pass", + message="No feature pair above 0.90 association", detail_link="#multivariate") + listed = ", ".join(f"{p['var1']}↔{p['var2']} ({p['value']:.2f})" for p in high) + return QualityCheck(name="Multicollinearity", status="warn", + message=f"Highly associated feature pairs: {listed}", detail_link="#multivariate") diff --git a/src/quality_assessment/utils.py b/src/quality_assessment/utils.py new file mode 100644 index 0000000..ede1ef4 --- /dev/null +++ b/src/quality_assessment/utils.py @@ -0,0 +1,10 @@ +def _worst(statuses: list) -> str: + if "fail" in statuses: + return "fail" + if "warn" in statuses: + return "warn" + return "pass" + + +def _rate(value, total): + return value / total if total else 0.0 diff --git a/src/static/css/cb.css b/src/static/css/cb.css index d8022ff..fec3910 100644 --- a/src/static/css/cb.css +++ b/src/static/css/cb.css @@ -5,11 +5,29 @@ * Licensed under MIT (https://github.com/twbs/bootstrap/blob/main/LICENSE) */ +/* ── Shared component styles ─────────────────────────────────────── */ #dup_table tbody td:first-child { font-weight: bold; } +/* ═══════════════════════════════════════════════════════════════════ + Theme colour palettes + Apply data-cb-theme="internal" (default) or "colorblind" to . + + Semantic aliases translate palette positions into component roles so + every component rule below is written once and works in both themes. + + Two alias families are provided: + --cb-action* → background / border use (buttons, active states) + --cb-action-text* → text / foreground use (links, labels, tab text) + + These must be separate because WCAG contrast requirements differ: + a colour used as a background only needs to contrast with its text, + while a colour used as text on white must reach ≥ 4.5:1 on its own. + ═══════════════════════════════════════════════════════════════════ */ + [data-cb-theme=internal] { + /* Raw palette */ --cb-blue-1: #00447B; --cb-blue-2: #0F65A0; --cb-blue-3: #4082C0; @@ -22,9 +40,24 @@ --cb-pink-2: #A83665; --cb-pink-3: #994564; --cb-pink-4: #D27897; + + /* Background / border aliases */ + --cb-action: var(--cb-pink-1); /* primary CTA bg #932263 */ + --cb-action-hover: var(--cb-pink-2); /* hover / active bg #A83665 */ + --cb-action-subtle: var(--cb-pink-4); /* border / focus ring #D27897 */ + --cb-surface: var(--cb-grey-1); /* dark navbar / footer #3D4756 */ + + /* Text-on-filled-action aliases */ + --cb-on-action: #fff; /* text on filled button — 7.9:1 */ + --cb-on-action-hover: #fff; + + /* Text-on-white aliases (internal theme: same as action bg) */ + --cb-action-text: var(--cb-pink-1); /* links, tabs, labels 7.9:1 ✓ */ + --cb-action-hover-text: var(--cb-pink-2); /* hover text 6.2:1 ✓ */ } [data-cb-theme=colorblind] { + /* Raw palette */ --cb-blue-1: #616D78; --cb-blue-2: #5B666A; --cb-blue-3: #6C7385; @@ -37,8 +70,59 @@ --cb-pink-2: #5EC8B6; --cb-pink-3: #5A969A; --cb-pink-4: #378979; + + /* Background / border aliases */ + --cb-action: var(--cb-pink-4); /* primary CTA bg #378979 */ + --cb-action-hover: var(--cb-pink-3); /* hover / active bg #5A969A */ + --cb-action-subtle: var(--cb-pink-2); /* border / focus ring #5EC8B6 */ + --cb-surface: var(--cb-grey-1); /* light navbar / footer #7DB3C6 */ + + /* + * Text-on-filled-action: #334A52 (L=0.062) has only 2.2:1 on #378979. + * Pure black (L=0) gives 5.0:1 on #378979 and 6.3:1 on #5A969A. ✓ + */ + --cb-on-action: #000; + --cb-on-action-hover: #000; + + /* + * Text-on-white: #378979 gives only 4.2:1 — below WCAG AA (4.5:1). + * --cb-grey-2 (#334A52) gives 9.4:1 and shares the teal character. ✓ + */ + --cb-action-text: var(--cb-grey-2); /* links, labels on white 9.4:1 ✓ */ + --cb-action-hover-text: var(--cb-grey-2); /* hover stays dark; underline signals state */ + + /* + * Bootstrap semantic colour overrides — Okabe-Ito palette. + * These drive .bg-danger / .bg-success / .bg-warning badges, + * .table-danger / -success / -warning row tints, and alerts. + * + * Red → Vermillion #D55E00 (safe for protanopia + deuteranopia) + * Green → Blue #0072B2 (clearly distinct from vermillion) + * Yellow → Orange-yellow #E69F00 (safe for all common types) + */ + --bs-danger: #D55E00; + --bs-danger-rgb: 213, 94, 0; + --bs-danger-bg-subtle: #FAE0CC; + --bs-danger-border-subtle: #ECA066; + --bs-danger-text-emphasis: #7F3700; + + --bs-success: #0072B2; + --bs-success-rgb: 0, 114, 178; + --bs-success-bg-subtle: #CCE5F5; + --bs-success-border-subtle: #80BEE0; + --bs-success-text-emphasis: #00416A; + + --bs-warning: #E69F00; + --bs-warning-rgb: 230, 159, 0; + --bs-warning-bg-subtle: #FBF0CC; + --bs-warning-border-subtle: #F3CF66; + --bs-warning-text-emphasis: #7F5700; } +/* ═══════════════════════════════════════════════════════════════════ + Bootstrap utility overrides + ═══════════════════════════════════════════════════════════════════ */ + .text-bg-dark { color: #fff !important; background-color: RGBA(var(--bs-dark-rgb), var(--bs-bg-opacity, 1)) !important; @@ -46,60 +130,134 @@ .bg-dark { --bs-bg-opacity: 1; - background-color: #3D4756 !important; + background-color: var(--cb-surface) !important; } .bg-secondary { --bs-bg-opacity: 1; - background-color: #3D4756 !important; + background-color: var(--cb-surface) !important; } +/* + * In colorblind mode the surface colour (#7DB3C6) is a light blue. + * Bootstrap's .text-white uses !important, so we must counter it here + * for any white text that appears on top of a surface-coloured container. + * Also invert the footer GitHub icon filter so it stays visible. + */ +/* + * #334A52 (cb-grey-2) gives only 4.1:1 on the CB surface (#7DB3C6), + * which is below AA for small body text. Use near-black instead (8.0:1 ✓). + */ +[data-cb-theme=colorblind] .bg-dark .text-white, +[data-cb-theme=colorblind] .bg-secondary .text-white { + color: #1a1a1a !important; +} + +[data-cb-theme=colorblind] #footer img { + filter: brightness(0) !important; /* keeps the SVG logo legible on light bg */ +} + +/* ── Primary filled button ────────────────────────────────────────── */ .btn-cb-primary { - --bs-btn-color: #fff; - --bs-btn-bg: #932263; - --bs-btn-border-color: #932263; - --bs-btn-hover-color: #fff; - --bs-btn-hover-bg: #A83665; - --bs-btn-hover-border-color: #D27897; - --bs-btn-focus-shadow-rgb: 49, 132, 253; - --bs-btn-active-color: #fff; - --bs-btn-active-bg: #A83665; - --bs-btn-active-border-color: #D27897; - --bs-btn-active-shadow: inset 0 3px 5px rgba(0, 0, 0, 0.125); - --bs-btn-disabled-color: #fff; - --bs-btn-disabled-bg: #932263; - --bs-btn-disabled-border-color: #932263; + --bs-btn-color: var(--cb-on-action); + --bs-btn-bg: var(--cb-action); + --bs-btn-border-color: var(--cb-action); + --bs-btn-hover-color: var(--cb-on-action-hover); + --bs-btn-hover-bg: var(--cb-action-hover); + --bs-btn-hover-border-color: var(--cb-action-subtle); + --bs-btn-focus-shadow-rgb: 49, 132, 253; + --bs-btn-active-color: var(--cb-on-action-hover); + --bs-btn-active-bg: var(--cb-action-hover); + --bs-btn-active-border-color: var(--cb-action-subtle); + --bs-btn-active-shadow: inset 0 3px 5px rgba(0, 0, 0, 0.125); + --bs-btn-disabled-color: var(--cb-on-action); + --bs-btn-disabled-bg: var(--cb-action); + --bs-btn-disabled-border-color: var(--cb-action); } +/* ── Primary outline button ──────────────────────────────────────── */ .btn-cb-outline-primary { - --bs-btn-color: #932263; - --bs-btn-border-color: #932263; - --bs-btn-hover-color: #fff; - --bs-btn-hover-bg: #A83665; - --bs-btn-hover-border-color: #D27897; - --bs-btn-focus-shadow-rgb: 13, 110, 253; - --bs-btn-active-color: #fff; - --bs-btn-active-bg: #A83665; - --bs-btn-active-border-color: #D27897; - --bs-btn-active-shadow: inset 0 3px 5px rgba(0, 0, 0, 0.125); - --bs-btn-disabled-color: #932263; - --bs-btn-disabled-bg: transparent; - --bs-btn-disabled-border-color: #932263; + --bs-btn-color: var(--cb-action-text); + --bs-btn-border-color: var(--cb-action); + --bs-btn-hover-color: var(--cb-on-action-hover); + --bs-btn-hover-bg: var(--cb-action-hover); + --bs-btn-hover-border-color: var(--cb-action-subtle); + --bs-btn-focus-shadow-rgb: 13, 110, 253; + --bs-btn-active-color: var(--cb-on-action-hover); + --bs-btn-active-bg: var(--cb-action-hover); + --bs-btn-active-border-color: var(--cb-action-subtle); + --bs-btn-active-shadow: inset 0 3px 5px rgba(0, 0, 0, 0.125); + --bs-btn-disabled-color: var(--cb-action-text); + --bs-btn-disabled-bg: transparent; + --bs-btn-disabled-border-color: var(--cb-action); --bs-gradient: none; } +/* ── Navbar ──────────────────────────────────────────────────────── */ .navbar-cb-dark, .navbar[data-cb-theme=internal] { - --bs-navbar-color: rgba(255, 255, 255, 0.55); - --bs-navbar-hover-color: rgba(255, 255, 255, 0.75); - --bs-navbar-disabled-color: rgba(255, 255, 255, 0.25); - --bs-navbar-active-color: #fff; - --bs-navbar-brand-color: #fff; - --bs-navbar-brand-hover-color: #fff; + --bs-navbar-color: rgba(255, 255, 255, 0.55); + --bs-navbar-hover-color: rgba(255, 255, 255, 0.75); + --bs-navbar-disabled-color: rgba(255, 255, 255, 0.25); + --bs-navbar-active-color: #fff; + --bs-navbar-brand-color: #fff; + --bs-navbar-brand-hover-color: #fff; --bs-navbar-toggler-border-color: rgba(255, 255, 255, 0.1); --bs-navbar-toggler-icon-bg: url("data:image/svg+xml,%3csvg xmlns='http://www.w3.org/2000/svg' viewBox='0 0 30 30'%3e%3cpath stroke='rgba%28255, 255, 255, 0.55%29' stroke-linecap='round' stroke-miterlimit='10' stroke-width='2' d='M4 7h22M4 15h22M4 23h22'/%3e%3c/svg%3e"); } +/* In colorblind mode the surface (#7DB3C6) is light — switch to dark text. */ +[data-cb-theme=colorblind] .navbar-cb-dark { + --bs-navbar-color: rgba(0, 0, 0, 0.65); + --bs-navbar-hover-color: rgba(0, 0, 0, 0.85); + --bs-navbar-disabled-color: rgba(0, 0, 0, 0.30); + --bs-navbar-active-color: #000; + --bs-navbar-brand-color: #000; + --bs-navbar-brand-hover-color: #000; + --bs-navbar-toggler-border-color: rgba(0, 0, 0, 0.2); + --bs-navbar-toggler-icon-bg: url("data:image/svg+xml,%3csvg xmlns='http://www.w3.org/2000/svg' viewBox='0 0 30 30'%3e%3cpath stroke='rgba%280, 0, 0, 0.65%29' stroke-linecap='round' stroke-miterlimit='10' stroke-width='2' d='M4 7h22M4 15h22M4 23h22'/%3e%3c/svg%3e"); +} + +/* + * ── Navbar toggle button (.btn-cb-navbar) ──────────────────────── + * Purpose-built for the navbar context so it is legible regardless + * of whether the palette has been switched: + * Internal mode (dark surface #3D4756) → white outline / white text + * Colorblind mode (light surface #7DB3C6) → dark outline / dark text + */ +.btn-cb-navbar { + --bs-btn-color: rgba(255, 255, 255, 0.85); + --bs-btn-bg: transparent; + --bs-btn-border-color: rgba(255, 255, 255, 0.55); + --bs-btn-hover-color: #fff; + --bs-btn-hover-bg: rgba(255, 255, 255, 0.15); + --bs-btn-hover-border-color: rgba(255, 255, 255, 0.85); + --bs-btn-focus-shadow-rgb: 255, 255, 255; + --bs-btn-active-color: #fff; + --bs-btn-active-bg: rgba(255, 255, 255, 0.20); + --bs-btn-active-border-color: #fff; + --bs-btn-active-shadow: inset 0 3px 5px rgba(0, 0, 0, 0.125); + --bs-btn-disabled-color: rgba(255, 255, 255, 0.35); + --bs-btn-disabled-bg: transparent; + --bs-btn-disabled-border-color: rgba(255, 255, 255, 0.25); + --bs-gradient: none; +} + +[data-cb-theme=colorblind] .btn-cb-navbar { + --bs-btn-color: rgba(0, 0, 0, 0.75); + --bs-btn-bg: transparent; + --bs-btn-border-color: rgba(0, 0, 0, 0.45); + --bs-btn-hover-color: #000; + --bs-btn-hover-bg: rgba(0, 0, 0, 0.10); + --bs-btn-hover-border-color: rgba(0, 0, 0, 0.75); + --bs-btn-active-color: #000; + --bs-btn-active-bg: rgba(0, 0, 0, 0.15); + --bs-btn-active-border-color: #000; + --bs-btn-disabled-color: rgba(0, 0, 0, 0.30); + --bs-btn-disabled-border-color: rgba(0, 0, 0, 0.20); +} + +/* ── Layout helpers ──────────────────────────────────────────────── */ #page-container { position: relative; min-height: 100vh; @@ -117,86 +275,124 @@ min-height: 5rem; } +/* ── Search input ────────────────────────────────────────────────── */ input.form-control[type="search"] { background-color: #fff; - color: #932263; - border: 1px solid #932263; + color: var(--cb-action-text); + border: 1px solid var(--cb-action); opacity: 0.65; box-shadow: inset 0 3px 5px rgba(0, 0, 0, 0.125); } +/* ── Pagination ──────────────────────────────────────────────────── */ ul.pagination li.page-item .page-link { background-color: #fff; - color: #932263; - border: 1px solid #932263; + color: var(--cb-action-text); + border: 1px solid var(--cb-action); } ul.pagination li.page-item .page-link:hover { - background-color: #A83665; - color: #fff; - border-color: #A83665; + background-color: var(--cb-action-hover); + color: var(--cb-on-action-hover); + border-color: var(--cb-action-hover); } ul.pagination li.page-item.active .page-link { - background-color: #A83665; - border-color: #A83665; - color: #fff; + background-color: var(--cb-action-hover); + border-color: var(--cb-action-hover); + color: var(--cb-on-action-hover); box-shadow: inset 0 3px 5px rgba(0, 0, 0, 0.125); } ul.pagination li.page-item.disabled .page-link { background-color: #fff; - color: #932263; - border-color: #932263; + color: var(--cb-action-text); + border-color: var(--cb-action); opacity: 0.45; cursor: not-allowed; } ul.pagination li.page-item .page-link.ellipsis { background-color: #fff; - color: #932263; - border: 1px solid #932263; + color: var(--cb-action-text); + border: 1px solid var(--cb-action); pointer-events: none; } +/* ── Internal cross-reference links ─────────────────────────────── */ .cb-pink-link { - color: var(--cb-pink-1) !important; + color: var(--cb-action-text) !important; text-decoration: none; } .cb-pink-link:hover { - color: var(--cb-pink-2) !important; + color: var(--cb-action-hover-text) !important; text-decoration: underline; } +/* ── Navbar dropdown ─────────────────────────────────────────────── */ .navbar-cb-dark .dropdown-menu { --bs-dropdown-bg: #fff; - --bs-dropdown-border-color: #932263; + --bs-dropdown-border-color: var(--cb-action); } .navbar-cb-dark .dropdown-item { - color: #932263; + color: var(--cb-action-text); } .navbar-cb-dark .dropdown-item.active, .navbar-cb-dark .dropdown-item:active { - background-color: #932263; - color: #fff; + background-color: var(--cb-action); + color: var(--cb-on-action); } .navbar-cb-dark .dropdown-item:active { - background-color: #932263; - color: #fff; + background-color: var(--cb-action); + color: var(--cb-on-action); box-shadow: inset 0 3px 5px rgba(0, 0, 0, 0.125); } -/* Tab Navigation */ +/* ── Tab navigation ──────────────────────────────────────────────── */ .nav-tabs .nav-link { - color: #932263; + color: var(--cb-action-text); } /* Tab Navigation - Active */ .nav-tabs .nav-link.active, .nav-tabs .nav-item.show .nav-link { - color: #932263; + color: var(--cb-action-text); } + +/* ── Table section header rows ───────────────────────────────────── */ +.table-secondary { + --bs-table-color: var(--cb-action-text); + --bs-table-bg: #ffffff; + --bs-table-border-color: #dee2e6; +} + +.table > :not(caption) > .table-secondary > * { + background-color: #ffffff !important; + box-shadow: none !important; + color: var(--cb-action-text) !important; + padding-top: 1rem !important; +} + +/* ═══════════════════════════════════════════════════════════════════ + Colorblind mode — badge / alert / table-row text contrast fixes. + + Bootstrap hardcodes white text on .bg-danger and .bg-success. + Vermillion (#D55E00) and orange-yellow (#E69F00) require dark text + for WCAG AA. Blue (#0072B2) passes at 5.2:1 with white text. + ═══════════════════════════════════════════════════════════════════ */ + +[data-cb-theme=colorblind] .bg-danger, +[data-cb-theme=colorblind] .badge.bg-danger, +[data-cb-theme=colorblind] .alert-danger { + color: #000 !important; +} + +[data-cb-theme=colorblind] .bg-warning, +[data-cb-theme=colorblind] .badge.bg-warning, +[data-cb-theme=colorblind] .alert-warning { + color: #000 !important; +} \ No newline at end of file diff --git a/src/templates/LandingPage.jinja b/src/templates/LandingPage.jinja index 2354628..2407e18 100644 --- a/src/templates/LandingPage.jinja +++ b/src/templates/LandingPage.jinja @@ -5,7 +5,7 @@ - BioProfileKit + BioProfileKit – Home @@ -14,7 +14,7 @@ - +