diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml new file mode 100644 index 0000000..5ab1019 --- /dev/null +++ b/.github/workflows/lint.yml @@ -0,0 +1,45 @@ +name: Lint + +on: + pull_request: + push: + branches: [main] + +jobs: + lint: + name: Lint + runs-on: ubuntu-latest + + permissions: + contents: read + packages: read + statuses: write + + steps: + - name: Checkout code + uses: actions/checkout@v5 + with: + fetch-depth: 0 + persist-credentials: false + + - name: Lint Code Base + uses: super-linter/super-linter@v8.5.0 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + VALIDATE_ALL_CODEBASE: true + DEFAULT_BRANCH: main + + VALIDATE_SNAKEMAKE_SNAKEFMT: true + VALIDATE_PYTHON_BLACK: true + VALIDATE_YAML_PRETTIER: true + + - name: Set up R + uses: r-lib/actions/setup-r@v2 + + - name: Set up R dependencies + uses: r-lib/actions/setup-r-dependencies@v2 + with: + packages: styler + + - name: Check R style + run: Rscript -e 'changed <- styler::style_dir(".", dry = "on"); if (any(changed$changed)) quit(status = 1)' diff --git a/.gitignore b/.gitignore index 8de48c1..f8e44a9 100644 --- a/.gitignore +++ b/.gitignore @@ -60,4 +60,7 @@ data .test/report.html snakejob_conductor.sh -settings.json \ No newline at end of file +settings.json + +!.github/ +!.prettierrc.yml diff --git a/.prettierrc.yml b/.prettierrc.yml new file mode 100644 index 0000000..250482b --- /dev/null +++ b/.prettierrc.yml @@ -0,0 +1,2 @@ +tabWidth: 2 +useTabs: false diff --git a/CITATION.cff b/CITATION.cff index bafdab9..fdf50dc 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -10,15 +10,15 @@ type: software authors: - given-names: Stephan family-names: Reichl - orcid: 'https://orcid.org/0000-0001-8555-7198' + orcid: "https://orcid.org/0000-0001-8555-7198" affiliation: CeMM Research Center for Molecular Medicine - given-names: Raphael family-names: Bednarsky - orcid: 'https://orcid.org/0009-0005-0404-3424' + orcid: "https://orcid.org/0009-0005-0404-3424" affiliation: CeMM Research Center for Molecular Medicine - given-names: Christoph family-names: Bock - orcid: 'https://orcid.org/0000-0001-6091-3088' + orcid: "https://orcid.org/0000-0001-6091-3088" affiliation: CeMM Research Center for Molecular Medicine identifiers: - type: doi @@ -26,8 +26,8 @@ identifiers: description: >- This DOI represents all versions, and will always resolve to the latest one. -repository-code: 'https://github.com/epigen/unsupervised_analysis' -url: 'https://epigen.github.io/unsupervised_analysis/' +repository-code: "https://github.com/epigen/unsupervised_analysis" +url: "https://epigen.github.io/unsupervised_analysis/" abstract: >- A general purpose Snakemake workflow to perform unsupervised analyses (dimensionality reduction & cluster diff --git a/README.md b/README.md index 8968264..f32f184 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,13 @@ [![MrBiomics](https://img.shields.io/badge/MrBiomics-red)](https://github.com/epigen/MrBiomics/) [![DOI](https://zenodo.org/badge/475465311.svg)](https://zenodo.org/badge/latestdoi/475465311) -[![](https://tokei.rs/b1/github/epigen/unsupervised_analysis?category=code)]() +[![](https://tokei.rs/b1/github/epigen/unsupervised_analysis?category=code)]() [![](https://tokei.rs/b1/github/epigen/unsupervised_analysis?category=files)]() [![GitHub license](https://img.shields.io/github/license/epigen/unsupervised_analysis)](https://github.com/epigen/unsupervised_analysis/blob/master/LICENSE) ![GitHub Release](https://img.shields.io/github/v/release/epigen/unsupervised_analysis) [![Snakemake](https://img.shields.io/badge/Snakemake->=8.20.1-green)](https://snakemake.readthedocs.io/en/stable/) - # Unsupervised Analysis Workflow + A general purpose [Snakemake 8](https://snakemake.readthedocs.io/en/stable/) workflow to perform unsupervised analyses (dimensionality reduction and cluster analysis) on and visualizations of high-dimensional data. > [!NOTE] @@ -21,37 +21,39 @@ A general purpose [Snakemake 8](https://snakemake.readthedocs.io/en/stable/) wor ![Workflow Rulegraph](./workflow/dags/rulegraph.svg) # 🖋️ Authors + - [Stephan Reichl](https://github.com/sreichl) - [Raphael Bednarsky](https://github.com/bednarsky) - [Christoph Bock](https://github.com/chrbock) # 💿 Software -This project wouldn't be possible without the following software and their dependencies -| Software | Reference (DOI) | -| :------------: | :-----------------------------------------------: | -| clusterCrit | https://CRAN.R-project.org/package=clusterCrit | -| clustree | https://doi.org/10.1093/gigascience/giy083 | -| ComplexHeatmap | https://doi.org/10.1093/bioinformatics/btw313 | -| densMAP | https://doi.org/10.1038/s41587-020-00801-7 | -| fastcluster | https://doi.org/10.18637/jss.v053.i09 | -| ggally | https://CRAN.R-project.org/package=GGally | -| ggplot2 | https://ggplot2.tidyverse.org/ | -| ggrepel | https://CRAN.R-project.org/package=ggrepel | -| igraph | https://doi.org/10.5281/zenodo.3630268 | -| leidenalg | https://doi.org/10.5281/zenodo.1469356 | -| pandas | https://doi.org/10.5281/zenodo.3509134 | -| patchwork | https://CRAN.R-project.org/package=patchwork | -| PCA | https://doi.org/10.1080/14786440109462720 | -| plotly express | https://plot.ly | -| pymcdm | https://doi.org/10.1016/j.softx.2023.101368 | -| scikit-learn | http://jmlr.org/papers/v12/pedregosa11a.html | -| scipy | https://doi.org/10.1038/s41592-019-0686-2 | -| Snakemake | https://doi.org/10.12688/f1000research.29032.2 | -| umap-learn | https://doi.org/10.21105/joss.00861 | +This project wouldn't be possible without the following software and their dependencies +| Software | Reference (DOI) | +| :------------: | :--------------------------------------------: | +| clusterCrit | https://CRAN.R-project.org/package=clusterCrit | +| clustree | https://doi.org/10.1093/gigascience/giy083 | +| ComplexHeatmap | https://doi.org/10.1093/bioinformatics/btw313 | +| densMAP | https://doi.org/10.1038/s41587-020-00801-7 | +| fastcluster | https://doi.org/10.18637/jss.v053.i09 | +| ggally | https://CRAN.R-project.org/package=GGally | +| ggplot2 | https://ggplot2.tidyverse.org/ | +| ggrepel | https://CRAN.R-project.org/package=ggrepel | +| igraph | https://doi.org/10.5281/zenodo.3630268 | +| leidenalg | https://doi.org/10.5281/zenodo.1469356 | +| pandas | https://doi.org/10.5281/zenodo.3509134 | +| patchwork | https://CRAN.R-project.org/package=patchwork | +| PCA | https://doi.org/10.1080/14786440109462720 | +| plotly express | https://plot.ly | +| pymcdm | https://doi.org/10.1016/j.softx.2023.101368 | +| scikit-learn | http://jmlr.org/papers/v12/pedregosa11a.html | +| scipy | https://doi.org/10.1038/s41592-019-0686-2 | +| Snakemake | https://doi.org/10.12688/f1000research.29032.2 | +| umap-learn | https://doi.org/10.21105/joss.00861 | # 🔬 Methods + This is a template for the Methods section of a scientific publication and is intended to serve as a starting point. Only retain paragraphs relevant to your analysis. References [ref] to the respective publications are curated in the software table above. Versions (ver) have to be read out from the respective conda environment specifications (`workflow/envs/*.yaml file`) or post-execution in the result directory (`unsupervised_analysis/envs/*.yaml`). Parameters that have to be adapted depending on the data or workflow configurations are denoted in squared brackets e.g. [X]. The outlined analyses were performed using the programming languages R (ver) [ref] and Python (ver) [ref] unless stated otherwise. We applied both linear and non-linear unsupervised analysis methods for dimensionality reduction on normalized data for downstream analyses (e.g., clustering) and to visualize emerging patterns in lower dimensional embeddings. @@ -91,19 +93,21 @@ We performed internal cluster validation using six complementary indices: Silhou **The analysis and visualizations described here were performed using a publicly available Snakemake [ver] (ref) workflow [10.5281/zenodo.8405360](https://doi.org/10.5281/zenodo.8405360).** - # 🚀 Features + The workflow perfroms the following analyses on each dataset provided in the annotation file. A result folder "unsupervised_analysis" is generated containing a folder for each dataset. ## Dimensionality Reduction + > _"High-dimensional spaces are where intuition goes to die and dimensionality reduction becomes the antidote to the curse of dimensionality."_ from Anonymous + - Principal Component Anlaysis (PCA) using [scikit-learn](https://scikit-learn.org) (.pickle and .CSV) - Supporting two configurations that greatly influence performance: [n_components] and [svd_solver] as described in the [documentation](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html). - diagnostics (.PNG): - - variance: scree-plot and cumulative explained variance-plot of all and top 10% principal components - - pairs: sequential pair-wise PCs for up to 10 PCs using scatter- and density-plots colored by [metadata_of_interest] - - loadings: showing the magnitude and direction of the 10 most influential features for each Principal Component combination (Biplot, but without the data) - - loadings lolliplot: showing the magnitude of the 10 most influential features for each Principal Component + - variance: scree-plot and cumulative explained variance-plot of all and top 10% principal components + - pairs: sequential pair-wise PCs for up to 10 PCs using scatter- and density-plots colored by [metadata_of_interest] + - loadings: showing the magnitude and direction of the 10 most influential features for each Principal Component combination (Biplot, but without the data) + - loadings lolliplot: showing the magnitude of the 10 most influential features for each Principal Component - Uniform Manifold Approximation & Projection (UMAP) - k-nearest-neighbor graph (.pickle): generated using the [n_neighbors] parameter together with the provided [metrics]. - fix any pickle load issue by specifying Python version to 3.9 (in case you want to use the graph downstream) @@ -112,69 +116,70 @@ The workflow perfroms the following analyses on each dataset provided in the ann - diagnostics (.PNG): 2D embedding colored by PCA coordinates, vector quantization coordinates, approximated local dimension, neighborhood Jaccard index - connectivity (.PNG): graph/network-connectivity plot with edge-bundling (hammer algorithm variant) - Hierarchically Clustered Heatmap (.PNG) - - Hierarchically clustered heatmaps of scaled data (z-score) with configured distances ([metrics]) and clustering methods ([hclust_methods]). - - Distance matrices of observations and features are precomputed using scipy's dist function. - - Hierarchical clustering is performed by the R implementation of fastcluster. - - The observations can be randomly downsampled by proportion or absolute number ([n_observations]) to reduce computational cost. - - The number of features can be reduced to a proportion or an absolute number of the top variable features ([n_features]) to reduce computational cost. - - All combinations are computed, and annotated with [metadata_of_interest]. + - Hierarchically clustered heatmaps of scaled data (z-score) with configured distances ([metrics]) and clustering methods ([hclust_methods]). + - Distance matrices of observations and features are precomputed using scipy's dist function. + - Hierarchical clustering is performed by the R implementation of fastcluster. + - The observations can be randomly downsampled by proportion or absolute number ([n_observations]) to reduce computational cost. + - The number of features can be reduced to a proportion or an absolute number of the top variable features ([n_features]) to reduce computational cost. + - All combinations are computed, and annotated with [metadata_of_interest]. - Visualization - - 2D metadata and feature plots (.PNG) of the first 2 principal components and all 2D embeddings, respectively. - - interactive 2D and 3D visualizations as self contained HTML files of all projections/embeddings. + - 2D metadata and feature plots (.PNG) of the first 2 principal components and all 2D embeddings, respectively. + - interactive 2D and 3D visualizations as self contained HTML files of all projections/embeddings. - Results directories for each dataset have the following structure: - - "method" (containing all the data as .pickle and/or .CSV files) - - plots (for all visualizations as .PNG files) + - "method" (containing all the data as .pickle and/or .CSV files) + - plots (for all visualizations as .PNG files) ## Cluster Analysis + > _"The validation of clustering structures is the most difficult and frustrating part of cluster analysis. Without a strong effort in this direction, cluster analysis will remain a black art accessible only to those true believers who have experience and great courage."_ from _Algorithms for Clustering Data (1988)_ by Jain & Dubes - Clustering - - Leiden algorithm - - Applied to the UMAP KNN graphs specified by the respective parameters (metric, n_neighbors). - - All algorithm specific parameters are supported: [partition_types], [resolutions], and [n_iterations]. - - Clustification: an ML-based clustering approach that iteratively merges clusters based on misclassification (beta) - 0. User: Specify a clustering method [method]. - 1. Chose the clustering with the most clusters as starting point (i.e., overclustered). - 2. Iterative classification using the cluster labels, to determine if the classifier can distinguish between clusters or if they should be merged. - - Stratified 5-fold CV - - RF with 100 trees (i.e., defaults) - - Retain predicted labels - 3. Merging of clusters. - - Build a normalized confusion matrix using the predicted labels. - - Make it symmetric and upper triangle, resulting in a similarity graph. - - Edge weight ranges from 0 to 1, where 0 means that the classifier was able to distinguish all observations between the two respective clusters. - - Check stopping criterion: if maximum edge weight < 2.5% (i.e., 0.025 – less than 5% of observations are misclassified between any two clusters). - - -> STOP and return current cluster labels - - Merge the two clusters connected by the maximum edge weight. - 4. Back to 2. using the new labels. + - Leiden algorithm + - Applied to the UMAP KNN graphs specified by the respective parameters (metric, n_neighbors). + - All algorithm specific parameters are supported: [partition_types], [resolutions], and [n_iterations]. + - Clustification: an ML-based clustering approach that iteratively merges clusters based on misclassification (beta) 0. User: Specify a clustering method [method]. + 1. Chose the clustering with the most clusters as starting point (i.e., overclustered). + 2. Iterative classification using the cluster labels, to determine if the classifier can distinguish between clusters or if they should be merged. + - Stratified 5-fold CV + - RF with 100 trees (i.e., defaults) + - Retain predicted labels + 3. Merging of clusters. + - Build a normalized confusion matrix using the predicted labels. + - Make it symmetric and upper triangle, resulting in a similarity graph. + - Edge weight ranges from 0 to 1, where 0 means that the classifier was able to distinguish all observations between the two respective clusters. + - Check stopping criterion: if maximum edge weight < 2.5% (i.e., 0.025 – less than 5% of observations are misclassified between any two clusters). + - -> STOP and return current cluster labels + - Merge the two clusters connected by the maximum edge weight. + 4. Back to 2. using the new labels. - Clustree analysis and visualization - - The following clustree specific parameters are supported: [count_filter], [prop_filter], and [layout]. - - default: produces the standard clustree visualization, ordered by number of clusters and annotated. - - custom: extends default by adding [metadata_of_interest] as additional "clusterings". - - metadata and features, specified in the config, are highlighted on top of the clusterings using aggregation functions - - numeric: available aggregation functions: mean, median, max, min - - categorical: available aggregation functions: "pure" or "majority" + - The following clustree specific parameters are supported: [count_filter], [prop_filter], and [layout]. + - default: produces the standard clustree visualization, ordered by number of clusters and annotated. + - custom: extends default by adding [metadata_of_interest] as additional "clusterings". + - metadata and features, specified in the config, are highlighted on top of the clusterings using aggregation functions + - numeric: available aggregation functions: mean, median, max, min + - categorical: available aggregation functions: "pure" or "majority" - Cluster Validation - - External cluster indices are determined comparing all clustering results with all categorical metadata - - all complementary indices from sklearn are used: [AMI](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.adjusted_mutual_info_score.html#sklearn.metrics.adjusted_mutual_info_score), [ARI](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.adjusted_rand_score.html#sklearn.metrics.adjusted_rand_score), [FMI](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.fowlkes_mallows_score.html#sklearn.metrics.fowlkes_mallows_score), [**Homogeneity** and **Completeness** and **V**-Measure](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.homogeneity_completeness_v_measure.html#sklearn.metrics.homogeneity_completeness_v_measure) - - Internal cluster indices are determined for each clustering and [metadata_of_interest] - - 6 complementary indices are used - - 5 from the package [clusterCrit](https://rdrr.io/cran/clusterCrit/man/intCriteria.html): Silhouette, Calinski-Harabasz, C-index, Dunn index, Davis-Bouldin Score. - - 1 weighted Bayesian Information Criterion (BIC) approach, previously described in [Reichl 2018 - Chapter 4.2.2 - Internal Indices](https://repositum.tuwien.at/handle/20.500.12708/3488) - - Due to the computational cost PCA results are used as input and a [sample_proportion] can be configured. - - Caveat: internal cluster indices are linear i.e., using Euclidean distance metrics. - - Multiple-criteria decision-making (MCDM) using TOPSIS for ranking clustering results - - The MCDM method TOPSIS is applied to the internal cluster indices to rank all clustering results (and [metadata_of_interest]) from best to worst. - - This approach has been described in [Reichl 2018 - Chapter 4.3.1 - The Favorite Approach](https://repositum.tuwien.at/handle/20.500.12708/3488) - - Caveat: Silhouette score sometimes generates NA due to a known [bug](https://github.com/cran/clusterCrit/pull/1/commits/b37a5e361d0a12f9d3900089aa03e3947d0d4ef7). Clusterings with NA scores are removed before TOPSIS is applied. + - External cluster indices are determined comparing all clustering results with all categorical metadata + - all complementary indices from sklearn are used: [AMI](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.adjusted_mutual_info_score.html#sklearn.metrics.adjusted_mutual_info_score), [ARI](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.adjusted_rand_score.html#sklearn.metrics.adjusted_rand_score), [FMI](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.fowlkes_mallows_score.html#sklearn.metrics.fowlkes_mallows_score), [**Homogeneity** and **Completeness** and **V**-Measure](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.homogeneity_completeness_v_measure.html#sklearn.metrics.homogeneity_completeness_v_measure) + - Internal cluster indices are determined for each clustering and [metadata_of_interest] + - 6 complementary indices are used + - 5 from the package [clusterCrit](https://rdrr.io/cran/clusterCrit/man/intCriteria.html): Silhouette, Calinski-Harabasz, C-index, Dunn index, Davis-Bouldin Score. + - 1 weighted Bayesian Information Criterion (BIC) approach, previously described in [Reichl 2018 - Chapter 4.2.2 - Internal Indices](https://repositum.tuwien.at/handle/20.500.12708/3488) + - Due to the computational cost PCA results are used as input and a [sample_proportion] can be configured. + - Caveat: internal cluster indices are linear i.e., using Euclidean distance metrics. + - Multiple-criteria decision-making (MCDM) using TOPSIS for ranking clustering results + - The MCDM method TOPSIS is applied to the internal cluster indices to rank all clustering results (and [metadata_of_interest]) from best to worst. + - This approach has been described in [Reichl 2018 - Chapter 4.3.1 - The Favorite Approach](https://repositum.tuwien.at/handle/20.500.12708/3488) + - Caveat: Silhouette score sometimes generates NA due to a known [bug](https://github.com/cran/clusterCrit/pull/1/commits/b37a5e361d0a12f9d3900089aa03e3947d0d4ef7). Clusterings with NA scores are removed before TOPSIS is applied. - Visualization - - all clustering results as 2D and interactive 2D & 3D plots for all available embedings/projections. - - external cluster indices as hierarchically clustered heatmaps, aggregated in one panel. - - internal cluster indices as one heatmap with clusterings (and [metadata_of_interest]) sorted by TOPSIS ranking from top to bottom and split cluster indices split by type (cost/benefit functions to be minimized/maximized). - + - all clustering results as 2D and interactive 2D & 3D plots for all available embedings/projections. + - external cluster indices as hierarchically clustered heatmaps, aggregated in one panel. + - internal cluster indices as one heatmap with clusterings (and [metadata_of_interest]) sorted by TOPSIS ranking from top to bottom and split cluster indices split by type (cost/benefit functions to be minimized/maximized). # 🛠️ Usage + Here are some tips for the usage of this workflow: + - Start with minimal parameter combinations and without UMAP diagnostics and connectivity plots (they are computational expensive and slow). - Heatmaps require **a lot** of memory, hence options to reduce computational cost are provided and the memory allocation is solved dynamically based on retries. If an out-of-memory exception occurs the flag `--retries X` can be used to trigger automatic resubmission X times upon failure with X times the memory. - Clustification performance scales with available cores, i.e., more cores faster internal parallelization of Random Forest training & testing. @@ -182,21 +187,26 @@ Here are some tips for the usage of this workflow: - Usage as a module and a selection of the results can be found on the [MrBiomics Wiki on "Module Usage in Projects".](https://github.com/epigen/MrBiomics/wiki/Module-Usage-in-Projects) # ⚙️ Configuration + Detailed specifications can be found here [./config/README.md](./config/README.md) # 📖 Examples + We provide a minimal example of the analysis of the [UCI ML hand-written digits datasets](https://archive.ics.uci.edu/ml/datasets/Optical+Recognition+of+Handwritten+Digits) imported from [sklearn](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_digits.html) in the [test folder](./test/): + - config - - configuration: config/config.yaml - - sample annotation: test/config/digits_unsupervised_analysis_annotation.csv + - configuration: config/config.yaml + - sample annotation: test/config/digits_unsupervised_analysis_annotation.csv - data - - dataset (1797 observations, 64 features): digits_data.csv - - metadata (consisting of the ground truth label "target"): digits_labels.csv + - dataset (1797 observations, 64 features): digits_data.csv + - metadata (consisting of the ground truth label "target"): digits_labels.csv - results will be generated in the configured subfolder `./test/results/` - performance: on an HPC it took less than 7 minutes to complete a full run (with up to 32GB of memory per task) ## 🧬 Bioinformatics examples + Explore detailed bioinformatics use cases showcasing module usage in our comprehensive end-to-end [MrBiomics Recipes](https://github.com/epigen/MrBiomics?tab=readme-ov-file#-recipes), including data, configuration, annotation and results: + - [ATAC-seq Analysis Recipe](https://github.com/epigen/MrBiomics/wiki/ATAC%E2%80%90seq-Analysis-Recipe) - [RNA-seq Analysis Recipe](https://github.com/epigen/MrBiomics/wiki/RNA%E2%80%90seq-Analysis-Recipe) - [Integrative Analysis Recipe](https://github.com/epigen/MrBiomics/wiki/Integrative-Analysis-Recipe) @@ -204,58 +214,64 @@ Explore detailed bioinformatics use cases showcasing module usage in our compreh - [scCRISPR-seq Analysis Recipe](https://github.com/epigen/MrBiomics/wiki/scCRISPR%E2%80%90seq-Analysis-Recipe) # 🧬 Single-cell RNA sequencing (scRNA-seq) data analysis + Unsupervised analyses, dimensionality reduction, and cluster analysis are cornerstones of scRNA-seq data analyses. A full run on a [published](https://www.nature.com/articles/s41588-020-0636-z) scRNA-seq [cancer dataset](https://www.weizmann.ac.il/sites/3CA/colorectal) with 21,657 cells and 18,245 genes took 2.5 hours to complete (without heatmaps, with 32GB memory and 8 cores for clustification). Below are configurations of the two most commonly used frameworks, [scanpy](https://scanpy.readthedocs.io/en/stable/index.html) (Python) and [Seurat](https://satijalab.org/seurat/) (R), and the original package's defaults as comparison and to facilitate reproducibility: UMAP for dimensionality reduction + - [umap-learn](https://umap-learn.readthedocs.io/en/latest/api.html) - - initialization: spectral - - metric: Euclidean - - neighbors: 15 - - min. distance: 0.1 + - initialization: spectral + - metric: Euclidean + - neighbors: 15 + - min. distance: 0.1 - [scanpy](https://scanpy.readthedocs.io/en/stable/generated/scanpy.pp.neighbors.html#scanpy.pp.neighbors) - - initialization: spectral - - metric: Euclidean - - neighbors: 15 - - min. distance: **0.5** + - initialization: spectral + - metric: Euclidean + - neighbors: 15 + - min. distance: **0.5** - [Seurat](https://satijalab.org/seurat/reference/runumap) - - initialization: **PCA** - - method: "uwot" (not umap-learn package) - - metric: **Cosine (or Correlation)** - - neighbors: **30** - - min. distance: **0.3** + - initialization: **PCA** + - method: "uwot" (not umap-learn package) + - metric: **Cosine (or Correlation)** + - neighbors: **30** + - min. distance: **0.3** Leiden algorithm for clustering + - [leidenalg](https://leidenalg.readthedocs.io/en/stable/reference.html) - - no defaults + - no defaults - [scanpy](https://scanpy.readthedocs.io/en/stable/generated/scanpy.tl.leiden.html) - - input: batch-balanced UMAP KNN graph - - partition type: RBConfigurationVertexPartition - - resolution: 1 + - input: batch-balanced UMAP KNN graph + - partition type: RBConfigurationVertexPartition + - resolution: 1 - [Seurat](https://github.com/satijalab/seurat/blob/763259d05991d40721dee99c9919ec6d4491d15e/R/clustering.R#L344) - - input: SNN graph - - partition type: RBConfigurationVertexPartition - - resolution: 0.8 - + - input: SNN graph + - partition type: RBConfigurationVertexPartition + - resolution: 0.8 # 🔗 Links + - [GitHub Repository](https://github.com/epigen/unsupervised_analysis/) - [GitHub Page](https://epigen.github.io/unsupervised_analysis/) - [Zenodo Repository](https://doi.org/10.5281/zenodo.8405360) - [Snakemake Workflow Catalog Entry](https://snakemake.github.io/snakemake-workflow-catalog?usage=epigen/unsupervised_analysis) # 📚 Resources + - Recommended compatible [MrBiomics Modules](https://github.com/epigen/MrBiomics/#-modules) for upstream analyses: - - [ATAC-seq Processing](https://github.com/epigen/atacseq_pipeline) to quantify chromatin accessibility. - - [scRNA-seq Data Processing & Visualization](https://github.com/epigen/scrnaseq_processing_seurat) for processing (multimodal) single-cell transcriptome data. - - [Split, Filter, Normalize and Integrate Sequencing Data](https://github.com/epigen/spilterlize_integrate/) after count quantification. - - [Differential Analysis with limma](https://github.com/epigen/dea_limma) to identify and visualize statistically significantly different features (e.g., genes or genomic regions) between sample groups. - - [Perturbation Analysis using Mixscape from Seurat](https://github.com/epigen/mixscape_seurat) to identify perturbed cells from pooled (multimodal) CRISPR screens with sc/snRNA-seq read-out (scCRISPR-seq). + - [ATAC-seq Processing](https://github.com/epigen/atacseq_pipeline) to quantify chromatin accessibility. + - [scRNA-seq Data Processing & Visualization](https://github.com/epigen/scrnaseq_processing_seurat) for processing (multimodal) single-cell transcriptome data. + - [Split, Filter, Normalize and Integrate Sequencing Data](https://github.com/epigen/spilterlize_integrate/) after count quantification. + - [Differential Analysis with limma](https://github.com/epigen/dea_limma) to identify and visualize statistically significantly different features (e.g., genes or genomic regions) between sample groups. + - [Perturbation Analysis using Mixscape from Seurat](https://github.com/epigen/mixscape_seurat) to identify perturbed cells from pooled (multimodal) CRISPR screens with sc/snRNA-seq read-out (scCRISPR-seq). - [Reichl, S. (2018). Mathematical methods in single cell RNA sequencing analysis with an emphasis on the validation of clustering results [Diploma Thesis, Technische Universität Wien]](https://doi.org/10.34726/hss.2018.49662) # 📑 Publications + The following publications successfully used this module for their analyses. + - [Traxler, Reichl et al. (2025) Cell Systems - Integrated time-series analysis and high-content CRISPR screening delineate the dynamics of macrophage immune regulation](https://doi.org/10.1016/j.cels.2025.101346) - ... diff --git a/_config.yml b/_config.yml index c419263..277f1f2 100644 --- a/_config.yml +++ b/_config.yml @@ -1 +1 @@ -theme: jekyll-theme-cayman \ No newline at end of file +theme: jekyll-theme-cayman diff --git a/config/README.md b/config/README.md index 0022eac..377fe7e 100644 --- a/config/README.md +++ b/config/README.md @@ -4,9 +4,9 @@ You need one configuration file to configure the analyses and one annotation fil - project configuration (`config/config.yaml`): Different for every project and configures the analyses to be performed. - sample annotation (annotation): CSV file consisting of four mandatory columns. - - name: A unique name for the dataset (tip: keep it short but descriptive). - - data: Path to the tabular data as a comma-separated table (CSV). - - metadata: Path to the metadata as a comma-separated table (CSV) with the first column being the index/identifier of each observation/sample and every other column metadata for the respective observation (either numeric or categorical, not mixed). **No NaN or empty values allowed, and no special characters (all except a-z, 0-9, `_`) in the index.** - - samples_by_features: Boolean indicator if the data matrix is observations/samples (rows) x features (columns): 0==no, 1==yes. + - name: A unique name for the dataset (tip: keep it short but descriptive). + - data: Path to the tabular data as a comma-separated table (CSV). + - metadata: Path to the metadata as a comma-separated table (CSV) with the first column being the index/identifier of each observation/sample and every other column metadata for the respective observation (either numeric or categorical, not mixed). **No NaN or empty values allowed, and no special characters (all except a-z, 0-9, `_`) in the index.** + - samples_by_features: Boolean indicator if the data matrix is observations/samples (rows) x features (columns): 0==no, 1==yes. Set workflow-specific `resources` or command line arguments (CLI) in the workflow profile `workflow/profiles/default.config.yaml`, which supersedes global Snakemake profiles. diff --git a/config/config.yaml b/config/config.yaml index 6c62b1f..4877bd3 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -2,7 +2,7 @@ ##### RESOURCES ##### # memory in MB -mem: '32000' +mem: 32000 threads: 2 ##### GENERAL ##### @@ -14,8 +14,8 @@ project_name: digits # https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html # especially relevant for large data pca: - n_components: 0.9 # variance as float (0-1], number of components as int e.g., 50, or 'mle' - svd_solver: 'auto' # options: ‘auto’, ‘full’, ‘covariance_eigh’, ‘arpack’, ‘randomized’ + n_components: 0.9 # variance as float (0-1], number of components as int e.g., 50, or 'mle' + svd_solver: "auto" # options: ‘auto’, ‘full’, ‘covariance_eigh’, ‘arpack’, ‘randomized’ ##### UMAP & densMAP ##### # https://umap-learn.readthedocs.io/en/latest/parameters.html @@ -25,13 +25,13 @@ pca: # densmap flag: perform densMAP (0==no, 1==yes) on top of UMAP with all parameter combinations (https://umap-learn.readthedocs.io/en/latest/densmap_demo.html) # connectivity visualization flag (0==no, 1==yes), computational expensive (slow), recommendation 0 for exploration and 1 for validation umap: - metrics: ['euclidean'] - n_neighbors: [15] - min_dist: [0.1] - n_components: [2] - densmap: 1 - connectivity: 1 - diagnostics: 1 + metrics: ["euclidean"] + n_neighbors: [15] + min_dist: [0.1] + n_components: [2] + densmap: 1 + connectivity: 1 + diagnostics: 1 ##### HEATMAP ##### # information on the ComplexHeatmap parameters: https://jokergoo.github.io/ComplexHeatmap-reference/book/index.html @@ -39,38 +39,39 @@ umap: # clustering methods: methods for hierarchical clustering that are supported by fastcluster's R implementation (https://danifold.net/fastcluster.html) # it is the most resource (memory) intensive method, leave empty [] if not required heatmap: - metrics: ['correlation','cosine'] - hclust_methods: ['complete'] - n_observations: 1 # random sampled proportion float (0-1] or absolute number as integer - n_features: 0.5 # highly variable features proportion float (0-1] or absolute number as integer - + metrics: ["correlation", "cosine"] + hclust_methods: ["complete"] + n_observations: 1 # random sampled proportion float (0-1] or absolute number as integer + n_features: 0.5 # highly variable features proportion float (0-1] or absolute number as integer + ##### LEIDEN ##### # Leiden clustering applied on UMAP KNN graphs specified by the respective parameters (metric, n_neighbors). # Leiden algorithm specific parameters (partition_types, resolutions, n_iterations) -> https://leidenalg.readthedocs.io/en/stable/index.html # To skip Leiden clustering, leave the "metrics" parameter empty i.e., [] leiden: - metrics: [] # has to be a subset of above's UMAP parameters - n_neighbors: [15] # has to be a subset of above's UMAP parameters - partition_types: ["RBConfigurationVertexPartition", "ModularityVertexPartition"] - resolutions: [0.5,1,1.5,2,4] # only used for relevant partition_types - n_iterations: 2 # default: 2; -1 until convergence + metrics: [] # has to be a subset of above's UMAP parameters + n_neighbors: [15] # has to be a subset of above's UMAP parameters + partition_types: + ["RBConfigurationVertexPartition", "ModularityVertexPartition"] + resolutions: [0.5, 1, 1.5, 2, 4] # only used for relevant partition_types + n_iterations: 2 # default: 2; -1 until convergence ##### clustification ##### # ML-based clustering approach that iteratively merges clusters based on misclassification. # Doesn't support externally provided clusterings. clustification: - method: "" # starting clustering result method, options: "Leiden" + method: "" # starting clustering result method, options: "Leiden" ##### clustree ##### # Cluster analysis and visualization using clustree: https://lazappi.github.io/clustree/index.html # Clustree specific parameters (count_filter, prop_filter, layout): https://lazappi.github.io/clustree/reference/clustree.html # Will be skipped if no clustering results are available clustree: - count_filter: 0 # default: 0 - prop_filter: 0.1 # default: 0.1 - layout: "tree" # options: "tree" or "sugiyama" - categorical_label_option: "majority" # aggregation function for categorical metadata, options: "pure" or "majority" - numerical_aggregation_option: "mean" # aggregation function for numerical metadata, options: mean, median, max, min + count_filter: 0 # default: 0 + prop_filter: 0.1 # default: 0.1 + layout: "tree" # options: "tree" or "sugiyama" + categorical_label_option: "majority" # aggregation function for categorical metadata, options: "pure" or "majority" + numerical_aggregation_option: "mean" # aggregation function for numerical metadata, options: mean, median, max, min ##### cluster validation #### # Cluster validation using internal cluster indices is computationally very expensive. @@ -85,21 +86,21 @@ sample_proportion: 1 # float [0-1], >500 samples should be included. # - internal cluster validation (categorical metadata interpreted as clusterings) # Can be empty [], then the first column of the metadata dataframe will be taken. metadata_of_interest: ["target"] - + ##### VISUALIZATION ##### -# flag if coordinates should be fixed in 2D plots by +coord_fixed() (0==no, 1==yes) +# flag if coordinates should be fixed in 2D plots by +coord_fixed() (0==no, 1==yes) # https://ggplot2.tidyverse.org/reference/coord_fixed.html coord_fixed: 0 # 2D/3D visualization with ggplot2 and plotly # for more than 200 datapoints it is recommended to reduce size to 0.1 scatterplot2d: - size: 1 - alpha: 1 + size: 1 + alpha: 1 # specify features of interest. these features from the data, will be highlighted in the 2D/3D plots # motivated by bioinformatics highlighting expression levels of marker genes (eg: ['PTPRC','STAT1','IRF8']) # use keyword ['ALL'] to plot all features. WARNING: Only useful for relatively low dimensional data, a plot is generated for each feature and method. # if not used leave empty [] -features_to_plot: ['ALL'] #['pixel_0_0','pixel_0_1','pixel_0_2','pixel_0_3'] +features_to_plot: ["ALL"] #['pixel_0_0','pixel_0_1','pixel_0_2','pixel_0_3'] diff --git a/workflow/Snakefile b/workflow/Snakefile index 09e3d36..25f93d7 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -1,4 +1,3 @@ - ##### global workflow dependencies ##### # conda: "envs/global.yaml" @@ -15,16 +14,19 @@ min_version("8.20.1") ##### module name ##### module_name = "unsupervised_analysis" + ##### setup report ##### report: os.path.join("report", "workflow.rst") + ##### load config and sample annotation sheets ##### configfile: os.path.join("config", "config.yaml") + # load annotations -annot = pd.read_csv(config['annotation'], index_col='name') +annot = pd.read_csv(config["annotation"], index_col="name") -result_path = os.path.join(config["result_path"],module_name) +result_path = os.path.join(config["result_path"], module_name) # Clustering: Leiden, clustification cluster_methods = [] @@ -32,151 +34,297 @@ clustree_content = [] clustree_metadata_content = [] ## Leiden -if len(config["leiden"]["metrics"])>0 and all([leiden_metric in config["umap"]["metrics"] for leiden_metric in config["leiden"]["metrics"]]): +if len(config["leiden"]["metrics"]) > 0 and all( + [ + leiden_metric in config["umap"]["metrics"] + for leiden_metric in config["leiden"]["metrics"] + ] +): cluster_methods.append("Leiden") ## clustification -if len(config["clustification"]["method"])!="" and config["clustification"]["method"] in cluster_methods: +if ( + len(config["clustification"]["method"]) != "" + and config["clustification"]["method"] in cluster_methods +): cluster_methods.append("clustification") # Cluster validation indices_external = [] indices_internal = [] - + # Dimensionality Reduction: apply PCA, UMAP and densMAP sepcific visualization flags -umap_content = ['metadata'] +umap_content = ["metadata"] umap_diagnostics = [] # densmap_content = ['metadata'] -pca_content = ['metadata'] -pca_diagnostics = ['variance','pairs','loadings','loadings_lollipop'] +pca_content = ["metadata"] +pca_diagnostics = ["variance", "pairs", "loadings", "loadings_lollipop"] -if config["umap"]["connectivity"]==1: - umap_diagnostics.append('connectivity') -# densmap_content.append('connectivity') +if config["umap"]["connectivity"] == 1: + umap_diagnostics.append("connectivity") + # densmap_content.append('connectivity') -if config["umap"]["diagnostics"]==1: - umap_diagnostics.append('diagnostics') -# densmap_content.append('diagnostics') +if config["umap"]["diagnostics"] == 1: + umap_diagnostics.append("diagnostics") + # densmap_content.append('diagnostics') if len(config["features_to_plot"]) > 0: - pca_content.append('features') - umap_content.append('features') -# densmap_content.append('features') + pca_content.append("features") + umap_content.append("features") + # densmap_content.append('features') if len(cluster_methods) > 0: - pca_content.append('clustering') - umap_content.append('clustering') + pca_content.append("clustering") + umap_content.append("clustering") clustree_content.append("default") clustree_metadata_content.append("metadata") - indices_external = ["AMI","ARI","FMI", "Homogeneity", "Completeness", "V"] - indices_internal = ["Silhouette", "Calinski_Harabasz", "Dunn", "C_index", "Davies_Bouldin", "BIC"] # removed "AIC" - + indices_external = ["AMI", "ARI", "FMI", "Homogeneity", "Completeness", "V"] + indices_internal = [ + "Silhouette", + "Calinski_Harabasz", + "Dunn", + "C_index", + "Davies_Bouldin", + "BIC", + ] # removed "AIC" + if len(config["features_to_plot"]) > 0: - clustree_metadata_content.append('features') - - if len(config["metadata_of_interest"])>0: - clustree_content.append('custom') + clustree_metadata_content.append("features") + + if len(config["metadata_of_interest"]) > 0: + clustree_content.append("custom") + ruleorder: clustification > aggregate_clustering_results ruleorder: clustree_analysis > clustree_analysis_metadata + rule all: input: # PCA - PCA_objects = expand(os.path.join(result_path,'{sample}','PCA','PCA_{svd_solver}_{n_components}_object.pickle'), - sample=list(annot.index), - svd_solver=config["pca"]["svd_solver"], - n_components=config["pca"]["n_components"], - ), - PCA_plots = expand(os.path.join(result_path,'{sample}','PCA','plots','PCA_{svd_solver}_{n_components}_2','{content}'), - sample=list(annot.index), - svd_solver=config["pca"]["svd_solver"], - n_components=config["pca"]["n_components"], - content=pca_content, - ), - PCA_diagnostics = expand(os.path.join(result_path,'{sample}','PCA','plots','PCA_{svd_solver}_{n_components}','diagnostics','{content}.png'), - sample=list(annot.index), - svd_solver=config["pca"]["svd_solver"], - n_components=config["pca"]["n_components"], - content=pca_diagnostics, - ), - PCA_interactive = expand(os.path.join(result_path,'{sample}','PCA','plots','PCA_{svd_solver}_{n_components}_{dim}','interactive.html'), - sample=list(annot.index), - svd_solver=config["pca"]["svd_solver"], - n_components=config["pca"]["n_components"], - dim=[2,3] - ), + PCA_objects=expand( + os.path.join( + result_path, + "{sample}", + "PCA", + "PCA_{svd_solver}_{n_components}_object.pickle", + ), + sample=list(annot.index), + svd_solver=config["pca"]["svd_solver"], + n_components=config["pca"]["n_components"], + ), + PCA_plots=expand( + os.path.join( + result_path, + "{sample}", + "PCA", + "plots", + "PCA_{svd_solver}_{n_components}_2", + "{content}", + ), + sample=list(annot.index), + svd_solver=config["pca"]["svd_solver"], + n_components=config["pca"]["n_components"], + content=pca_content, + ), + PCA_diagnostics=expand( + os.path.join( + result_path, + "{sample}", + "PCA", + "plots", + "PCA_{svd_solver}_{n_components}", + "diagnostics", + "{content}.png", + ), + sample=list(annot.index), + svd_solver=config["pca"]["svd_solver"], + n_components=config["pca"]["n_components"], + content=pca_diagnostics, + ), + PCA_interactive=expand( + os.path.join( + result_path, + "{sample}", + "PCA", + "plots", + "PCA_{svd_solver}_{n_components}_{dim}", + "interactive.html", + ), + sample=list(annot.index), + svd_solver=config["pca"]["svd_solver"], + n_components=config["pca"]["n_components"], + dim=[2, 3], + ), # UMAP & densMAP - MAP_graphs = expand(os.path.join(result_path,'{sample}','UMAP','UMAP_{metric}_{n_neighbors}_graph.pickle'), - sample=list(annot.index), - metric=config["umap"]["metrics"], - n_neighbors=config["umap"]["n_neighbors"] #max(config["umap"]["n_neighbors"]) - ), - MAP_objects = expand(os.path.join(result_path,'{sample}','{method}','{method}_{metric}_{n_neighbors}_{min_dist}_{n_components}_object.pickle'), - sample=list(annot.index), - method=["UMAP","densMAP"] if config["umap"]["densmap"]==1 else ["UMAP"], - metric=config["umap"]["metrics"], - n_neighbors=config["umap"]["n_neighbors"], - min_dist=config["umap"]["min_dist"], - n_components=config["umap"]["n_components"] - ), - MAP_plots = expand(os.path.join(result_path,'{sample}','{method}','plots','{method}_{metric}_{n_neighbors}_{min_dist}_{n_components}','{content}'), - sample=list(annot.index), - method=["UMAP","densMAP"] if config["umap"]["densmap"]==1 else ["UMAP"], - metric=config["umap"]["metrics"], - n_neighbors=config["umap"]["n_neighbors"], - min_dist=config["umap"]["min_dist"], - n_components=2, - content=umap_content - ) if 2 in config["umap"]["n_components"] else [], - MAP_diagnostics = expand(os.path.join(result_path,'{sample}','{method}','plots','{method}_{metric}_{n_neighbors}_{min_dist}_{n_components}','diagnostics','{content}.png'), - sample=list(annot.index), - method=["UMAP","densMAP"] if config["umap"]["densmap"]==1 else ["UMAP"], - metric=config["umap"]["metrics"], - n_neighbors=config["umap"]["n_neighbors"], - min_dist=config["umap"]["min_dist"], - n_components=2, - content=umap_diagnostics - ) if 2 in config["umap"]["n_components"] else [], - MAP_interactive = expand(os.path.join(result_path,'{sample}','{method}','plots','{method}_{metric}_{n_neighbors}_{min_dist}_{n_components}','interactive.html'), - sample=list(annot.index), - method=["UMAP","densMAP"] if config["umap"]["densmap"]==1 else ["UMAP"], - metric=config["umap"]["metrics"], - n_neighbors=config["umap"]["n_neighbors"], - min_dist=config["umap"]["min_dist"], - n_components=[dims for dims in config["umap"]["n_components"] if dims in [2,3]] - ) if 2 in config["umap"]["n_components"] or 3 in config["umap"]["n_components"] else [], + MAP_graphs=expand( + os.path.join( + result_path, + "{sample}", + "UMAP", + "UMAP_{metric}_{n_neighbors}_graph.pickle", + ), + sample=list(annot.index), + metric=config["umap"]["metrics"], + n_neighbors=config["umap"][ + "n_neighbors" + ], # max(config["umap"]["n_neighbors"]) + ), + MAP_objects=expand( + os.path.join( + result_path, + "{sample}", + "{method}", + "{method}_{metric}_{n_neighbors}_{min_dist}_{n_components}_object.pickle", + ), + sample=list(annot.index), + method=["UMAP", "densMAP"] if config["umap"]["densmap"] == 1 else ["UMAP"], + metric=config["umap"]["metrics"], + n_neighbors=config["umap"]["n_neighbors"], + min_dist=config["umap"]["min_dist"], + n_components=config["umap"]["n_components"], + ), + MAP_plots=( + expand( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{metric}_{n_neighbors}_{min_dist}_{n_components}", + "{content}", + ), + sample=list(annot.index), + method=( + ["UMAP", "densMAP"] if config["umap"]["densmap"] == 1 else ["UMAP"] + ), + metric=config["umap"]["metrics"], + n_neighbors=config["umap"]["n_neighbors"], + min_dist=config["umap"]["min_dist"], + n_components=2, + content=umap_content, + ) + if 2 in config["umap"]["n_components"] + else [] + ), + MAP_diagnostics=( + expand( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{metric}_{n_neighbors}_{min_dist}_{n_components}", + "diagnostics", + "{content}.png", + ), + sample=list(annot.index), + method=( + ["UMAP", "densMAP"] if config["umap"]["densmap"] == 1 else ["UMAP"] + ), + metric=config["umap"]["metrics"], + n_neighbors=config["umap"]["n_neighbors"], + min_dist=config["umap"]["min_dist"], + n_components=2, + content=umap_diagnostics, + ) + if 2 in config["umap"]["n_components"] + else [] + ), + MAP_interactive=( + expand( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{metric}_{n_neighbors}_{min_dist}_{n_components}", + "interactive.html", + ), + sample=list(annot.index), + method=( + ["UMAP", "densMAP"] if config["umap"]["densmap"] == 1 else ["UMAP"] + ), + metric=config["umap"]["metrics"], + n_neighbors=config["umap"]["n_neighbors"], + min_dist=config["umap"]["min_dist"], + n_components=[ + dims for dims in config["umap"]["n_components"] if dims in [2, 3] + ], + ) + if 2 in config["umap"]["n_components"] + or 3 in config["umap"]["n_components"] + else [] + ), # Heatmap - heatmap_plots = expand(os.path.join(result_path,'{sample}','Heatmap','plots','Heatmap_{metric}_{method}.png'), - sample=list(annot.index), - method=config["heatmap"]["hclust_methods"], - metric=config["heatmap"]["metrics"], - ), + heatmap_plots=expand( + os.path.join( + result_path, + "{sample}", + "Heatmap", + "plots", + "Heatmap_{metric}_{method}.png", + ), + sample=list(annot.index), + method=config["heatmap"]["hclust_methods"], + metric=config["heatmap"]["metrics"], + ), # clustree - clustree_plots = expand(os.path.join(result_path,'{sample}','clustree','{content}.png'), - sample=list(annot.index), - content=clustree_content - ), - clustree_metadata_plots = expand(os.path.join(result_path,'{sample}','clustree','{content}'), - sample=list(annot.index), - content=clustree_metadata_content - ), + clustree_plots=expand( + os.path.join(result_path, "{sample}", "clustree", "{content}.png"), + sample=list(annot.index), + content=clustree_content, + ), + clustree_metadata_plots=expand( + os.path.join(result_path, "{sample}", "clustree", "{content}"), + sample=list(annot.index), + content=clustree_metadata_content, + ), # cluster validation with internal and external indices - cluster_validation = expand(os.path.join(result_path,'{sample}','cluster_validation','plots','{index_type}'), - sample=list(annot.index), - index_type = ["external", "internal"] if config["sample_proportion"]>0 else ["external"], - ) if len(cluster_methods)>0 else [], - envs = expand(os.path.join(result_path,'envs','{env}.yaml'),env=['clusterCrit','clustree','ComplexHeatmap','ggplot','umap_leiden','plotly','pymcdm']), - configs = os.path.join(result_path,'configs','{}_config.yaml'.format(config["project_name"])), - annotations = os.path.join(result_path,'configs','{}_annot.csv'.format(config["project_name"])), + cluster_validation=( + expand( + os.path.join( + result_path, + "{sample}", + "cluster_validation", + "plots", + "{index_type}", + ), + sample=list(annot.index), + index_type=( + ["external", "internal"] + if config["sample_proportion"] > 0 + else ["external"] + ), + ) + if len(cluster_methods) > 0 + else [] + ), + envs=expand( + os.path.join(result_path, "envs", "{env}.yaml"), + env=[ + "clusterCrit", + "clustree", + "ComplexHeatmap", + "ggplot", + "umap_leiden", + "plotly", + "pymcdm", + ], + ), + configs=os.path.join( + result_path, "configs", "{}_config.yaml".format(config["project_name"]) + ), + annotations=os.path.join( + result_path, "configs", "{}_annot.csv".format(config["project_name"]) + ), resources: mem_mb=1000, threads: config.get("threads", 1) log: - os.path.join("logs","rules","all.log"), + os.path.join("logs", "rules", "all.log"), - ##### load rules ##### include: os.path.join("rules", "common.smk") include: os.path.join("rules", "dimred.smk") diff --git a/workflow/envs/ComplexHeatmap.yaml b/workflow/envs/ComplexHeatmap.yaml index c60bd93..f792079 100644 --- a/workflow/envs/ComplexHeatmap.yaml +++ b/workflow/envs/ComplexHeatmap.yaml @@ -8,4 +8,4 @@ dependencies: - r-rcolorbrewer=1.1_3 - r-fastcluster=1.2.3 - r-magick=2.7.3 - - r-data.table=1.14.8 \ No newline at end of file + - r-data.table=1.14.8 diff --git a/workflow/envs/clusterCrit.yaml b/workflow/envs/clusterCrit.yaml index 6461e21..b03a756 100644 --- a/workflow/envs/clusterCrit.yaml +++ b/workflow/envs/clusterCrit.yaml @@ -3,4 +3,4 @@ channels: - nodefaults dependencies: - r-clustercrit=1.2.8 - - r-data.table=1.15.4 \ No newline at end of file + - r-data.table=1.15.4 diff --git a/workflow/envs/clustree.yaml b/workflow/envs/clustree.yaml index ccf1eac..e0b1f5b 100644 --- a/workflow/envs/clustree.yaml +++ b/workflow/envs/clustree.yaml @@ -5,4 +5,4 @@ channels: dependencies: - r-clustree=0.5.0 - r-patchwork=1.1.2 - - r-data.table=1.15.4 \ No newline at end of file + - r-data.table=1.15.4 diff --git a/workflow/envs/ggplot.yaml b/workflow/envs/ggplot.yaml index 0732d20..b457b72 100644 --- a/workflow/envs/ggplot.yaml +++ b/workflow/envs/ggplot.yaml @@ -9,4 +9,4 @@ dependencies: - r-ggrepel=0.9.1 - r-reshape2=1.4.4 - r-stringi=1.8.4 - - r-data.table=1.15.4 \ No newline at end of file + - r-data.table=1.15.4 diff --git a/workflow/envs/global.yaml b/workflow/envs/global.yaml index 482aa8e..55e5263 100644 --- a/workflow/envs/global.yaml +++ b/workflow/envs/global.yaml @@ -4,4 +4,4 @@ channels: - nodefaults dependencies: - numpy=2.0.1 - - pandas=2.2.2 \ No newline at end of file + - pandas=2.2.2 diff --git a/workflow/envs/plotly.yaml b/workflow/envs/plotly.yaml index 0d45fad..0dc3368 100644 --- a/workflow/envs/plotly.yaml +++ b/workflow/envs/plotly.yaml @@ -7,4 +7,4 @@ dependencies: - plotly=5.10 - plotly_express=0.4.1 - pandas=1.5.0 - - seaborn-base=0.12.0 \ No newline at end of file + - seaborn-base=0.12.0 diff --git a/workflow/envs/pymcdm.yaml b/workflow/envs/pymcdm.yaml index d93d36c..d9e24cb 100644 --- a/workflow/envs/pymcdm.yaml +++ b/workflow/envs/pymcdm.yaml @@ -4,4 +4,4 @@ channels: dependencies: - pymcdm=1.1.0 - pandas=1.5.0 - - numpy=1.23.3 \ No newline at end of file + - numpy=1.23.3 diff --git a/workflow/envs/umap_leiden.yaml b/workflow/envs/umap_leiden.yaml index 0c77608..f8d0b34 100644 --- a/workflow/envs/umap_leiden.yaml +++ b/workflow/envs/umap_leiden.yaml @@ -16,5 +16,5 @@ dependencies: - dask=2024.2.1 - scikit-image=0.24.0 - pip: - - umap-learn==0.5.6 - - umap-learn[plot] + - umap-learn==0.5.6 + - umap-learn[plot] diff --git a/workflow/profiles/default/config.yaml b/workflow/profiles/default/config.yaml index 29bebdd..856b930 100644 --- a/workflow/profiles/default/config.yaml +++ b/workflow/profiles/default/config.yaml @@ -1,3 +1,3 @@ default-resources: - slurm_partition: shortq - slurm_extra: "'--qos=shortq'" \ No newline at end of file + slurm_partition: shortq + slurm_extra: "'--qos=shortq'" diff --git a/workflow/rules/cluster_validation.smk b/workflow/rules/cluster_validation.smk index 06d421d..6c9c7ae 100644 --- a/workflow/rules/cluster_validation.smk +++ b/workflow/rules/cluster_validation.smk @@ -1,124 +1,175 @@ - # clustree analysis rule clustree_analysis: input: unpack(get_clustree_paths), output: - plot = report(os.path.join(config["result_path"],'unsupervised_analysis','{sample}','clustree','{content}.png'), - caption="../report/clustree.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "clustree", - "parameters": "-", - "dimensions": "-", - "type": "clustree", - "content": "{content}", - } - ), + plot=report( + os.path.join( + config["result_path"], + "unsupervised_analysis", + "{sample}", + "clustree", + "{content}.png", + ), + caption="../report/clustree.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "clustree", + "parameters": "-", + "dimensions": "-", + "type": "clustree", + "content": "{content}", + }, + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/clustree.yaml" log: - os.path.join("logs","rules","clustree_{sample}_{content}.log"), + os.path.join("logs", "rules", "clustree_{sample}_{content}.log"), params: - content = lambda w: "{}".format(w.content), - count_filter = config["clustree"]["count_filter"], - prop_filter = config["clustree"]["prop_filter"], - layout = config["clustree"]["layout"], - categorical_label_option = config["clustree"]["categorical_label_option"], - numerical_aggregation_option = config["clustree"]["numerical_aggregation_option"], - custom_metadata = config["metadata_of_interest"], + content=lambda w: "{}".format(w.content), + count_filter=config["clustree"]["count_filter"], + prop_filter=config["clustree"]["prop_filter"], + layout=config["clustree"]["layout"], + categorical_label_option=config["clustree"]["categorical_label_option"], + numerical_aggregation_option=config["clustree"]["numerical_aggregation_option"], + custom_metadata=config["metadata_of_interest"], script: "../scripts/clustree.R" - + + # clustree analysis for highlighting individual metadata and features rule clustree_analysis_metadata: input: unpack(get_clustree_paths), output: - plot = report(directory(os.path.join(config["result_path"],'unsupervised_analysis','{sample}','clustree','{content}')), - patterns=["{metadata}.png"], - caption="../report/clustree.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "clustree", - "parameters": "-", - "dimensions": "-", - "type": "{content}", - "content": "{metadata}", - } - ), + plot=report( + directory( + os.path.join( + config["result_path"], + "unsupervised_analysis", + "{sample}", + "clustree", + "{content}", + ) + ), + patterns=["{metadata}.png"], + caption="../report/clustree.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "clustree", + "parameters": "-", + "dimensions": "-", + "type": "{content}", + "content": "{metadata}", + }, + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/clustree.yaml" log: - os.path.join("logs","rules","clustree_{sample}_{content}.log"), + os.path.join("logs", "rules", "clustree_{sample}_{content}.log"), params: - content = lambda w: "{}".format(w.content), - count_filter = config["clustree"]["count_filter"], - prop_filter = config["clustree"]["prop_filter"], - layout = config["clustree"]["layout"], - categorical_label_option = config["clustree"]["categorical_label_option"], - numerical_aggregation_option = config["clustree"]["numerical_aggregation_option"], - custom_metadata = config["metadata_of_interest"], + content=lambda w: "{}".format(w.content), + count_filter=config["clustree"]["count_filter"], + prop_filter=config["clustree"]["prop_filter"], + layout=config["clustree"]["layout"], + categorical_label_option=config["clustree"]["categorical_label_option"], + numerical_aggregation_option=config["clustree"]["numerical_aggregation_option"], + custom_metadata=config["metadata_of_interest"], script: "../scripts/clustree.R" + # determine external cluster indices rule validation_external: input: unpack(get_external_validation_paths), output: - expand(os.path.join(config["result_path"],"unsupervised_analysis","{{sample}}", "cluster_validation", "external_index_{index}.csv"), index=indices_external), + expand( + os.path.join( + config["result_path"], + "unsupervised_analysis", + "{{sample}}", + "cluster_validation", + "external_index_{index}.csv", + ), + index=indices_external, + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","validation_external_{sample}.log"), + os.path.join("logs", "rules", "validation_external_{sample}.log"), script: "../scripts/validation_external.py" - + + # determine internal cluster indices rule validation_internal: input: unpack(get_internal_validation_paths), output: - internal_indices = os.path.join(config["result_path"],"unsupervised_analysis","{sample}", "cluster_validation", "internal_index_{internal_index}.csv"), + internal_indices=os.path.join( + config["result_path"], + "unsupervised_analysis", + "{sample}", + "cluster_validation", + "internal_index_{internal_index}.csv", + ), resources: - mem_mb=2*int(config.get("mem", "16000")), + mem_mb=2 * int(config.get("mem", "16000")), threads: config.get("threads", 1) conda: "../envs/clusterCrit.yaml" log: - os.path.join("logs","rules","validation_internal_{internal_index}_{sample}.log"), + os.path.join( + "logs", "rules", "validation_internal_{internal_index}_{sample}.log" + ), params: - internal_index = lambda w: "{}".format(w.internal_index), - sample_proportion = config["sample_proportion"], - metadata_of_interest = config["metadata_of_interest"], + internal_index=lambda w: "{}".format(w.internal_index), + sample_proportion=config["sample_proportion"], + metadata_of_interest=config["metadata_of_interest"], script: "../scripts/validation_internal.R" - + + # rank internal cluster indices using MCDM method TOPSIS rule aggregate_rank_internal: input: - expand(os.path.join(config["result_path"],"unsupervised_analysis","{{sample}}", "cluster_validation", "internal_index_{index}.csv"), index=indices_internal), + expand( + os.path.join( + config["result_path"], + "unsupervised_analysis", + "{{sample}}", + "cluster_validation", + "internal_index_{index}.csv", + ), + index=indices_internal, + ), output: - internal_indices_ranked = os.path.join(config["result_path"],"unsupervised_analysis","{sample}", "cluster_validation", "internal_indices_ranked.csv"), + internal_indices_ranked=os.path.join( + config["result_path"], + "unsupervised_analysis", + "{sample}", + "cluster_validation", + "internal_indices_ranked.csv", + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/pymcdm.yaml" log: - os.path.join("logs","rules","rank_internal_{sample}.log"), + os.path.join("logs", "rules", "rank_internal_{sample}.log"), script: "../scripts/mcdm_topsis.py" @@ -128,25 +179,35 @@ rule plot_indices: input: unpack(get_validation_paths), output: - plot = report(directory(os.path.join(config["result_path"],'unsupervised_analysis','{sample}','cluster_validation','plots','{type}')), - patterns=["{index}.png"], - caption="../report/cluster_validation.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "cluster validation", - "parameters": "-", - "dimensions": "-", - "type": "{type}", - "content": "{index}", - } - ), + plot=report( + directory( + os.path.join( + config["result_path"], + "unsupervised_analysis", + "{sample}", + "cluster_validation", + "plots", + "{type}", + ) + ), + patterns=["{index}.png"], + caption="../report/cluster_validation.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "cluster validation", + "parameters": "-", + "dimensions": "-", + "type": "{type}", + "content": "{index}", + }, + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/ggplot.yaml" log: - os.path.join("logs","rules","plot_{type}_indices_{sample}.log"), + os.path.join("logs", "rules", "plot_{type}_indices_{sample}.log"), script: - "../scripts/plot_indices.R" \ No newline at end of file + "../scripts/plot_indices.R" diff --git a/workflow/rules/clustering.smk b/workflow/rules/clustering.smk index a51355e..65250f0 100644 --- a/workflow/rules/clustering.smk +++ b/workflow/rules/clustering.smk @@ -3,77 +3,111 @@ rule leiden_cluster: input: get_umap_sample_paths, output: - clustering = os.path.join(config["result_path"],'unsupervised_analysis','{sample}','Leiden','Leiden_{metric}_{n_neighbors}_{partition_type}_{resolution}_clustering.csv'), + clustering=os.path.join( + config["result_path"], + "unsupervised_analysis", + "{sample}", + "Leiden", + "Leiden_{metric}_{n_neighbors}_{partition_type}_{resolution}_clustering.csv", + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","leiden_{sample}_{metric}_{n_neighbors}_{partition_type}_{resolution}_clustering.log"), + os.path.join( + "logs", + "rules", + "leiden_{sample}_{metric}_{n_neighbors}_{partition_type}_{resolution}_clustering.log", + ), params: - samples_by_features = get_data_orientation, - metric = lambda w: "{}".format(w.metric), - n_neighbors = lambda w: "{}".format(w.n_neighbors), - partition_type = lambda w: "{}".format(w.partition_type), - resolution = lambda w: "{}".format(w.resolution), - n_iterations = config["leiden"]["n_iterations"] + samples_by_features=get_data_orientation, + metric=lambda w: "{}".format(w.metric), + n_neighbors=lambda w: "{}".format(w.n_neighbors), + partition_type=lambda w: "{}".format(w.partition_type), + resolution=lambda w: "{}".format(w.resolution), + n_iterations=config["leiden"]["n_iterations"], script: "../scripts/leiden_cluster.py" - + + # perform clustification based on initial clustering rule clustification: input: get_clustification_paths, output: - clustering = os.path.join(config["result_path"],'unsupervised_analysis','{sample}','clustification','clustification_clusterings.csv'), + clustering=os.path.join( + config["result_path"], + "unsupervised_analysis", + "{sample}", + "clustification", + "clustification_clusterings.csv", + ), resources: mem_mb=config.get("mem", "16000"), - threads: 8 #config.get("threads", 1) + threads: 8 #config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","clustification_{sample}_clusterings.log"), + os.path.join("logs", "rules", "clustification_{sample}_clusterings.log"), params: - samples_by_features = get_data_orientation, + samples_by_features=get_data_orientation, script: "../scripts/clustification.py" - + + # aggregate clustering results per method rule aggregate_clustering_results: input: get_clustering_paths, output: - aggregated_clusterings = os.path.join(config["result_path"],'unsupervised_analysis','{sample}','{method}','{method}_clusterings.csv'), + aggregated_clusterings=os.path.join( + config["result_path"], + "unsupervised_analysis", + "{sample}", + "{method}", + "{method}_clusterings.csv", + ), resources: mem_mb=config.get("mem", "16000"), log: - os.path.join("logs","rules","aggregate_clustering_results_{sample}_{method}.log"), + os.path.join( + "logs", "rules", "aggregate_clustering_results_{sample}_{method}.log" + ), run: # list to hold the individual clusterings agg_clust = [] # read each clustering result and add to list for filename in input[1:]: - clust_tmp = pd.read_csv(filename, header=0, index_col=0)#.squeeze("columns") + clust_tmp = pd.read_csv( + filename, header=0, index_col=0 + ) # .squeeze("columns") agg_clust.append(clust_tmp) - # convert list to dataframe + # convert list to dataframe agg_clust_df = pd.concat(agg_clust, axis=1) - + # Write the DataFrame to a CSV file agg_clust_df.to_csv(output.aggregated_clusterings, index=True) - + + # aggregate clustering results across methods rule aggregate_all_clustering_results: input: get_aggregated_clustering_paths, output: - metadata_clusterings = os.path.join(config["result_path"],'unsupervised_analysis','{sample}','metadata_clusterings.csv'), + metadata_clusterings=os.path.join( + config["result_path"], + "unsupervised_analysis", + "{sample}", + "metadata_clusterings.csv", + ), resources: mem_mb=config.get("mem", "16000"), log: - os.path.join("logs","rules","aggregate_all_clustering_results_{sample}.log"), + os.path.join("logs", "rules", "aggregate_all_clustering_results_{sample}.log"), run: # list to hold the data agg_clust = [] @@ -81,9 +115,9 @@ rule aggregate_all_clustering_results: # read each clustering result and add to data dict for filename in input: agg_clust.append(pd.read_csv(filename, header=0, index_col=0)) - - # convert the dictionary to a DataFrame + + # convert the dictionary to a DataFrame agg_clust_df = pd.concat(agg_clust, axis=1) - + # Write the DataFrame to a CSV file agg_clust_df.to_csv(output.metadata_clusterings, index=True) diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 477f86c..47e340b 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -2,150 +2,348 @@ ########## GENERAL ########## + def get_sample_paths(wildcards): - return {'data': annot.loc[wildcards.sample,'data'], - 'metadata': annot.loc[wildcards.sample,"metadata"] - } + return { + "data": annot.loc[wildcards.sample, "data"], + "metadata": annot.loc[wildcards.sample, "metadata"], + } + def get_data_orientation(wildcards): - return int(annot.loc[wildcards.sample,'samples_by_features']) + return int(annot.loc[wildcards.sample, "samples_by_features"]) + ########## DIMENSIONALITY REDUCTION ########## + def get_umap_sample_paths(wildcards): - return [annot.loc[wildcards.sample,'data'], - os.path.join(result_path,'{}'.format(wildcards.sample),'UMAP','UMAP_{}_'.format(wildcards.metric)+'{}'.format(max(config["umap"]["n_neighbors"]))+'_graph.pickle')] + return [ + annot.loc[wildcards.sample, "data"], + os.path.join( + result_path, + "{}".format(wildcards.sample), + "UMAP", + "UMAP_{}_".format(wildcards.metric) + + "{}".format(max(config["umap"]["n_neighbors"])) + + "_graph.pickle", + ), + ] + def get_dimred_paths(wildcards): path_dict = {} - - if wildcards.method=="PCA": - path_dict['dimred_data'] = os.path.join(result_path,wildcards.sample,'PCA','PCA_{wildcards.parameters}_data_small.csv'.format(wildcards=wildcards)) - path_dict['dimred_axes'] = os.path.join(result_path,wildcards.sample,'PCA','PCA_{wildcards.parameters}_axes.csv'.format(wildcards=wildcards)) - path_dict['dimred_var'] = os.path.join(result_path,wildcards.sample,'PCA','PCA_{wildcards.parameters}_var.csv'.format(wildcards=wildcards)) - path_dict['dimred_loadings'] = os.path.join(result_path,wildcards.sample,'PCA','PCA_{wildcards.parameters}_loadings_small.csv'.format(wildcards=wildcards)) + + if wildcards.method == "PCA": + path_dict["dimred_data"] = os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{wildcards.parameters}_data_small.csv".format(wildcards=wildcards), + ) + path_dict["dimred_axes"] = os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{wildcards.parameters}_axes.csv".format(wildcards=wildcards), + ) + path_dict["dimred_var"] = os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{wildcards.parameters}_var.csv".format(wildcards=wildcards), + ) + path_dict["dimred_loadings"] = os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{wildcards.parameters}_loadings_small.csv".format(wildcards=wildcards), + ) else: - path_dict['dimred_data'] = os.path.join(result_path,wildcards.sample,wildcards.method,'{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_data.csv'.format(wildcards=wildcards)) - path_dict['dimred_axes'] = os.path.join(result_path,wildcards.sample,wildcards.method,'{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_axes.csv'.format(wildcards=wildcards)) - + path_dict["dimred_data"] = os.path.join( + result_path, + wildcards.sample, + wildcards.method, + "{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_data.csv".format( + wildcards=wildcards + ), + ) + path_dict["dimred_axes"] = os.path.join( + result_path, + wildcards.sample, + wildcards.method, + "{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_axes.csv".format( + wildcards=wildcards + ), + ) + # add metadata - path_dict['metadata'] = annot.loc[wildcards.sample,"metadata"] + path_dict["metadata"] = annot.loc[wildcards.sample, "metadata"] # add features - path_dict['metadata_features'] = os.path.join(result_path,wildcards.sample,'metadata_features.csv') + path_dict["metadata_features"] = os.path.join( + result_path, wildcards.sample, "metadata_features.csv" + ) # add clustering results if len(cluster_methods) > 0: - path_dict['metadata_clusterings'] = os.path.join(result_path,wildcards.sample,'metadata_clusterings.csv') - + path_dict["metadata_clusterings"] = os.path.join( + result_path, wildcards.sample, "metadata_clusterings.csv" + ) + return path_dict - + + def get_dimred_features_paths(wildcards): - - if wildcards.method=="PCA": + + if wildcards.method == "PCA": return { - 'dimred_data': os.path.join(result_path,wildcards.sample,'PCA','PCA_{wildcards.parameters}_data_small.csv'.format(wildcards=wildcards)), - 'dimred_axes': os.path.join(result_path,wildcards.sample,'PCA','PCA_{wildcards.parameters}_axes.csv'.format(wildcards=wildcards)), - 'metadata': os.path.join(result_path,wildcards.sample,'metadata_features.csv') - } + "dimred_data": os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{wildcards.parameters}_data_small.csv".format(wildcards=wildcards), + ), + "dimred_axes": os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{wildcards.parameters}_axes.csv".format(wildcards=wildcards), + ), + "metadata": os.path.join( + result_path, wildcards.sample, "metadata_features.csv" + ), + } else: return { - 'dimred_data': os.path.join(result_path,wildcards.sample,wildcards.method,'{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_data.csv'.format(wildcards=wildcards)), - 'dimred_axes': os.path.join(result_path,wildcards.sample,wildcards.method,'{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_axes.csv'.format(wildcards=wildcards)), - 'metadata': os.path.join(result_path,wildcards.sample,'metadata_features.csv') - } + "dimred_data": os.path.join( + result_path, + wildcards.sample, + wildcards.method, + "{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_data.csv".format( + wildcards=wildcards + ), + ), + "dimred_axes": os.path.join( + result_path, + wildcards.sample, + wildcards.method, + "{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_axes.csv".format( + wildcards=wildcards + ), + ), + "metadata": os.path.join( + result_path, wildcards.sample, "metadata_features.csv" + ), + } + ########## HEATMAPS ########## def get_heatmap_paths(wildcards): - return {'data': annot.loc[wildcards.sample,'data'], - 'metadata': annot.loc[wildcards.sample,"metadata"], - 'observations_distance': os.path.join(result_path,wildcards.sample,'Heatmap','DistanceMatrix_{wildcards.metric}_observations.csv'.format(wildcards=wildcards)), - 'features_distance': os.path.join(result_path,wildcards.sample,'Heatmap','DistanceMatrix_{wildcards.metric}_features.csv'.format(wildcards=wildcards)), - } + return { + "data": annot.loc[wildcards.sample, "data"], + "metadata": annot.loc[wildcards.sample, "metadata"], + "observations_distance": os.path.join( + result_path, + wildcards.sample, + "Heatmap", + "DistanceMatrix_{wildcards.metric}_observations.csv".format( + wildcards=wildcards + ), + ), + "features_distance": os.path.join( + result_path, + wildcards.sample, + "Heatmap", + "DistanceMatrix_{wildcards.metric}_features.csv".format( + wildcards=wildcards + ), + ), + } + ########## CLUSTERING ########## + # get paths for clustification def get_clustification_paths(wildcards): - return [annot.loc[wildcards.sample,'data'], - os.path.join(result_path,wildcards.sample,'{}'.format(config["clustification"]["method"]),'{}_clusterings.csv'.format(config["clustification"]["method"])) - ] + return [ + annot.loc[wildcards.sample, "data"], + os.path.join( + result_path, + wildcards.sample, + "{}".format(config["clustification"]["method"]), + "{}_clusterings.csv".format(config["clustification"]["method"]), + ), + ] + # get all clustering results of one method to be aggregated into {method}/{method}_clusterings.csv def get_clustering_paths(wildcards): - path_list = [annot.loc[wildcards.sample,"metadata"]] - - if wildcards.method=="Leiden": + path_list = [annot.loc[wildcards.sample, "metadata"]] + + if wildcards.method == "Leiden": leiden_parameters = [] - + # add resolution parameter only to relevant partition_type algorithms, otherwise NA for partition_type in config["leiden"]["partition_types"]: - if partition_type in ["RBConfigurationVertexPartition", "RBERVertexPartition", "CPMVertexPartition"]: - leiden_parameters = leiden_parameters + ["{}_{}".format(partition_type, res) for res in config["leiden"]["resolutions"]] + if partition_type in [ + "RBConfigurationVertexPartition", + "RBERVertexPartition", + "CPMVertexPartition", + ]: + leiden_parameters = leiden_parameters + [ + "{}_{}".format(partition_type, res) + for res in config["leiden"]["resolutions"] + ] else: leiden_parameters.append("{}_NA".format(partition_type)) - - path_list = path_list + expand(os.path.join(result_path,wildcards.sample,'Leiden','Leiden_{metric}_{n_neighbors}_{leiden_parameters}_clustering.csv'), - metric=config["leiden"]["metrics"], - n_neighbors=config["leiden"]["n_neighbors"], - leiden_parameters=leiden_parameters, - ) + + path_list = path_list + expand( + os.path.join( + result_path, + wildcards.sample, + "Leiden", + "Leiden_{metric}_{n_neighbors}_{leiden_parameters}_clustering.csv", + ), + metric=config["leiden"]["metrics"], + n_neighbors=config["leiden"]["n_neighbors"], + leiden_parameters=leiden_parameters, + ) return path_list + # get all aggregated clustering results across methods to be aggregated into {sample}/metadata_clusterings.csv def get_aggregated_clustering_paths(wildcards): - return expand(os.path.join(result_path,wildcards.sample,'{method}','{method}_clusterings.csv'), method=cluster_methods) + return expand( + os.path.join( + result_path, wildcards.sample, "{method}", "{method}_clusterings.csv" + ), + method=cluster_methods, + ) + # get the aggregated clustering results across methods for visualization def get_metadata_clustering_paths(wildcards): - - if wildcards.method=="PCA": + + if wildcards.method == "PCA": return { - 'dimred_data': os.path.join(result_path,wildcards.sample,'PCA','PCA_{wildcards.parameters}_data_small.csv'.format(wildcards=wildcards)), - 'dimred_axes': os.path.join(result_path,wildcards.sample,'PCA','PCA_{wildcards.parameters}_axes.csv'.format(wildcards=wildcards)), - 'metadata': os.path.join(result_path,wildcards.sample,'metadata_clusterings.csv') - } + "dimred_data": os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{wildcards.parameters}_data_small.csv".format(wildcards=wildcards), + ), + "dimred_axes": os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{wildcards.parameters}_axes.csv".format(wildcards=wildcards), + ), + "metadata": os.path.join( + result_path, wildcards.sample, "metadata_clusterings.csv" + ), + } else: return { - 'dimred_data': os.path.join(result_path,wildcards.sample,wildcards.method,'{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_data.csv'.format(wildcards=wildcards)), - 'dimred_axes': os.path.join(result_path,wildcards.sample,wildcards.method,'{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_axes.csv'.format(wildcards=wildcards)), - 'metadata': os.path.join(result_path,wildcards.sample,'metadata_clusterings.csv') - } + "dimred_data": os.path.join( + result_path, + wildcards.sample, + wildcards.method, + "{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_data.csv".format( + wildcards=wildcards + ), + ), + "dimred_axes": os.path.join( + result_path, + wildcards.sample, + wildcards.method, + "{wildcards.method}_{wildcards.parameters}_{wildcards.n_components}_axes.csv".format( + wildcards=wildcards + ), + ), + "metadata": os.path.join( + result_path, wildcards.sample, "metadata_clusterings.csv" + ), + } + ########## CLUSTER VALIDATION ########## + # get the input paths for clustree analysis depending on requested content type def get_clustree_paths(wildcards): - - if wildcards.content=="features": + + if wildcards.content == "features": return { - 'metadata_clustering': os.path.join(result_path,wildcards.sample, "metadata_clusterings.csv"), - 'metadata': os.path.join(result_path,wildcards.sample,'metadata_features.csv') + "metadata_clustering": os.path.join( + result_path, wildcards.sample, "metadata_clusterings.csv" + ), + "metadata": os.path.join( + result_path, wildcards.sample, "metadata_features.csv" + ), } else: return { - 'metadata_clustering': os.path.join(result_path,wildcards.sample, "metadata_clusterings.csv"), - 'metadata': annot.loc[wildcards.sample,"metadata"] + "metadata_clustering": os.path.join( + result_path, wildcards.sample, "metadata_clusterings.csv" + ), + "metadata": annot.loc[wildcards.sample, "metadata"], } + # get paths to determine external cluster indices def get_external_validation_paths(wildcards): - return {'clusterings': os.path.join(result_path,wildcards.sample, "metadata_clusterings.csv"), - 'metadata': annot.loc[wildcards.sample,"metadata"] - } + return { + "clusterings": os.path.join( + result_path, wildcards.sample, "metadata_clusterings.csv" + ), + "metadata": annot.loc[wildcards.sample, "metadata"], + } + # get paths to determine internal cluster indices def get_internal_validation_paths(wildcards): - return {'metadata': annot.loc[wildcards.sample,"metadata"], - 'clusterings': os.path.join(result_path,wildcards.sample, "metadata_clusterings.csv"), - 'pca': os.path.join(result_path,wildcards.sample,'PCA','PCA_{}_{}_data.csv'.format(config["pca"]["svd_solver"],config["pca"]["n_components"])), - 'pca_var': os.path.join(result_path,wildcards.sample,'PCA','PCA_{}_{}_var.csv'.format(config["pca"]["svd_solver"],config["pca"]["n_components"])) - } + return { + "metadata": annot.loc[wildcards.sample, "metadata"], + "clusterings": os.path.join( + result_path, wildcards.sample, "metadata_clusterings.csv" + ), + "pca": os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{}_{}_data.csv".format( + config["pca"]["svd_solver"], config["pca"]["n_components"] + ), + ), + "pca_var": os.path.join( + result_path, + wildcards.sample, + "PCA", + "PCA_{}_{}_var.csv".format( + config["pca"]["svd_solver"], config["pca"]["n_components"] + ), + ), + } + # for plotting heatmaps of cluster indices def get_validation_paths(wildcards): - if wildcards.type=="external": + if wildcards.type == "external": return { - idx: os.path.join(result_path,wildcards.sample, "cluster_validation", "external_index_{}.csv".format(idx)) for idx in indices_external + idx: os.path.join( + result_path, + wildcards.sample, + "cluster_validation", + "external_index_{}.csv".format(idx), + ) + for idx in indices_external } else: - return {"ranked_internal_indices": os.path.join(result_path,wildcards.sample, "cluster_validation", "internal_indices_ranked.csv")} - + return { + "ranked_internal_indices": os.path.join( + result_path, + wildcards.sample, + "cluster_validation", + "internal_indices_ranked.csv", + ) + } diff --git a/workflow/rules/dimred.smk b/workflow/rules/dimred.smk index bf90b63..bf185f0 100644 --- a/workflow/rules/dimred.smk +++ b/workflow/rules/dimred.smk @@ -3,114 +3,174 @@ rule pca: input: unpack(get_sample_paths), output: - result_object = os.path.join(result_path,'{sample}','PCA','PCA_{parameters}_object.pickle'), - result_data = os.path.join(result_path,'{sample}','PCA','PCA_{parameters}_data.csv'), - result_data_small = os.path.join(result_path,'{sample}','PCA','PCA_{parameters}_data_small.csv'), - result_loadings = os.path.join(result_path,'{sample}','PCA','PCA_{parameters}_loadings.csv'), - result_loadings_small = os.path.join(result_path,'{sample}','PCA','PCA_{parameters}_loadings_small.csv'), - result_var = os.path.join(result_path,'{sample}','PCA','PCA_{parameters}_var.csv'), - result_axes = os.path.join(result_path,'{sample}','PCA','PCA_{parameters}_axes.csv'), + result_object=os.path.join( + result_path, "{sample}", "PCA", "PCA_{parameters}_object.pickle" + ), + result_data=os.path.join( + result_path, "{sample}", "PCA", "PCA_{parameters}_data.csv" + ), + result_data_small=os.path.join( + result_path, "{sample}", "PCA", "PCA_{parameters}_data_small.csv" + ), + result_loadings=os.path.join( + result_path, "{sample}", "PCA", "PCA_{parameters}_loadings.csv" + ), + result_loadings_small=os.path.join( + result_path, "{sample}", "PCA", "PCA_{parameters}_loadings_small.csv" + ), + result_var=os.path.join( + result_path, "{sample}", "PCA", "PCA_{parameters}_var.csv" + ), + result_axes=os.path.join( + result_path, "{sample}", "PCA", "PCA_{parameters}_axes.csv" + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","PCA_{sample}_{parameters}.log"), + os.path.join("logs", "rules", "PCA_{sample}_{parameters}.log"), params: - samples_by_features = get_data_orientation, + samples_by_features=get_data_orientation, script: "../scripts/pca.py" - - + + ####### perform Uniform Manifold Approximation and Projection (UMAP) ####### + # generate parametrized knn graphs using the UMAP package rule umap_graph: input: unpack(get_sample_paths), output: - result_object = os.path.join(result_path,'{sample}','UMAP','UMAP_{metric}_{n_neighbors}_graph.pickle'), + result_object=os.path.join( + result_path, "{sample}", "UMAP", "UMAP_{metric}_{n_neighbors}_graph.pickle" + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","umap_{sample}_{metric}_{n_neighbors}.log"), + os.path.join("logs", "rules", "umap_{sample}_{metric}_{n_neighbors}.log"), params: - samples_by_features = get_data_orientation, - metric = lambda w: "{}".format(w.metric), - n_neighbors = lambda w: "{}".format(w.n_neighbors), + samples_by_features=get_data_orientation, + metric=lambda w: "{}".format(w.metric), + n_neighbors=lambda w: "{}".format(w.n_neighbors), script: "../scripts/umap_graph.py" - - + + # embed parametrized knn graphs using the UMAP package into lower dimensional space rule umap_embed: input: get_umap_sample_paths, output: - result_object = os.path.join(result_path,'{sample}','UMAP','UMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_object.pickle'), - result_data = os.path.join(result_path,'{sample}','UMAP','UMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_data.csv'), - result_axes = os.path.join(result_path,'{sample}','UMAP','UMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_axes.csv'), + result_object=os.path.join( + result_path, + "{sample}", + "UMAP", + "UMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_object.pickle", + ), + result_data=os.path.join( + result_path, + "{sample}", + "UMAP", + "UMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_data.csv", + ), + result_axes=os.path.join( + result_path, + "{sample}", + "UMAP", + "UMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_axes.csv", + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","umap_{sample}_{metric}_{n_neighbors}_{min_dist}_{n_components}.log"), + os.path.join( + "logs", + "rules", + "umap_{sample}_{metric}_{n_neighbors}_{min_dist}_{n_components}.log", + ), params: - samples_by_features = get_data_orientation, - metric = lambda w: "{}".format(w.metric), - n_neighbors = lambda w: "{}".format(w.n_neighbors), - min_dist = lambda w: "{}".format(w.min_dist), - n_components = lambda w: "{}".format(w.n_components), - densmap = 0, + samples_by_features=get_data_orientation, + metric=lambda w: "{}".format(w.metric), + n_neighbors=lambda w: "{}".format(w.n_neighbors), + min_dist=lambda w: "{}".format(w.min_dist), + n_components=lambda w: "{}".format(w.n_components), + densmap=0, script: "../scripts/umap_embed.py" - + + # embed parametrized knn graphs using the UMAP package with densMAP flag into lower dimensional space rule densmap_embed: input: get_umap_sample_paths, output: - result_object = os.path.join(result_path,'{sample}','densMAP','densMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_object.pickle'), - result_data = os.path.join(result_path,'{sample}','densMAP','densMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_data.csv'), - result_axes = os.path.join(result_path,'{sample}','densMAP','densMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_axes.csv'), + result_object=os.path.join( + result_path, + "{sample}", + "densMAP", + "densMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_object.pickle", + ), + result_data=os.path.join( + result_path, + "{sample}", + "densMAP", + "densMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_data.csv", + ), + result_axes=os.path.join( + result_path, + "{sample}", + "densMAP", + "densMAP_{metric}_{n_neighbors}_{min_dist}_{n_components}_axes.csv", + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","densmap_{sample}_{metric}_{n_neighbors}_{min_dist}_{n_components}.log"), + os.path.join( + "logs", + "rules", + "densmap_{sample}_{metric}_{n_neighbors}_{min_dist}_{n_components}.log", + ), params: - samples_by_features = get_data_orientation, - metric = lambda w: "{}".format(w.metric), - n_neighbors = lambda w: "{}".format(w.n_neighbors), - min_dist = lambda w: "{}".format(w.min_dist), - n_components = lambda w: "{}".format(w.n_components), - densmap = 1, + samples_by_features=get_data_orientation, + metric=lambda w: "{}".format(w.metric), + n_neighbors=lambda w: "{}".format(w.n_neighbors), + min_dist=lambda w: "{}".format(w.min_dist), + n_components=lambda w: "{}".format(w.n_components), + densmap=1, script: "../scripts/umap_embed.py" + ####### determine distance matrices ####### + rule distance_matrix: input: unpack(get_sample_paths), output: - distance_matrix = os.path.join(result_path,'{sample}','Heatmap','DistanceMatrix_{metric}_{type}.csv'), + distance_matrix=os.path.join( + result_path, "{sample}", "Heatmap", "DistanceMatrix_{metric}_{type}.csv" + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","DistanceMatrix_{sample}_{metric}_{type}.log"), + os.path.join("logs", "rules", "DistanceMatrix_{sample}_{metric}_{type}.log"), params: - samples_by_features = get_data_orientation, + samples_by_features=get_data_orientation, script: - "../scripts/distance_matrix.py" \ No newline at end of file + "../scripts/distance_matrix.py" diff --git a/workflow/rules/envs_export.smk b/workflow/rules/envs_export.smk index b6ff371..abecf85 100644 --- a/workflow/rules/envs_export.smk +++ b/workflow/rules/envs_export.smk @@ -1,55 +1,64 @@ # conda environment export rule to document the exact versions and builds of the used software rule env_export: output: - report(os.path.join(result_path,'envs','{env}.yaml'), - caption="../report/software.rst", - category="Software", - subcategory="{}_{}".format(config["project_name"], module_name) - ), + report( + os.path.join(result_path, "envs", "{env}.yaml"), + caption="../report/software.rst", + category="Software", + subcategory="{}_{}".format(config["project_name"], module_name), + ), conda: "../envs/{env}.yaml" resources: mem_mb=1000, threads: config.get("threads", 1) log: - os.path.join("logs","rules","env_{env}.log"), + os.path.join("logs", "rules", "env_{env}.log"), shell: """ conda env export > {output} """ - + + # export used config file for documentation and reproducibility rule config_export: output: - configs = report(os.path.join(result_path,'configs','{}_config.yaml'.format(config["project_name"])), - caption="../report/configs.rst", - category="Configuration", - subcategory="{}_{}".format(config["project_name"], module_name) - ) + configs=report( + os.path.join( + result_path, "configs", "{}_config.yaml".format(config["project_name"]) + ), + caption="../report/configs.rst", + category="Configuration", + subcategory="{}_{}".format(config["project_name"], module_name), + ), resources: mem_mb=1000, threads: config.get("threads", 1) log: - os.path.join("logs","rules","config_export.log"), + os.path.join("logs", "rules", "config_export.log"), run: - with open(output["configs"], 'w') as outfile: + with open(output["configs"], "w") as outfile: yaml.dump(config, outfile) -# export used annotation file for documentation and reproducibility + +# export used annotation file for documentation and reproducibility rule annot_export: input: config["annotation"], output: - annot = report(os.path.join(result_path,'configs','{}_annot.csv'.format(config["project_name"])), - caption="../report/configs.rst", - category="Configuration", - subcategory="{}_{}".format(config["project_name"], module_name) - ) + annot=report( + os.path.join( + result_path, "configs", "{}_annot.csv".format(config["project_name"]) + ), + caption="../report/configs.rst", + category="Configuration", + subcategory="{}_{}".format(config["project_name"], module_name), + ), resources: mem_mb=1000, threads: config.get("threads", 1) log: - os.path.join("logs","rules","annot_export.log"), + os.path.join("logs", "rules", "annot_export.log"), shell: """ cp {input} {output} diff --git a/workflow/rules/visualization.smk b/workflow/rules/visualization.smk index 72c01e2..95a6297 100644 --- a/workflow/rules/visualization.smk +++ b/workflow/rules/visualization.smk @@ -1,222 +1,339 @@ ########## FEATURE PLOTS ########## + # prepare "metadata" for feature plot rule prep_feature_plot: input: unpack(get_sample_paths), output: - metadata_features = os.path.join(result_path,'{sample}','metadata_features.csv'), + metadata_features=os.path.join(result_path, "{sample}", "metadata_features.csv"), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","prep_feature_plot_{sample}.log"), + os.path.join("logs", "rules", "prep_feature_plot_{sample}.log"), params: - samples_by_features = get_data_orientation, - features_to_plot = config["features_to_plot"], + samples_by_features=get_data_orientation, + features_to_plot=config["features_to_plot"], script: "../scripts/subset_data.py" - + + # dimred feature scatter plots rule plot_dimred_features: input: unpack(get_dimred_features_paths), output: - plot = report(directory(os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}_{n_components}','features')), - patterns=["{feature}.png"], - caption="../report/dimred_2d_features.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{parameters}", - "dimensions": "{n_components}", - "type": "features", - "content": "{feature}", - } - ), + plot=report( + directory( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}_{n_components}", + "features", + ) + ), + patterns=["{feature}.png"], + caption="../report/dimred_2d_features.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{parameters}", + "dimensions": "{n_components}", + "type": "features", + "content": "{feature}", + }, + ), resources: - mem_mb=config.get("mem", "32000"), + mem_mb=config.get("mem", 32000), threads: config.get("threads", 1) conda: "../envs/ggplot.yaml" log: - os.path.join("logs","rules","plot_features_{sample}_{method}_{parameters}_{n_components}.log"), + os.path.join( + "logs", + "rules", + "plot_features_{sample}_{method}_{parameters}_{n_components}.log", + ), params: - size = config["scatterplot2d"]["size"], - alpha = config["scatterplot2d"]["alpha"], + size=config["scatterplot2d"]["size"], + alpha=config["scatterplot2d"]["alpha"], script: "../scripts/plot_2d.R" + ########## METADATA PLOTS ########## - + + # dimred scatter plot panel by metadata rule plot_dimred_metadata: input: unpack(get_dimred_paths), output: - plot = report(directory(os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}_{n_components}','metadata')), - patterns=["{metadata}.png"], - caption="../report/dimred_2d_metadata.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{parameters}", - "dimensions": "{n_components}", - "type": "metadata", - "content": "{metadata}", - } - ), + plot=report( + directory( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}_{n_components}", + "metadata", + ) + ), + patterns=["{metadata}.png"], + caption="../report/dimred_2d_metadata.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{parameters}", + "dimensions": "{n_components}", + "type": "metadata", + "content": "{metadata}", + }, + ), resources: - mem_mb=config.get("mem", "32000"), + mem_mb=config.get("mem", 32000), threads: config.get("threads", 1) conda: "../envs/ggplot.yaml" log: - os.path.join("logs","rules","plot_metadata_{sample}_{method}_{parameters}_{n_components}.log"), + os.path.join( + "logs", + "rules", + "plot_metadata_{sample}_{method}_{parameters}_{n_components}.log", + ), params: - size = config["scatterplot2d"]["size"], - alpha = config["scatterplot2d"]["alpha"], + size=config["scatterplot2d"]["size"], + alpha=config["scatterplot2d"]["alpha"], script: "../scripts/plot_2d.R" + ########## DIAGNOSTIC PLOTS ########## - + + # PCA scree plot, cumulative variance plot, pairs plot, and loadings plot rule plot_pca_diagnostics: input: - unpack(get_dimred_paths) + unpack(get_dimred_paths), output: - diagnostics_plot = report(os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}','diagnostics','variance.png'), - caption="../report/pca_diagnostics.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{parameters}", - "dimensions": "-", - "type": "diagnostics", - "content": "variance", - }), - pairs_plot = report(os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}','diagnostics','pairs.png'), - caption="../report/pca_diagnostics.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{parameters}", - "dimensions": "-", - "type": "diagnostics", - "content": "pairs", - }), - loadings_plot = report(os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}','diagnostics','loadings.png'), - caption="../report/pca_diagnostics.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{parameters}", - "dimensions": "-", - "type": "diagnostics", - "content": "loadings", - }), - loadings_lollipop_plot = report(os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}','diagnostics','loadings_lollipop.png'), - caption="../report/pca_diagnostics.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{parameters}", - "dimensions": "-", - "type": "diagnostics", - "content": "loadings lollipop", - }), + diagnostics_plot=report( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}", + "diagnostics", + "variance.png", + ), + caption="../report/pca_diagnostics.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{parameters}", + "dimensions": "-", + "type": "diagnostics", + "content": "variance", + }, + ), + pairs_plot=report( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}", + "diagnostics", + "pairs.png", + ), + caption="../report/pca_diagnostics.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{parameters}", + "dimensions": "-", + "type": "diagnostics", + "content": "pairs", + }, + ), + loadings_plot=report( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}", + "diagnostics", + "loadings.png", + ), + caption="../report/pca_diagnostics.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{parameters}", + "dimensions": "-", + "type": "diagnostics", + "content": "loadings", + }, + ), + loadings_lollipop_plot=report( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}", + "diagnostics", + "loadings_lollipop.png", + ), + caption="../report/pca_diagnostics.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{parameters}", + "dimensions": "-", + "type": "diagnostics", + "content": "loadings lollipop", + }, + ), resources: mem_mb=config.get("mem", "8000"), threads: config.get("threads", 1) conda: "../envs/ggplot.yaml" log: - os.path.join("logs","rules","plot_{method}_diagnostics_{sample}_{parameters}.log"), + os.path.join( + "logs", "rules", "plot_{method}_diagnostics_{sample}_{parameters}.log" + ), script: "../scripts/plot_pca.R" - + # plot UMAP & densMAP diagnostic visualizations rule plot_umap_diagnostics: input: - umap_object = os.path.join(result_path,'{sample}','{method}','{method}_{parameters}_object.pickle'), + umap_object=os.path.join( + result_path, "{sample}", "{method}", "{method}_{parameters}_object.pickle" + ), output: - diagnostics_plot = report(os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}','diagnostics','diagnostics.png'), - caption="../report/umap_diagnostics.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{parameters}", - "dimensions": "-", - "type": "diagnostics", - "content": "diagnostics", - }), + diagnostics_plot=report( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}", + "diagnostics", + "diagnostics.png", + ), + caption="../report/umap_diagnostics.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{parameters}", + "dimensions": "-", + "type": "diagnostics", + "content": "diagnostics", + }, + ), resources: - mem_mb=config.get("mem", "32000"), + mem_mb=config.get("mem", 32000), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","plot_diagnostics_{sample}_{method}_{parameters}.log"), + os.path.join( + "logs", "rules", "plot_diagnostics_{sample}_{method}_{parameters}.log" + ), script: "../scripts/plot_umap_diagnostics.py" - - + + # plot UMAP & densMAP connectivity visualizations rule plot_umap_connectivity: input: - umap_object = os.path.join(result_path,'{sample}','{method}','{method}_{parameters}_object.pickle'), + umap_object=os.path.join( + result_path, "{sample}", "{method}", "{method}_{parameters}_object.pickle" + ), output: - connectivity_plot = report(os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}','diagnostics','connectivity.png'), - caption="../report/umap_connectivity.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{parameters}", - "dimensions": "-", - "type": "diagnostics", - "content": "connectivity", - }), + connectivity_plot=report( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}", + "diagnostics", + "connectivity.png", + ), + caption="../report/umap_connectivity.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{parameters}", + "dimensions": "-", + "type": "diagnostics", + "content": "connectivity", + }, + ), resources: mem_mb=config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/umap_leiden.yaml" log: - os.path.join("logs","rules","plot_connectivity_{sample}_{method}_{parameters}.log"), + os.path.join( + "logs", "rules", "plot_connectivity_{sample}_{method}_{parameters}.log" + ), script: "../scripts/plot_umap_connectivity.py" + ########## INTERACTIVE PLOTS ########## + # plot interactive HTML plots of 2D and 3D embeddings using plotly rule plot_dimred_interactive: input: unpack(get_dimred_paths), output: - plot = os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}_{n_components}','interactive.html'), + plot=os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}_{n_components}", + "interactive.html", + ), resources: mem_mb=config.get("mem", "8000"), threads: config.get("threads", 1) conda: "../envs/plotly.yaml" log: - os.path.join("logs","rules","plot_interactive_{sample}_{method}_{parameters}_{n_components}.log"), + os.path.join( + "logs", + "rules", + "plot_interactive_{sample}_{method}_{parameters}_{n_components}.log", + ), params: - n_components = lambda w: "{}".format(w.n_components), - size = config["scatterplot2d"]["size"], - alpha = config["scatterplot2d"]["alpha"] + n_components=lambda w: "{}".format(w.n_components), + size=config["scatterplot2d"]["size"], + alpha=config["scatterplot2d"]["alpha"], script: "../scripts/plot_interactive.py" @@ -226,60 +343,83 @@ rule plot_heatmap: input: unpack(get_heatmap_paths), output: - plot = report(os.path.join(result_path,'{sample}','Heatmap','plots','Heatmap_{metric}_{method}.png'), - caption="../report/heatmap.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{metric}", - "dimensions": "-", - "type": "Heatmap", - "content": "Heatmap", - }), + plot=report( + os.path.join( + result_path, + "{sample}", + "Heatmap", + "plots", + "Heatmap_{metric}_{method}.png", + ), + caption="../report/heatmap.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{metric}", + "dimensions": "-", + "type": "Heatmap", + "content": "Heatmap", + }, + ), resources: # dynamic memory allocation based on input size and attempts (multiple attempts can be triggered with --retries X) - mem_mb=lambda wildcards, attempt: attempt*int(config.get("mem", "16000")),#lambda wildcards, input, attempt: max(int(config.get("mem", "16000")),((input.size//1000000) * attempt * 70)),#config.get("mem", "16000"), + mem_mb=lambda wildcards, attempt: attempt * int(config.get("mem", "16000")), #lambda wildcards, input, attempt: max(int(config.get("mem", "16000")),((input.size//1000000) * attempt * 70)),#config.get("mem", "16000"), threads: config.get("threads", 1) conda: "../envs/ComplexHeatmap.yaml" log: - os.path.join("logs","rules","plot_heatmap_{sample}_{metric}_{method}.log"), + os.path.join("logs", "rules", "plot_heatmap_{sample}_{metric}_{method}.log"), params: - samples_by_features = get_data_orientation, + samples_by_features=get_data_orientation, script: "../scripts/plot_heatmap.R" - - + + ########## CLUSTERING PLOTS ########## - + + # dimred 2D scatter plot panel by clustering rule plot_dimred_clustering: input: unpack(get_metadata_clustering_paths), output: - plot = report(directory(os.path.join(result_path,'{sample}','{method}','plots','{method}_{parameters}_{n_components}','clustering')), - patterns=["{clustering}.png"], - caption="../report/dimred_2d_clusterings.rst", - category="{}_{}".format(config["project_name"], module_name), - subcategory="{sample}", - labels={ - "method": "{method}", - "parameters": "{parameters}", - "dimensions": "{n_components}", - "type": "clustering", - "content": "{clustering}", - } - ), + plot=report( + directory( + os.path.join( + result_path, + "{sample}", + "{method}", + "plots", + "{method}_{parameters}_{n_components}", + "clustering", + ) + ), + patterns=["{clustering}.png"], + caption="../report/dimred_2d_clusterings.rst", + category="{}_{}".format(config["project_name"], module_name), + subcategory="{sample}", + labels={ + "method": "{method}", + "parameters": "{parameters}", + "dimensions": "{n_components}", + "type": "clustering", + "content": "{clustering}", + }, + ), resources: - mem_mb=config.get("mem", "32000"), + mem_mb=config.get("mem", 32000), threads: config.get("threads", 1) conda: "../envs/ggplot.yaml" log: - os.path.join("logs","rules","plot_clustering_{sample}_{method}_{parameters}_{n_components}.log"), + os.path.join( + "logs", + "rules", + "plot_clustering_{sample}_{method}_{parameters}_{n_components}.log", + ), params: - size = config["scatterplot2d"]["size"], - alpha = config["scatterplot2d"]["alpha"] + size=config["scatterplot2d"]["size"], + alpha=config["scatterplot2d"]["alpha"], script: "../scripts/plot_2d.R" diff --git a/workflow/scripts/clustification.py b/workflow/scripts/clustification.py index 64093d0..78ef843 100644 --- a/workflow/scripts/clustification.py +++ b/workflow/scripts/clustification.py @@ -7,8 +7,10 @@ from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import StratifiedKFold from sklearn.metrics import confusion_matrix + # from sklearn.metrics import f1_score, accuracy_score + # clustification function def iterative_classification(data, labels, n_trees=100, max_iterations=100): @@ -18,32 +20,36 @@ def iterative_classification(data, labels, n_trees=100, max_iterations=100): new_labels = np.zeros_like(labels) new_data = np.zeros_like(data) skf = StratifiedKFold(n_splits=5, shuffle=True) - + for train_index, test_index in skf.split(data, labels): X_train, X_test = data[train_index], data[test_index] y_train, y_test = labels[train_index], labels[test_index] - clf = RandomForestClassifier(n_estimators=n_trees, - random_state = 42, - n_jobs = -1) + clf = RandomForestClassifier( + n_estimators=n_trees, random_state=42, n_jobs=-1 + ) clf.fit(X_train, y_train) y_pred = clf.predict(X_test) new_labels[test_index] = y_pred new_data[test_index] = X_test # determine symmetric confusiton matrix interpreted as weighted graph - cm = confusion_matrix(y_true=labels, y_pred=new_labels, normalize='true') + cm = confusion_matrix(y_true=labels, y_pred=new_labels, normalize="true") cm_symmetric = (cm + cm.T) / 2 cm_symmetric = np.triu(cm_symmetric, 1) max_weight = np.max(cm_symmetric) print("max edge weight: {}".format(max_weight)) - + # check stopping criterion if max_weight < 0.025: print("Misclassification threshold met, stopping.") break - - merge_indices = np.unravel_index(np.argmax(cm_symmetric, axis=None), cm_symmetric.shape) - labels = np.array([merge_indices[0] if num == merge_indices[1] else num for num in labels]) + + merge_indices = np.unravel_index( + np.argmax(cm_symmetric, axis=None), cm_symmetric.shape + ) + labels = np.array( + [merge_indices[0] if num == merge_indices[1] else num for num in labels] + ) labels, _ = pd.factorize(labels) data = new_data print("#clusters: {}".format(len(set(labels)))) @@ -57,7 +63,7 @@ def iterative_classification(data, labels, n_trees=100, max_iterations=100): clusterings_path = os.path.join(snakemake.input[1]) # parameters -samples_by_features = int(snakemake.params['samples_by_features']) #0 +samples_by_features = int(snakemake.params["samples_by_features"]) # 0 # outputs result_path = os.path.join(snakemake.output["clustering"]) @@ -69,12 +75,16 @@ def iterative_classification(data, labels, n_trees=100, max_iterations=100): else: data = pd.read_csv(data_path, index_col=0).T - + clusterings = pd.read_csv(clusterings_path, index_col=0) clustering_init = clusterings[clusterings.nunique().idxmax()].to_numpy().ravel() # run clustification -clustering_new = iterative_classification(data.to_numpy(), clustering_init)#, n_trees=100, max_iterations=100) +clustering_new = iterative_classification( + data.to_numpy(), clustering_init # , n_trees=100, max_iterations=100 +) # save clustering as CSV -pd.DataFrame({"clustification_clustering": clustering_new}, index=data.index).to_csv(result_path, index=True) +pd.DataFrame({"clustification_clustering": clustering_new}, index=data.index).to_csv( + result_path, index=True +) diff --git a/workflow/scripts/clustree.R b/workflow/scripts/clustree.R index 484dfa4..0ea84ca 100644 --- a/workflow/scripts/clustree.R +++ b/workflow/scripts/clustree.R @@ -14,133 +14,142 @@ plot_path <- snakemake@output[["plot"]] # parameters content <- snakemake@wildcards[["content"]] # wildcard, one of: metadata, features, default, custom -categorical_label_option <- as.character(snakemake@params[["categorical_label_option"]]) #"pure" or "majority" +categorical_label_option <- as.character(snakemake@params[["categorical_label_option"]]) # "pure" or "majority" numerical_aggregation_option <- as.character(snakemake@params[["numerical_aggregation_option"]]) # mean, median, max, min custom_metadata <- c(snakemake@params[["custom_metadata"]]) # c("target") count_filter <- as.numeric(snakemake@params[["count_filter"]]) # 0 prop_filter <- as.numeric(snakemake@params[["prop_filter"]]) # 0.1 -layout <- as.character(snakemake@params[["layout"]]) # "tree" or "sugiyama" +layout <- as.character(snakemake@params[["layout"]]) # "tree" or "sugiyama" # helper function for labeling catgorical metadata categorical_labeler <- function(labels) { - - if (categorical_label_option == "majority"){ - label <- as.character(names(which.max(table(labels)))) - } else if (categorical_label_option == "pure"){ - if (length(unique(labels)) == 1) { - label <- as.character(unique(labels)) - } else { - label <- "mixed" - } + if (categorical_label_option == "majority") { + label <- as.character(names(which.max(table(labels)))) + } else if (categorical_label_option == "pure") { + if (length(unique(labels)) == 1) { + label <- as.character(unique(labels)) + } else { + label <- "mixed" } - return(label) + } + return(label) } # helper function for aggregation of numerical metadata -numerical_aggregation <- function(values){ - if(numerical_aggregation_option=="mean"){ - label <- mean(values) - } else if (numerical_aggregation_option=="median"){ - label <- median(values) - } else if (numerical_aggregation_option=="min"){ - label <- min(values) - } else if (numerical_aggregation_option=="max"){ - label <- max(values) - } - - return(round(label,1)) +numerical_aggregation <- function(values) { + if (numerical_aggregation_option == "mean") { + label <- mean(values) + } else if (numerical_aggregation_option == "median") { + label <- median(values) + } else if (numerical_aggregation_option == "min") { + label <- min(values) + } else if (numerical_aggregation_option == "max") { + label <- max(values) + } + + return(round(label, 1)) } plot_clustree <- function(data, col) { - clustree_plot <- clustree( - x = data, - prefix = "X_", - suffix = NULL, - #metadata = features, # not used in case of dataframe as input - count_filter = count_filter, - prop_filter = prop_filter, - layout = layout, - use_core_edges = TRUE, - highlight_core = FALSE, # check effect -> TODO - node_colour = col, # depending on col - node_colour_aggr = if(col=="X_") NULL else if(is.numeric(data[[col]])) "numerical_aggregation" else "categorical_labeler", # depending on col - node_size = "size", - node_size_aggr = NULL, - node_size_range = c(4, 15), - node_alpha = 1, - node_alpha_aggr = NULL, - node_text_size = 3, - scale_node_text = FALSE, - node_text_colour = "black", - node_label = if(col=="X_") NULL else col, # depending on col - node_label_aggr = if(col=="X_") NULL else if(is.numeric(data[[col]])) "numerical_aggregation" else "categorical_labeler", # depending on col - node_label_size = 3, - node_label_nudge = -0.2, - edge_width = 1.5, - edge_arrow = FALSE, # default: TRUE, but makes plot unnecessary busy - edge_arrow_ends = "last", - show_axis = FALSE, - return = "plot" - ) + - {if(col=="X_") guides(color="none")} + # remove legend in case of color being clusterings - {if(is.factor(data[[col]]) & length(unique(data[[col]]))>25) guides(color="none")} + # remove legend in case of color being categorical with more than 25 categories - scale_edge_color_continuous(low = "grey80", high = "firebrick2") + - {if(is.numeric(data[[col]])) scale_color_gradient2(midpoint=0, low="royalblue4", mid="grey80", high="firebrick2", space ="Lab")} + - {if(is.numeric(data[[col]])) scale_fill_gradient2(midpoint=0, low="royalblue4", mid="grey80", high="firebrick2", space ="Lab")} + - labs(color=col) + guides(fill="none")+ - geom_label(data = clustering_names, aes(x = Inf, y = index-1+0.25, label = clustering), hjust = 1, vjust = 0.5, size = 3, inherit.aes = FALSE) + # add clustering names - {if(col=="X_") ggtitle(content) else ggtitle(col)} + theme(plot.title = element_text(size = 10)) - - return(clustree_plot) + clustree_plot <- clustree( + x = data, + prefix = "X_", + suffix = NULL, + # metadata = features, # not used in case of dataframe as input + count_filter = count_filter, + prop_filter = prop_filter, + layout = layout, + use_core_edges = TRUE, + highlight_core = FALSE, # check effect -> TODO + node_colour = col, # depending on col + node_colour_aggr = if (col == "X_") NULL else if (is.numeric(data[[col]])) "numerical_aggregation" else "categorical_labeler", # depending on col + node_size = "size", + node_size_aggr = NULL, + node_size_range = c(4, 15), + node_alpha = 1, + node_alpha_aggr = NULL, + node_text_size = 3, + scale_node_text = FALSE, + node_text_colour = "black", + node_label = if (col == "X_") NULL else col, # depending on col + node_label_aggr = if (col == "X_") NULL else if (is.numeric(data[[col]])) "numerical_aggregation" else "categorical_labeler", # depending on col + node_label_size = 3, + node_label_nudge = -0.2, + edge_width = 1.5, + edge_arrow = FALSE, # default: TRUE, but makes plot unnecessary busy + edge_arrow_ends = "last", + show_axis = FALSE, + return = "plot" + ) + + { + if (col == "X_") guides(color = "none") + } + # remove legend in case of color being clusterings + { + if (is.factor(data[[col]]) & length(unique(data[[col]])) > 25) guides(color = "none") + } + # remove legend in case of color being categorical with more than 25 categories + scale_edge_color_continuous(low = "grey80", high = "firebrick2") + + { + if (is.numeric(data[[col]])) scale_color_gradient2(midpoint = 0, low = "royalblue4", mid = "grey80", high = "firebrick2", space = "Lab") + } + + { + if (is.numeric(data[[col]])) scale_fill_gradient2(midpoint = 0, low = "royalblue4", mid = "grey80", high = "firebrick2", space = "Lab") + } + + labs(color = col) + guides(fill = "none") + + geom_label(data = clustering_names, aes(x = Inf, y = index - 1 + 0.25, label = clustering), hjust = 1, vjust = 0.5, size = 3, inherit.aes = FALSE) + # add clustering names + { + if (col == "X_") ggtitle(content) else ggtitle(col) + } + theme(plot.title = element_text(size = 10)) + + return(clustree_plot) } ### load data # clusterings <- read.csv(file=file.path(clustering_path), row.names=1, header=TRUE) # metadata <- read.csv(file=file.path(metadata_path), row.names=1, header=TRUE) -clusterings <- data.frame(fread(file.path(clustering_path), header=TRUE), row.names=1) -metadata <- data.frame(fread(file.path(metadata_path), header=TRUE), row.names=1) +clusterings <- data.frame(fread(file.path(clustering_path), header = TRUE), row.names = 1) +metadata <- data.frame(fread(file.path(metadata_path), header = TRUE), row.names = 1) ### transform data # sort (rows by clustering rownames and) columns by alphabet -> rows should be sorted correctly in all cases -metadata <- metadata[,sort(colnames(metadata)),drop=FALSE] # rownames(clusterings) +metadata <- metadata[, sort(colnames(metadata)), drop = FALSE] # rownames(clusterings) # transform metadata -if (content=="metadata"){ - na_cols <- c() - for (col in colnames(metadata)){ - # if NA -> remove column and move on - if (any(is.na(metadata[[col]]))){ - na_cols <- c(na_cols, col) - next - } - # if a metadata class is empty ("") fill with "unknown" - if (!any(is.na(metadata[[col]]))){ - if (any(metadata[[col]]=="")){ - metadata[metadata[[col]]=="", col] <- "unknown" - } - } - # convert metadata to categorical if less than 25 unique integer values - if (is.numeric(metadata[[col]]) & length(unique(metadata[[col]]))<=25){ - if(all(metadata[[col]] == round(metadata[[col]]))){ - metadata[col] <- as.factor(metadata[[col]]) - } - } +if (content == "metadata") { + na_cols <- c() + for (col in colnames(metadata)) { + # if NA -> remove column and move on + if (any(is.na(metadata[[col]]))) { + na_cols <- c(na_cols, col) + next + } + # if a metadata class is empty ("") fill with "unknown" + if (!any(is.na(metadata[[col]]))) { + if (any(metadata[[col]] == "")) { + metadata[metadata[[col]] == "", col] <- "unknown" + } } - # remove columns with NA - metadata <- metadata[, !(colnames(metadata) %in% na_cols),drop=FALSE] -} else if (content=="custom"){# for custom plot with selected categorical metadata as clustering - # subset to selected metadata & add prefix to colnames "metadata_" - metadata <- metadata[,custom_metadata, drop=FALSE] - colnames(metadata) <- paste0("metadata_", colnames(metadata)) - # convert to factor - for (col in colnames(metadata)){ + # convert metadata to categorical if less than 25 unique integer values + if (is.numeric(metadata[[col]]) & length(unique(metadata[[col]])) <= 25) { + if (all(metadata[[col]] == round(metadata[[col]]))) { metadata[col] <- as.factor(metadata[[col]]) + } } - - # add to clusterings - clusterings <- cbind(clusterings, metadata) + } + # remove columns with NA + metadata <- metadata[, !(colnames(metadata) %in% na_cols), drop = FALSE] +} else if (content == "custom") { # for custom plot with selected categorical metadata as clustering + # subset to selected metadata & add prefix to colnames "metadata_" + metadata <- metadata[, custom_metadata, drop = FALSE] + colnames(metadata) <- paste0("metadata_", colnames(metadata)) + # convert to factor + for (col in colnames(metadata)) { + metadata[col] <- as.factor(metadata[[col]]) + } + + # add to clusterings + clusterings <- cbind(clusterings, metadata) } # sort by number of clusters @@ -155,46 +164,45 @@ clustering_names <- data.frame(clustering = colnames(clusterings), index = ncol( # rename columns to X_{n} colnames(clusterings) <- paste0("X_", 1:ncol(clusterings)) - + # plot specifications -width <- min(50,ceiling(length(unique(clusterings[,ncol(clusterings)]))/3) + 2) # per clustering 1/3 inch -height <- min(50,ncol(clusterings) + 2) - +width <- min(50, ceiling(length(unique(clusterings[, ncol(clusterings)])) / 3) + 2) # per clustering 1/3 inch +height <- min(50, ncol(clusterings) + 2) + # add metadata and features to clusterings -data <- cbind(clusterings, metadata)#, features) +data <- cbind(clusterings, metadata) # , features) + - ### clustree analysis -if (content=="default" | content=="custom"){ - # make default or custom plot without any metadata highlighted - clustree_final <- plot_clustree(data, "X_") - - ggsave(basename(plot_path), - plot = clustree_final, - device = 'png', - path = file.path(dirname(plot_path)), - scale = 1, - dpi = 300, - width = width, - height = height, - limitsize = FALSE, - ) -} else{ - dir.create(plot_path, recursive = TRUE) - # generate & save metadata/features plots - for (col in colnames(metadata)){ - - clustree_final <- plot_clustree(data, col) - - ggsave(paste0(col,".png"), - plot = clustree_final, - device = 'png', - path = file.path(plot_path), - scale = 1, - dpi = 300, - width = width, - height = height, - limitsize = FALSE, - ) - } +if (content == "default" | content == "custom") { + # make default or custom plot without any metadata highlighted + clustree_final <- plot_clustree(data, "X_") + + ggsave(basename(plot_path), + plot = clustree_final, + device = "png", + path = file.path(dirname(plot_path)), + scale = 1, + dpi = 300, + width = width, + height = height, + limitsize = FALSE, + ) +} else { + dir.create(plot_path, recursive = TRUE) + # generate & save metadata/features plots + for (col in colnames(metadata)) { + clustree_final <- plot_clustree(data, col) + + ggsave(paste0(col, ".png"), + plot = clustree_final, + device = "png", + path = file.path(plot_path), + scale = 1, + dpi = 300, + width = width, + height = height, + limitsize = FALSE, + ) + } } diff --git a/workflow/scripts/distance_matrix.py b/workflow/scripts/distance_matrix.py index 5188666..063ba9d 100644 --- a/workflow/scripts/distance_matrix.py +++ b/workflow/scripts/distance_matrix.py @@ -4,6 +4,7 @@ import pandas as pd from scipy.spatial.distance import pdist, squareform import math + # from fastdist import fastdist #### configurations @@ -15,8 +16,8 @@ distance_matrix_path = snakemake.output["distance_matrix"] # parameters -samples_by_features = int(snakemake.params['samples_by_features']) -data_or_feature = snakemake.wildcards["type"] # "observations" or "features" +samples_by_features = int(snakemake.params["samples_by_features"]) +data_or_feature = snakemake.wildcards["type"] # "observations" or "features" metric = snakemake.wildcards["metric"] n_observations = snakemake.config["heatmap"]["n_observations"] n_features = snakemake.config["heatmap"]["n_features"] @@ -25,25 +26,27 @@ # check data orientation to fit: samples/observations x features data = pd.read_csv(data_path, index_col=0) -if (samples_by_features == 0 and data_or_feature == "observations") or (samples_by_features == 1 and data_or_feature == "features"): +if (samples_by_features == 0 and data_or_feature == "observations") or ( + samples_by_features == 1 and data_or_feature == "features" +): data = data.T # retain only highly variable features if data_or_feature == "features": variances = data.var(axis=1) - if isinstance(n_features, float) or n_features==1: + if isinstance(n_features, float) or n_features == 1: n_features = int(math.floor(n_features * data.shape[0])) top_features = variances.nlargest(n_features).index - data = data.loc[top_features,:] + data = data.loc[top_features, :] # downsample observations if data_or_feature == "observations": - if isinstance(n_observations, float) or n_observations==1: + if isinstance(n_observations, float) or n_observations == 1: n_observations = int(math.floor(n_observations * data.shape[0])) if n_observations < data.shape[0]: data = data.sample(n=n_observations, random_state=42) - + # Convert DataFrame to NumPy array # data_np = data.to_numpy() @@ -55,7 +58,7 @@ # fastdist # metric_function = getattr(fastdist, metric) -# dist_mtx = fastdist.matrix_pairwise_distance(data_np, metric_function, metric, return_matrix=True) +# dist_mtx = fastdist.matrix_pairwise_distance(data_np, metric_function, metric, return_matrix=True) # convert to squareform dataframe dist_mtx_df = pd.DataFrame(squareform(dist_mtx), index=data.index, columns=data.index) diff --git a/workflow/scripts/leiden_cluster.py b/workflow/scripts/leiden_cluster.py index 6764b2c..a9fea92 100644 --- a/workflow/scripts/leiden_cluster.py +++ b/workflow/scripts/leiden_cluster.py @@ -4,6 +4,7 @@ import os import pickle import numpy as np + # Leiden algorithm import umap import leidenalg as la @@ -11,6 +12,7 @@ from scipy.sparse import csr_matrix import pandas as pd + # helper function adapted from here: https://dynamo-release.readthedocs.io/en/latest/_modules/dynamo/tools/connectivity.html def knn_to_adj(knn_indices: np.ndarray, knn_weights: np.ndarray) -> csr_matrix: """Convert a knn graph's indices and weights to an adjacency matrix of the corresponding nearest neighbor graph. @@ -38,6 +40,7 @@ def knn_to_adj(knn_indices: np.ndarray, knn_weights: np.ndarray) -> csr_matrix: return adj + #### configurations # inputs @@ -48,17 +51,25 @@ def knn_to_adj(knn_indices: np.ndarray, knn_weights: np.ndarray) -> csr_matrix: result_path = os.path.join(snakemake.output["clustering"]) # UMAP parameters for small data (<11 observations) -samples_by_features = int(snakemake.params['samples_by_features']) #0 -metric = snakemake.params['metric'] #"correlation" -n_neighbors = int(snakemake.params['n_neighbors']) #100 +samples_by_features = int(snakemake.params["samples_by_features"]) # 0 +metric = snakemake.params["metric"] # "correlation" +n_neighbors = int(snakemake.params["n_neighbors"]) # 100 # Leiden algorithm parameters -n_iterations = int(snakemake.params["n_iterations"]) # 2 +n_iterations = int(snakemake.params["n_iterations"]) # 2 # Get the partition method from the leidenalg module -partition_type = getattr(la, str(snakemake.params["partition_type"])) #"RBConfigurationVertexPartition" +partition_type = getattr( + la, str(snakemake.params["partition_type"]) +) # "RBConfigurationVertexPartition" # Get the kwargs from the config dict -if str(snakemake.params["partition_type"]) in ["RBConfigurationVertexPartition", "RBERVertexPartition", "CPMVertexPartition"]: - la_kwargs = {"resolution_parameter": float(snakemake.params["resolution"])} #0.05 #1.0 +if str(snakemake.params["partition_type"]) in [ + "RBConfigurationVertexPartition", + "RBERVertexPartition", + "CPMVertexPartition", +]: + la_kwargs = { + "resolution_parameter": float(snakemake.params["resolution"]) # 0.05 #1.0 + } else: la_kwargs = {} @@ -75,53 +86,54 @@ def knn_to_adj(knn_indices: np.ndarray, knn_weights: np.ndarray) -> csr_matrix: data = pd.read_csv(data_path, index_col=0).T # if less than 11 datapoints there is no pre-computed KNN graph -if data.shape[0]<11: +if data.shape[0] < 11: # run UMAP - umap_obj = umap.umap_.UMAP(n_neighbors=n_neighbors, - n_components=2, - metric=metric, - metric_kwds=None, - output_metric='euclidean', - output_metric_kwds=None, - n_epochs=None, - learning_rate=1.0, - init='spectral', - min_dist=0.1, - spread=1.0, - low_memory=True, - n_jobs=-1, - set_op_mix_ratio=1.0, - local_connectivity=1.0, - repulsion_strength=1.0, - negative_sample_rate=5, - transform_queue_size=4.0, - a=None, - b=None, - random_state=42, - angular_rp_forest=False, - target_n_neighbors=-1, - target_metric='categorical', - target_metric_kwds=None, - target_weight=0.5, - transform_seed=42, - transform_mode='embedding', - force_approximation_algorithm=False, - verbose=False, - tqdm_kwds=None, - unique=False, - densmap=False, - dens_lambda=2.0, - dens_frac=0.3, - dens_var_shift=0.1, - output_dens=False, - disconnection_distance=None, - precomputed_knn=(None, None, None) - ).fit(data) + umap_obj = umap.umap_.UMAP( + n_neighbors=n_neighbors, + n_components=2, + metric=metric, + metric_kwds=None, + output_metric="euclidean", + output_metric_kwds=None, + n_epochs=None, + learning_rate=1.0, + init="spectral", + min_dist=0.1, + spread=1.0, + low_memory=True, + n_jobs=-1, + set_op_mix_ratio=1.0, + local_connectivity=1.0, + repulsion_strength=1.0, + negative_sample_rate=5, + transform_queue_size=4.0, + a=None, + b=None, + random_state=42, + angular_rp_forest=False, + target_n_neighbors=-1, + target_metric="categorical", + target_metric_kwds=None, + target_weight=0.5, + transform_seed=42, + transform_mode="embedding", + force_approximation_algorithm=False, + verbose=False, + tqdm_kwds=None, + unique=False, + densmap=False, + dens_lambda=2.0, + dens_frac=0.3, + dens_var_shift=0.1, + output_dens=False, + disconnection_distance=None, + precomputed_knn=(None, None, None), + ).fit(data) # extract graph adj_coo = umap_obj.graph_.tocoo() else: # load pre-computed KNN graph - with open(graph_path, 'rb') as f: + with open(graph_path, "rb") as f: # The protocol version used is detected automatically, so we do not # have to specify it. knn = pickle.load(f) @@ -133,20 +145,23 @@ def knn_to_adj(knn_indices: np.ndarray, knn_weights: np.ndarray) -> csr_matrix: # Convert the adjacency matrix in COO format to igraph object (weighted undirected graph) edges = np.column_stack((adj_coo.row, adj_coo.col)) graph = Graph(edges.tolist(), directed=False) -graph.es['weight'] = adj_coo.data.tolist() +graph.es["weight"] = adj_coo.data.tolist() # Perform Leiden clustering -partition = la.find_partition(graph=graph, - partition_type=partition_type, - initial_membership=None, # default - weights='weight', # default: None - n_iterations=n_iterations, - max_comm_size=0, - seed=42, - **la_kwargs - ) +partition = la.find_partition( + graph=graph, + partition_type=partition_type, + initial_membership=None, # default + weights="weight", # default: None + n_iterations=n_iterations, + max_comm_size=0, + seed=42, + **la_kwargs +) # save clustering as CSV clustering_name = os.path.splitext(os.path.basename(result_path))[0] -pd.DataFrame({clustering_name: partition.membership}, index=data.index).to_csv(result_path, index=True) +pd.DataFrame({clustering_name: partition.membership}, index=data.index).to_csv( + result_path, index=True +) diff --git a/workflow/scripts/mcdm_topsis.py b/workflow/scripts/mcdm_topsis.py index 30146af..4e9b3cb 100644 --- a/workflow/scripts/mcdm_topsis.py +++ b/workflow/scripts/mcdm_topsis.py @@ -5,7 +5,6 @@ import numpy as np import pandas as pd from pymcdm.methods import TOPSIS -from pymcdm.helpers import rrankdata #### configurations @@ -28,7 +27,7 @@ # perform TOPSIS # give all scores equal weights as they are supposed to be complementary -weights = np.full(indices.shape[1], 1.0/indices.shape[1]) +weights = np.full(indices.shape[1], 1.0 / indices.shape[1]) # set 1 for benefit (max) and -1 for cost (min) functions (hard coded for internal indices in order) types = np.array([1, 1, 1, -1, -1, -1]) # create object with defaults diff --git a/workflow/scripts/pca.py b/workflow/scripts/pca.py index 33cfb9b..4ca8429 100644 --- a/workflow/scripts/pca.py +++ b/workflow/scripts/pca.py @@ -5,6 +5,7 @@ import os import pickle import pandas as pd + # dimensionality reduction from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA @@ -24,7 +25,7 @@ result_axes_path = snakemake.output["result_axes"] # parameters -samples_by_features = int(snakemake.params['samples_by_features']) +samples_by_features = int(snakemake.params["samples_by_features"]) n_components = snakemake.config["pca"]["n_components"] svd_solver = snakemake.config["pca"]["svd_solver"] @@ -35,45 +36,54 @@ data = pd.read_csv(data_path, index_col=0) else: data = pd.read_csv(data_path, index_col=0).T - + ### transform data - + # unsupervised PCA -pca_obj = PCA(n_components=n_components, - copy=True, - whiten=False, - svd_solver=svd_solver, - tol=0.0, - iterated_power='auto', - random_state=42 - ) +pca_obj = PCA( + n_components=n_components, + copy=True, + whiten=False, + svd_solver=svd_solver, + tol=0.0, + iterated_power="auto", + random_state=42, +) data_pca = pca_obj.fit_transform(StandardScaler().fit_transform(data)) -data_df = pd.DataFrame(data_pca, index=data.index,) +data_df = pd.DataFrame( + data_pca, + index=data.index, +) data_df = data_df.rename_axis(("sample_name")) -data_df.columns = ["PC_{}".format(str(idx+1)) for idx in data_df.columns] +data_df.columns = ["PC_{}".format(str(idx + 1)) for idx in data_df.columns] ### save data # save pca object -with open(result_object_path, 'wb') as f: +with open(result_object_path, "wb") as f: pickle.dump(pca_obj, f, pickle.HIGHEST_PROTOCOL) - + # save transformed data data_df.to_csv(result_data_path) -data_df.iloc[:,:min(10,data_df.shape[1])].to_csv(result_data_small_path) +data_df.iloc[:, : min(10, data_df.shape[1])].to_csv(result_data_small_path) # save loadings -loadings = pd.DataFrame(pca_obj.components_.T, columns = data_df.columns, index=data.columns) +loadings = pd.DataFrame( + pca_obj.components_.T, columns=data_df.columns, index=data.columns +) loadings.to_csv(result_loadings_path) -loadings.iloc[:,:min(10,data_df.shape[1])].to_csv(result_loadings_small_path) +loadings.iloc[:, : min(10, data_df.shape[1])].to_csv(result_loadings_small_path) # save explained variance axes_info_df = pd.DataFrame(pca_obj.explained_variance_ratio_) axes_info_df.to_csv(result_var_path) # save axes information for visualization -axes_info_df.columns = ['label'] -axes_info_df['label'] = ["PC{} ({}%)".format(idx+1, round(100*var,2)) for idx, var in axes_info_df['label'].items()] +axes_info_df.columns = ["label"] +axes_info_df["label"] = [ + "PC{} ({}%)".format(idx + 1, round(100 * var, 2)) + for idx, var in axes_info_df["label"].items() +] axes_info_df.to_csv(result_axes_path) diff --git a/workflow/scripts/plot_2d.R b/workflow/scripts/plot_2d.R index 2fe69aa..739fdfb 100644 --- a/workflow/scripts/plot_2d.R +++ b/workflow/scripts/plot_2d.R @@ -5,11 +5,13 @@ library("data.table") # utility function to adapt legend according to metadata addSmallLegend <- function(myPlot, pointSize = 2, textSize = 3, spaceLegend = 0, alpha = 1) { - new_plot <- myPlot + - guides(color = guide_legend(override.aes = list(size = pointSize, alpha = alpha))) + - theme(legend.text = element_text(size = textSize), - legend.key.size = unit(spaceLegend, "lines")) - return (new_plot) + new_plot <- myPlot + + guides(color = guide_legend(override.aes = list(size = pointSize, alpha = alpha))) + + theme( + legend.text = element_text(size = textSize), + legend.key.size = unit(spaceLegend, "lines") + ) + return(new_plot) } ### configurations @@ -23,105 +25,109 @@ metadata_path <- snakemake@input[["metadata"]] plot_path <- snakemake@output[["plot"]] # parameters -size <- snakemake@params[["size"]]# 0.5 -alpha <- snakemake@params[["alpha"]]# 1 -coord_fixed_flag <- if(as.integer(snakemake@config[["coord_fixed"]])==1) TRUE else FALSE +size <- snakemake@params[["size"]] # 0.5 +alpha <- snakemake@params[["alpha"]] # 1 +coord_fixed_flag <- if (as.integer(snakemake@config[["coord_fixed"]]) == 1) TRUE else FALSE dir.create(plot_path, recursive = TRUE) ### load data -data <- data.frame(fread(file.path(data_path), header=TRUE), row.names=1)[,1:2] -axes <- data.frame(fread(file.path(axes_path), header=TRUE), row.names=1)[1:2,] -metadata <- data.frame(fread(file.path(metadata_path), header=TRUE), row.names=1) +data <- data.frame(fread(file.path(data_path), header = TRUE), row.names = 1)[, 1:2] +axes <- data.frame(fread(file.path(axes_path), header = TRUE), row.names = 1)[1:2, ] +metadata <- data.frame(fread(file.path(metadata_path), header = TRUE), row.names = 1) # plot specifications width <- 4 height <- 3 -shape <- if(nrow(data)>50000) '.' else 16 +shape <- if (nrow(data) > 50000) "." else 16 ### make plots -for (col in sort(colnames(metadata))){ - print(col) - - # check if metadata column is only NA - if(all(is.na(metadata[[col]]))){ - next - } - - # convert to categorical if the substring "cluster" is present - if (grepl("cluster", col)) { - metadata[col] <- as.factor(metadata[[col]]) +for (col in sort(colnames(metadata))) { + print(col) + + # check if metadata column is only NA + if (all(is.na(metadata[[col]]))) { + next + } + + # convert to categorical if the substring "cluster" is present + if (grepl("cluster", col)) { + metadata[col] <- as.factor(metadata[[col]]) + } + + # convert metadata to categorical if less than 25 unique integer values + if (grepl("metadata", basename(plot_path)) & is.numeric(metadata[[col]]) & length(unique(metadata[[col]])) <= 25) { + if (all(metadata[[col]] == round(metadata[[col]]))) { + metadata[col] <- as.factor(metadata[[col]]) } - - # convert metadata to categorical if less than 25 unique integer values - if (grepl("metadata", basename(plot_path)) & is.numeric(metadata[[col]]) & length(unique(metadata[[col]]))<=25){ - if(all(metadata[[col]] == round(metadata[[col]]))){ - metadata[col] <- as.factor(metadata[[col]]) - } + } + + # if a metadata class is empty ("") fill with "unknown" + if (!any(is.na(metadata[[col]]))) { + if (any(metadata[[col]] == "")) { + metadata[metadata[[col]] == "", col] <- "unknown" } - - # if a metadata class is empty ("") fill with "unknown" - if (!any(is.na(metadata[[col]]))){ - if (any(metadata[[col]]=="")){ - metadata[metadata[[col]]=="", col] <- "unknown" - } + } + + # prepare data for plotting + tmp_data <- cbind(data, metadata[col]) + + # make 2D scatter plots + if (!is.numeric(metadata[[col]])) { + # plot categorical data + tmp_plot <- ggplot(tmp_data, aes_string(x = colnames(tmp_data)[1], y = colnames(tmp_data)[2])) + + geom_point(aes_string(color = col), size = size, stroke = 0, alpha = alpha, shape = shape) + + { + if (coord_fixed_flag) coord_fixed() + } + + xlab(axes[1]) + + ylab(axes[2]) + + ggtitle(col) + + theme_linedraw() + + theme(plot.title = element_text(size = 10), legend.title = element_blank()) + + if (length(unique(metadata[[col]])) > 25) { + tmp_plot <- tmp_plot + theme(legend.position = "none") # addSmallLegend(myPlot = tmp_plot, pointSize = 0, textSize = 0) + } else if (length(unique(metadata[[col]])) > 15) { + tmp_plot <- addSmallLegend(myPlot = tmp_plot, textSize = 4) + } else { + tmp_plot <- addSmallLegend(myPlot = tmp_plot, textSize = 10) } - - # prepare data for plotting - tmp_data <- cbind(data, metadata[col]) - - # make 2D scatter plots - if (!is.numeric(metadata[[col]])){ - # plot categorical data - tmp_plot <- ggplot(tmp_data, aes_string(x=colnames(tmp_data)[1], y=colnames(tmp_data)[2])) + - geom_point(aes_string(color=col), size=size, stroke=0, alpha=alpha, shape=shape) + - {if(coord_fixed_flag) coord_fixed()} + - xlab(axes[1]) + - ylab(axes[2]) + - ggtitle(col) + - theme_linedraw() + - theme(plot.title = element_text(size = 10), legend.title = element_blank()) - - if (length(unique(metadata[[col]]))>25){ - tmp_plot <- tmp_plot + theme(legend.position="none") #addSmallLegend(myPlot = tmp_plot, pointSize = 0, textSize = 0) - }else if (length(unique(metadata[[col]]))>15){ - tmp_plot <- addSmallLegend(myPlot = tmp_plot, textSize = 4) - }else{ - tmp_plot <- addSmallLegend(myPlot = tmp_plot, textSize = 10) - } - }else{ - # plot numerical data - tmp_plot <- ggplot(tmp_data, aes_string(x=colnames(tmp_data)[1], y=colnames(tmp_data)[2])) + - geom_point(aes_string(color=col), size=size, stroke=0, alpha=alpha, shape=shape) + - {if(coord_fixed_flag) coord_fixed()} + - xlab(axes[1]) + - ylab(axes[2]) + - ggtitle(col) + - theme_linedraw() + - theme(plot.title = element_text(size = 10), legend.title = element_blank()) - - # if the data is ranging from negative to positive values, use a diverging color scale - if (min(tmp_data[[col]], na.rm = TRUE) < 0 & max(tmp_data[[col]], na.rm = TRUE) > 0){ - tmp_plot <- tmp_plot + - scale_color_gradient2(midpoint=0, low="royalblue4", mid="grey80", high="firebrick2", space ="Lab") - }else{ - # otherwise use a sequential color scale (mostly because then fixing the midpoint to 0 shifts the scale so - # that if you have a range from e.g., 12-14, the values only range in the high red part of the scale) - tmp_plot <- tmp_plot + - scale_color_viridis_c(option = "magma", direction = -1) - } + } else { + # plot numerical data + tmp_plot <- ggplot(tmp_data, aes_string(x = colnames(tmp_data)[1], y = colnames(tmp_data)[2])) + + geom_point(aes_string(color = col), size = size, stroke = 0, alpha = alpha, shape = shape) + + { + if (coord_fixed_flag) coord_fixed() + } + + xlab(axes[1]) + + ylab(axes[2]) + + ggtitle(col) + + theme_linedraw() + + theme(plot.title = element_text(size = 10), legend.title = element_blank()) + + # if the data is ranging from negative to positive values, use a diverging color scale + if (min(tmp_data[[col]], na.rm = TRUE) < 0 & max(tmp_data[[col]], na.rm = TRUE) > 0) { + tmp_plot <- tmp_plot + + scale_color_gradient2(midpoint = 0, low = "royalblue4", mid = "grey80", high = "firebrick2", space = "Lab") + } else { + # otherwise use a sequential color scale (mostly because then fixing the midpoint to 0 shifts the scale so + # that if you have a range from e.g., 12-14, the values only range in the high red part of the scale) + tmp_plot <- tmp_plot + + scale_color_viridis_c(option = "magma", direction = -1) } - - # save plot - ggsave(paste0(col,".png"), - plot = tmp_plot, - device = 'png', - path = plot_path, - scale = 1, - dpi = 300, - width = width, - height = height, - limitsize = FALSE, - ) + } + + # save plot + ggsave(paste0(col, ".png"), + plot = tmp_plot, + device = "png", + path = plot_path, + scale = 1, + dpi = 300, + width = width, + height = height, + limitsize = FALSE, + ) } diff --git a/workflow/scripts/plot_heatmap.R b/workflow/scripts/plot_heatmap.R index 8b05261..4f986a4 100644 --- a/workflow/scripts/plot_heatmap.R +++ b/workflow/scripts/plot_heatmap.R @@ -4,7 +4,7 @@ library("circlize") library("RColorBrewer") library("data.table") library("fastcluster") -#library("dendsort") +# library("dendsort") ### configurations set.seed(42) @@ -20,44 +20,44 @@ features_distance_path <- snakemake@input[["features_distance"]] plot_path <- snakemake@output[["plot"]] # parameters -samples_by_features <- as.integer(snakemake@params['samples_by_features']) +samples_by_features <- as.integer(snakemake@params["samples_by_features"]) metric <- snakemake@wildcards[["metric"]] cluster_method <- snakemake@wildcards[["method"]] metadata_col <- c(snakemake@config[["metadata_of_interest"]])[1] ### load data -data <- data.frame(fread(file.path(data_path), header=TRUE), row.names=1) -metadata <- data.frame(fread(file.path(metadata_path), header=TRUE), row.names=1) -observations_distance <- fread(file.path(observations_distance_path), header=TRUE) -features_distance <- fread(file.path(features_distance_path), header=TRUE) +data <- data.frame(fread(file.path(data_path), header = TRUE), row.names = 1) +metadata <- data.frame(fread(file.path(metadata_path), header = TRUE), row.names = 1) +observations_distance <- fread(file.path(observations_distance_path), header = TRUE) +features_distance <- fread(file.path(features_distance_path), header = TRUE) ### prepare distance matrices # Remove rownames and convert to matrix -observations_distance <- as.matrix(observations_distance[ , -1, with = FALSE]) -features_distance <- as.matrix(features_distance[ , -1, with = FALSE]) +observations_distance <- as.matrix(observations_distance[, -1, with = FALSE]) +features_distance <- as.matrix(features_distance[, -1, with = FALSE]) # add rownames rownames(observations_distance) <- colnames(observations_distance) rownames(features_distance) <- colnames(features_distance) # check and fix orientation -if(samples_by_features==0){ - data <- t(data) +if (samples_by_features == 0) { + data <- t(data) } # filter data and metadata in case of downsampled observations and features data <- data[colnames(observations_distance), colnames(features_distance)] -metadata <- metadata[colnames(observations_distance),, drop=FALSE] +metadata <- metadata[colnames(observations_distance), , drop = FALSE] # scale data data <- scale(data) # remove complete NA columns & replace remaining NA values -data <- data[,colSums(is.na(data))100, FALSE, TRUE) -show_column_names <- ifelse(ncol(data)>100, FALSE, TRUE) +show_row_names <- ifelse(nrow(data) > 100, FALSE, TRUE) +show_column_names <- ifelse(ncol(data) > 100, FALSE, TRUE) # # alternative way to: determine distance, hierarchical clustering and order dendrograms # # would be used in the parameters cluster_rows=row_dend and cluster_columns=col_dend as arguments @@ -159,7 +158,7 @@ show_column_names <- ifelse(ncol(data)>100, FALSE, TRUE) # row_data[] <- sapply(row_data, function(x) 1-abs(x)) # row_data <- as.dist(row_data) # row_dend <- dendsort(hclust(row_data, method=cluster_method)) - + # #cols # col_data <- cor(data, method = metric) # col_data[] <- sapply(col_data, function(x) 1-abs(x)) @@ -171,25 +170,25 @@ show_column_names <- ifelse(ncol(data)>100, FALSE, TRUE) ### make & save heatmap # options(repr.plot.width=10, repr.plot.height=10) -png(filename=plot_path, width=20, height=20, units = "in", res=300) +png(filename = plot_path, width = 20, height = 20, units = "in", res = 300) Heatmap(data, - name = "z-scores", - column_title = paste0("Heatmap of data scaled by features (z-scores), hierarchically clustered using method ",cluster_method," with distance metric ",metric,", and colorscape limited to the top percentiles."), - col = col_fun, - left_annotation = row_annot, - show_row_names = show_row_names, - show_column_names = show_column_names, - cluster_rows = as.dendrogram(obs_hc), - cluster_columns = as.dendrogram(feat_hc), -# clustering_distance_rows = metric, -# clustering_method_rows = cluster_method, -# clustering_distance_columns = metric, -# clustering_method_columns = cluster_method, - row_dend_reorder = TRUE, - column_dend_reorder = TRUE, - use_raster = TRUE, - raster_quality = 9 - ) + name = "z-scores", + column_title = paste0("Heatmap of data scaled by features (z-scores), hierarchically clustered using method ", cluster_method, " with distance metric ", metric, ", and colorscape limited to the top percentiles."), + col = col_fun, + left_annotation = row_annot, + show_row_names = show_row_names, + show_column_names = show_column_names, + cluster_rows = as.dendrogram(obs_hc), + cluster_columns = as.dendrogram(feat_hc), + # clustering_distance_rows = metric, + # clustering_method_rows = cluster_method, + # clustering_distance_columns = metric, + # clustering_method_columns = cluster_method, + row_dend_reorder = TRUE, + column_dend_reorder = TRUE, + use_raster = TRUE, + raster_quality = 9 +) dev.off() diff --git a/workflow/scripts/plot_indices.R b/workflow/scripts/plot_indices.R index d64e3d2..18966ce 100644 --- a/workflow/scripts/plot_indices.R +++ b/workflow/scripts/plot_indices.R @@ -14,7 +14,7 @@ index_paths <- snakemake@input plot_path <- snakemake@output[["plot"]] indices <- names(index_paths) -indices <- indices[indices!=""] +indices <- indices[indices != ""] # make result directory dir.create(plot_path, recursive = TRUE) @@ -22,62 +22,69 @@ dir.create(plot_path, recursive = TRUE) # heatmaps <- list() # loop through all proivided index paths and create (clustered) heatmap -for (idx in indices){ - -# scores <- read.csv(file.path(index_paths[[idx]]), row.names = 1) - scores <- data.frame(fread(file.path(index_paths[[idx]]), header=TRUE), row.names=1) - - # scale internal indices, add seperator and reorder columns (max/sep/min) - if(length(indices)==1){ - scores <- as.data.frame(scale(scores)) - scores$sep <- NA - scores <- scores[,c("Silhouette", "Calinski_Harabasz", "Dunn", "sep", "C_index", "Davies_Bouldin", "BIC")] - } - - height <- nrow(scores) * 0.5 - width <- ncol(scores)*0.5 + 5 - - # Perform hierarchical clustering on the rows and columns of external indices - if(nrow(scores)>2 & ncol(scores)>2 & length(indices)!=1){ - cluster_rows <- hclust(dist(scores)) - cluster_cols <- hclust(dist(t(scores))) - scores <- scores[order.dendrogram(as.dendrogram(cluster_rows)), order.dendrogram(as.dendrogram(cluster_cols))] - } - - # prepare plotting dataframe - scores$clusterings <- rownames(scores) - scores_plot <- melt(scores, id.vars = "clusterings", variable.name = "metadata") - - # retain MCDM ranking - if(length(indices)==1){ - scores_plot$clusterings <- factor(scores_plot$clusterings, levels = rev(unique(scores_plot$clusterings))) - } - - # Create a heatmap - heatmap_tmp <- ggplot(scores_plot, aes(x=metadata, y=clusterings, fill=value)) + - geom_tile(color = "white", size = 1) + - geom_text(aes(label=round(value, 2)), size=3) + - theme_minimal() + - labs(title=gsub("_", " ", idx)) + - {if(length(indices)==1) xlab("internal cluster indices (benefit/max || cost/min)") } + - {if(length(indices)==1) scale_x_discrete(labels = function(x) gsub(".*\\.", "", x), breaks = function(x) {x[x!="sep"]}) } + - scale_fill_gradient2(midpoint=0, low="royalblue4", mid="white", high="firebrick2", space ="Lab") + - theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotate x-axis labels by 45 degrees - -# options(repr.plot.width=width, repr.plot.height=height) -# print(heatmap_tmp) - - # save plot - ggsave(paste0(idx,".png"), - plot = heatmap_tmp, - device = 'png', - path = plot_path, - scale = 1, - dpi = 300, - width = width, - height = height, - limitsize = FALSE, - ) +for (idx in indices) { + # scores <- read.csv(file.path(index_paths[[idx]]), row.names = 1) + scores <- data.frame(fread(file.path(index_paths[[idx]]), header = TRUE), row.names = 1) + + # scale internal indices, add seperator and reorder columns (max/sep/min) + if (length(indices) == 1) { + scores <- as.data.frame(scale(scores)) + scores$sep <- NA + scores <- scores[, c("Silhouette", "Calinski_Harabasz", "Dunn", "sep", "C_index", "Davies_Bouldin", "BIC")] + } + + height <- nrow(scores) * 0.5 + width <- ncol(scores) * 0.5 + 5 + + # Perform hierarchical clustering on the rows and columns of external indices + if (nrow(scores) > 2 & ncol(scores) > 2 & length(indices) != 1) { + cluster_rows <- hclust(dist(scores)) + cluster_cols <- hclust(dist(t(scores))) + scores <- scores[order.dendrogram(as.dendrogram(cluster_rows)), order.dendrogram(as.dendrogram(cluster_cols))] + } + + # prepare plotting dataframe + scores$clusterings <- rownames(scores) + scores_plot <- melt(scores, id.vars = "clusterings", variable.name = "metadata") + + # retain MCDM ranking + if (length(indices) == 1) { + scores_plot$clusterings <- factor(scores_plot$clusterings, levels = rev(unique(scores_plot$clusterings))) + } + + # Create a heatmap + heatmap_tmp <- ggplot(scores_plot, aes(x = metadata, y = clusterings, fill = value)) + + geom_tile(color = "white", size = 1) + + geom_text(aes(label = round(value, 2)), size = 3) + + theme_minimal() + + labs(title = gsub("_", " ", idx)) + + { + if (length(indices) == 1) xlab("internal cluster indices (benefit/max || cost/min)") + } + + { + if (length(indices) == 1) { + scale_x_discrete(labels = function(x) gsub(".*\\.", "", x), breaks = function(x) { + x[x != "sep"] + }) + } + } + + scale_fill_gradient2(midpoint = 0, low = "royalblue4", mid = "white", high = "firebrick2", space = "Lab") + + theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotate x-axis labels by 45 degrees + + # options(repr.plot.width=width, repr.plot.height=height) + # print(heatmap_tmp) + + # save plot + ggsave(paste0(idx, ".png"), + plot = heatmap_tmp, + device = "png", + path = plot_path, + scale = 1, + dpi = 300, + width = width, + height = height, + limitsize = FALSE, + ) } # differentiate between internal (one heatmap) and external (panel of six heatmaps) indices diff --git a/workflow/scripts/plot_interactive.py b/workflow/scripts/plot_interactive.py index ed32629..be4b2dc 100644 --- a/workflow/scripts/plot_interactive.py +++ b/workflow/scripts/plot_interactive.py @@ -2,49 +2,53 @@ #### libraries # general -import os import numpy as np import pandas as pd + # plotting -import seaborn as sns import matplotlib.pyplot as plt import plotly.express as px + # helper function for button to determine colors of categorical metadata def config_button_cat(fig, data, cat_var): - + unique_vals = data[cat_var].unique() unique_vals = unique_vals[pd.notna(unique_vals)] - - if len(unique_vals)<10: - cm = plt.get_cmap('tab10') - elif len(unique_vals)<20: - cm = plt.get_cmap('tab20') + + if len(unique_vals) < 10: + cm = plt.get_cmap("tab10") + elif len(unique_vals) < 20: + cm = plt.get_cmap("tab20") else: - cm = plt.get_cmap('gist_ncar') + cm = plt.get_cmap("gist_ncar") -# cm = plt.get_cmap('gist_ncar') - colors = [cm(1.*i/(len(unique_vals)+1)) for i in range(len(unique_vals)+1)] -# colors = [cm(1.*i/(len(unique_vals))) for i in range(len(unique_vals))] + # cm = plt.get_cmap('gist_ncar') + colors = [cm(1.0 * i / (len(unique_vals) + 1)) for i in range(len(unique_vals) + 1)] + # colors = [cm(1.*i/(len(unique_vals))) for i in range(len(unique_vals))] color_map = dict(zip(unique_vals, colors)) - - tmp_idx = list(data.columns).index(cat_var)-dimensions - data_colors = [color_map[datapoint[tmp_idx]] for datapoint in fig["data"][0]["customdata"]] + + tmp_idx = list(data.columns).index(cat_var) - dimensions + data_colors = [ + color_map[datapoint[tmp_idx]] for datapoint in fig["data"][0]["customdata"] + ] return [data_colors] + # helper function to put both plotly figures into the same HTML next to each other (independently) def figures_to_html(figs, filename="dashboard.html"): - with open(filename, 'w') as dashboard: + with open(filename, "w") as dashboard: dashboard.write("" + "\n") dashboard.write("
" + "\n") for fig in figs: dashboard.write("
" + "\n") - inner_html = fig.to_html().split('')[1].split('')[0] + inner_html = fig.to_html().split("")[1].split("")[0] dashboard.write(inner_html) dashboard.write("
" + "\n") dashboard.write("
" + "\n") dashboard.write("" + "\n") + #### configurations # inputs @@ -61,9 +65,11 @@ def figures_to_html(figs, filename="dashboard.html"): plot_path = snakemake.output["plot"] # parameters -dimensions = int(snakemake.params["n_components"]) #2 -point_size = 2*snakemake.params["size"] if dimensions==3 else 5*snakemake.params["size"] # 2 -point_alpha = snakemake.params["alpha"] #1 +dimensions = int(snakemake.params["n_components"]) # 2 +point_size = ( + 2 * snakemake.params["size"] if dimensions == 3 else 5 * snakemake.params["size"] +) # 2 +point_alpha = snakemake.params["alpha"] # 1 width = 750 height = 750 @@ -73,96 +79,103 @@ def figures_to_html(figs, filename="dashboard.html"): metadata = pd.read_csv(metadata_path, index_col=0) # fix metadata indices if they do not agree with data as they come from outside the workflow (e.g., R) -if not(all(data.index==metadata.index)): - #metadata.index = metadata.index.map(str) - if metadata.index.inferred_type=='string': - metadata.index = [idx.replace('-','.') for idx in metadata.index] +if not (all(data.index == metadata.index)): + # metadata.index = metadata.index.map(str) + if metadata.index.inferred_type == "string": + metadata.index = [idx.replace("-", ".") for idx in metadata.index] metadata_features = pd.read_csv(metadata_features_path, index_col=0) -data_all = pd.concat([data.iloc[:,:dimensions], metadata, metadata_features], axis=1) -data_all = data_all.fillna('') +data_all = pd.concat([data.iloc[:, :dimensions], metadata, metadata_features], axis=1) +data_all = data_all.fillna("") -# sort metadata by data type +# sort metadata by data type meta_num = list() meta_cat = list() for variable in data_all.columns[dimensions:]: unique_vals = list(data_all[variable].unique()) - #unique_vals = unique_vals[~np.isnan(unique_vals)] #unique_vals[pd.notna(unique_vals)] - + # unique_vals = unique_vals[~np.isnan(unique_vals)] #unique_vals[pd.notna(unique_vals)] + # check if integer AND less than 25 unique values -> categorical metadata - if all([isinstance(i, (int, np.int64)) for i in unique_vals]) and len(unique_vals)<25: - #data_all[variable] = data_all[variable].values.astype(str) + if ( + all([isinstance(i, (int, np.int64)) for i in unique_vals]) + and len(unique_vals) < 25 + ): + # data_all[variable] = data_all[variable].values.astype(str) meta_cat.append(variable) continue - + if all([isinstance(i, (str, bool, np.bool_)) for i in unique_vals]): -# print('discrete variable ', variable) + # print('discrete variable ', variable) meta_cat.append(variable) - + elif all([isinstance(i, (int, float, np.int64)) for i in unique_vals]): -# print('continous variable ', variable) + # print('continous variable ', variable) meta_num.append(variable) - + else: print("variable type not-detected for {}".format(variable)) # if clustering results are provided add them as categorical data -if metadata_clusterings_path!="": +if metadata_clusterings_path != "": metadata_clusterings = pd.read_csv(metadata_clusterings_path, index_col=0) data_all = pd.concat([data_all, metadata_clusterings], axis=1) meta_cat = meta_cat + metadata_clusterings.columns.tolist() - + # plotting the interactive scatter plot -#in 2D -if dimensions==2: - fig_num = px.scatter(data_all, - x=data_all.columns[0], - y=data_all.columns[1], - hover_data=meta_cat, - custom_data=list(data_all.columns)[dimensions:], - width=width, - height=height, - opacity=point_alpha, - title="Numerical Metadata", - ) - fig_cat = px.scatter(data_all, - x=data_all.columns[0], - y=data_all.columns[1], - hover_data=meta_cat, - custom_data=list(data_all.columns)[dimensions:], - width=width, - height=height, - opacity=point_alpha, - title="Categorical Metadata", - render_mode = "webgl" # required for less than 1000 datapoints, otherwise metadata selection does not work - ) +# in 2D +if dimensions == 2: + fig_num = px.scatter( + data_all, + x=data_all.columns[0], + y=data_all.columns[1], + hover_data=meta_cat, + custom_data=list(data_all.columns)[dimensions:], + width=width, + height=height, + opacity=point_alpha, + title="Numerical Metadata", + ) + fig_cat = px.scatter( + data_all, + x=data_all.columns[0], + y=data_all.columns[1], + hover_data=meta_cat, + custom_data=list(data_all.columns)[dimensions:], + width=width, + height=height, + opacity=point_alpha, + title="Categorical Metadata", + render_mode="webgl", # required for less than 1000 datapoints, otherwise metadata selection does not work + ) # in 3D -elif dimensions==3: - fig_num = px.scatter_3d(data_all, - x=data_all.columns[0], - y=data_all.columns[1], - z=data_all.columns[2], - hover_data=meta_cat, - custom_data=list(data_all.columns)[dimensions:], - width=width, - height=height, - opacity=point_alpha, - title="Numerical Metadata", - ) - fig_cat = px.scatter_3d(data_all, - x=data_all.columns[0], - y=data_all.columns[1], - z=data_all.columns[2], - hover_data=meta_cat, - custom_data=list(data_all.columns)[dimensions:], - width=width, - height=height, - opacity=point_alpha, - title="Categorical Metadata", - ) - -# set point size +elif dimensions == 3: + fig_num = px.scatter_3d( + data_all, + x=data_all.columns[0], + y=data_all.columns[1], + z=data_all.columns[2], + hover_data=meta_cat, + custom_data=list(data_all.columns)[dimensions:], + width=width, + height=height, + opacity=point_alpha, + title="Numerical Metadata", + ) + fig_cat = px.scatter_3d( + data_all, + x=data_all.columns[0], + y=data_all.columns[1], + z=data_all.columns[2], + hover_data=meta_cat, + custom_data=list(data_all.columns)[dimensions:], + width=width, + height=height, + opacity=point_alpha, + title="Categorical Metadata", + ) + +# set point size fig_num.update_traces(marker=dict(size=point_size)) fig_cat.update_traces(marker=dict(size=point_size)) @@ -182,24 +195,28 @@ def figures_to_html(figs, filename="dashboard.html"): "label": variable, "method": "update", "args": [ - {'legendgroup': '', - 'marker': {'color': data_all[variable].to_numpy(), - 'coloraxis': 'coloraxis', - 'symbol': 'circle', - 'size':point_size}, - 'mode': 'markers', - 'name': '', - 'showlegend': False, - } + { + "legendgroup": "", + "marker": { + "color": data_all[variable].to_numpy(), + "coloraxis": "coloraxis", + "symbol": "circle", + "size": point_size, + }, + "mode": "markers", + "name": "", + "showlegend": False, + } ], - } for variable in meta_num + } + for variable in meta_num ], "direction": "down", "showactive": True, "x": 1, "xanchor": "right", "y": 1, - "yanchor": "top" + "yanchor": "top", }, ], ) @@ -216,20 +233,24 @@ def figures_to_html(figs, filename="dashboard.html"): "label": variable, "method": "update", "args": [ - {"marker.color": config_button_cat(fig_cat, data_all, variable), - 'showlegend': False, #[data_all.shape[0]*[True]], - 'legendgroup': '',#data_all[variable], - 'name': '',#data_all[variable] + { + "marker.color": config_button_cat( + fig_cat, data_all, variable + ), + "showlegend": False, # [data_all.shape[0]*[True]], + "legendgroup": "", # data_all[variable], + "name": "", # data_all[variable] } ], - } for variable in meta_cat + } + for variable in meta_cat ], "direction": "down", "showactive": True, "x": 1, "xanchor": "right", "y": 1, - "yanchor": "top" + "yanchor": "top", }, ], ) diff --git a/workflow/scripts/plot_pca.R b/workflow/scripts/plot_pca.R index 4c40717..3a6c4cf 100644 --- a/workflow/scripts/plot_pca.R +++ b/workflow/scripts/plot_pca.R @@ -20,47 +20,47 @@ pairs_path <- snakemake@output[["pairs_plot"]] loadingsplot_path <- snakemake@output[["loadings_plot"]] loadings_lollipop_plot_path <- snakemake@output[["loadings_lollipop_plot"]] -pairs_size <- snakemake@config[["scatterplot2d"]][["size"]]/10 -pairs_alpha <- snakemake@config[["scatterplot2d"]][["alpha"]]/2 +pairs_size <- snakemake@config[["scatterplot2d"]][["size"]] / 10 +pairs_alpha <- snakemake@config[["scatterplot2d"]][["alpha"]] / 2 metadata_col <- c(snakemake@config[["metadata_of_interest"]])[1] ### load data -data <- data.frame(fread(file.path(data_path), header=TRUE), row.names=1) -metadata <- data.frame(fread(file.path(metadata_path), header=TRUE), row.names=1) +data <- data.frame(fread(file.path(data_path), header = TRUE), row.names = 1) +metadata <- data.frame(fread(file.path(metadata_path), header = TRUE), row.names = 1) # make rownames (R) syntactically valid rownames(data) <- make.names(rownames(data)) rownames(metadata) <- make.names(rownames(metadata)) # prepare metadata -if(is.null(metadata_col)){ - metadata_col <- colnames(metadata)[1] +if (is.null(metadata_col)) { + metadata_col <- colnames(metadata)[1] } # check if metadata column is only NA and switch to the first that is not -if(all(is.na(metadata[[metadata_col]]))){ - for(col in colnames(metadata)){ - if(all(is.na(metadata[[col]]))){ - next - }else{ - metadata_col <- col - break - } +if (all(is.na(metadata[[metadata_col]]))) { + for (col in colnames(metadata)) { + if (all(is.na(metadata[[col]]))) { + next + } else { + metadata_col <- col + break } + } } # make metadata rownames R "compatible" -rownames(metadata) <- gsub(pattern= '-' ,replacement = '.', x = rownames(metadata)) +rownames(metadata) <- gsub(pattern = "-", replacement = ".", x = rownames(metadata)) # align rows -data <- data[rownames(metadata),] +data <- data[rownames(metadata), ] -data_axes <- data.frame(fread(file.path(axes_path), header=TRUE), row.names=1) -data_loadings <- data.frame(fread(file.path(loadings_path), header=TRUE), row.names=1) +data_axes <- data.frame(fread(file.path(axes_path), header = TRUE), row.names = 1) +data_loadings <- data.frame(fread(file.path(loadings_path), header = TRUE), row.names = 1) -data_var <- data.frame(fread(file.path(var_path), header=TRUE), row.names=1) -colnames(data_var) <- c('var') -data_var$PC <- as.numeric(rownames(data_var))+1 +data_var <- data.frame(fread(file.path(var_path), header = TRUE), row.names = 1) +colnames(data_var) <- c("var") +data_var$PC <- as.numeric(rownames(data_var)) + 1 ### variance plot print("Variance plots") @@ -75,45 +75,49 @@ width_panel <- n_col * width height_panel <- 2 * height # number of top 10% of PCs -top_n <- ceiling(nrow(data_var)*0.1) +top_n <- ceiling(nrow(data_var) * 0.1) # make plots # options(repr.plot.width=width_panel, repr.plot.height=height_panel) pca_plots <- list() -pca_plots[["scree_all"]] <- ggplot(data_var, aes(x=PC,y=var, group=1))+ - geom_point(size=point_size)+ - geom_line(size=line_size)+ - xlab("Principal Component") + +pca_plots[["scree_all"]] <- ggplot(data_var, aes(x = PC, y = var, group = 1)) + + geom_point(size = point_size) + + geom_line(size = line_size) + + xlab("Principal Component") + ylab("Variance Explained") + - ggtitle("Scree Plot of all Principal Components")+ -theme_linedraw()+ theme(plot.title = element_text(size = 10)) - -pca_plots[["cum_all"]] <- ggplot(data_var, aes(x=PC,y=cumsum(var), group=1))+ - geom_point(size=point_size)+ - geom_line(size=line_size)+ - xlab("Principal Component") + + ggtitle("Scree Plot of all Principal Components") + + theme_linedraw() + + theme(plot.title = element_text(size = 10)) + +pca_plots[["cum_all"]] <- ggplot(data_var, aes(x = PC, y = cumsum(var), group = 1)) + + geom_point(size = point_size) + + geom_line(size = line_size) + + xlab("Principal Component") + ylab("Variance Explained") + - ggtitle("Cumulative Explained Variance of all Principal Components")+ -theme_linedraw()+ theme(plot.title = element_text(size = 10)) - -pca_plots[["scree_top"]] <- ggplot(data_var[1:top_n,], aes(x=PC,y=var, group=1))+ - geom_point(size=point_size)+ - geom_line(size=line_size)+ - xlab("Principal Component") + + ggtitle("Cumulative Explained Variance of all Principal Components") + + theme_linedraw() + + theme(plot.title = element_text(size = 10)) + +pca_plots[["scree_top"]] <- ggplot(data_var[1:top_n, ], aes(x = PC, y = var, group = 1)) + + geom_point(size = point_size) + + geom_line(size = line_size) + + xlab("Principal Component") + ylab("Variance Explained") + - ggtitle(paste0("Scree Plot of top ", top_n," Principal Components"))+ -theme_linedraw()+ theme(plot.title = element_text(size = 10)) + ggtitle(paste0("Scree Plot of top ", top_n, " Principal Components")) + + theme_linedraw() + + theme(plot.title = element_text(size = 10)) -pca_plots[["cum_top"]] <- ggplot(data_var[1:top_n,], aes(x=PC,y=cumsum(var), group=1))+ - geom_point(size=point_size)+ - geom_line(size=line_size)+ - xlab("Principal Component") + +pca_plots[["cum_top"]] <- ggplot(data_var[1:top_n, ], aes(x = PC, y = cumsum(var), group = 1)) + + geom_point(size = point_size) + + geom_line(size = line_size) + + xlab("Principal Component") + ylab("Variance Explained") + - ggtitle(paste0("Cumulative Explained Variance of top ", top_n," Principal Components"))+ -theme_linedraw()+ theme(plot.title = element_text(size = 10)) + ggtitle(paste0("Cumulative Explained Variance of top ", top_n, " Principal Components")) + + theme_linedraw() + + theme(plot.title = element_text(size = 10)) pca_plot_panel <- wrap_plots(pca_plots, ncol = n_col) @@ -122,135 +126,136 @@ pca_plot_panel <- wrap_plots(pca_plots, ncol = n_col) # pca_plot_panel ggsave(basename(diagnostics_path), - plot = pca_plot_panel, - device = 'png', - path = dirname(diagnostics_path), - scale = 1, - dpi = 300, - width = width_panel, - height = height_panel, - limitsize = FALSE, - ) + plot = pca_plot_panel, + device = "png", + path = dirname(diagnostics_path), + scale = 1, + dpi = 300, + width = width_panel, + height = height_panel, + limitsize = FALSE, +) ### pairs plot print("Pairs plot") # convert to categorical if less than 25 unique integer values -if (is.numeric(metadata[[metadata_col]]) & length(unique(metadata[[metadata_col]]))<=25){ - if(all(metadata[[metadata_col]] == round(metadata[[metadata_col]]))){ - metadata[metadata_col] <- as.factor(metadata[[metadata_col]]) - } +if (is.numeric(metadata[[metadata_col]]) & length(unique(metadata[[metadata_col]])) <= 25) { + if (all(metadata[[metadata_col]] == round(metadata[[metadata_col]]))) { + metadata[metadata_col] <- as.factor(metadata[[metadata_col]]) + } } # if a metadata class is empty ("") fill with "unknown" -if (!any(is.na(metadata[[metadata_col]]))){ - if (any(metadata[[metadata_col]]=="")){ - metadata[metadata[[metadata_col]]=="", metadata_col] <- "unknown" - } +if (!any(is.na(metadata[[metadata_col]]))) { + if (any(metadata[[metadata_col]] == "")) { + metadata[metadata[[metadata_col]] == "", metadata_col] <- "unknown" + } } # remove groups with less than 3 members from metadata and data # and set legend parameter according to data type legend <- NULL -if (!is.numeric(metadata[[metadata_col]])){ - - keep_groups <- names(table(metadata[[metadata_col]]))[table(metadata[[metadata_col]])>2] - keep_idx <- metadata[[metadata_col]] %in% keep_groups - metadata <- metadata[keep_idx,,drop=FALSE] - data <- data[rownames(metadata),] - - # only add legend in case of less than 10 groups within metadata - if(length(unique(metadata[[metadata_col]]))<11){ - legend <- 1 - } +if (!is.numeric(metadata[[metadata_col]])) { + keep_groups <- names(table(metadata[[metadata_col]]))[table(metadata[[metadata_col]]) > 2] + keep_idx <- metadata[[metadata_col]] %in% keep_groups + metadata <- metadata[keep_idx, , drop = FALSE] + data <- data[rownames(metadata), ] + + # only add legend in case of less than 10 groups within metadata + if (length(unique(metadata[[metadata_col]])) < 11) { + legend <- 1 + } } # check if one PC is only zeros (yes, that's apparently possible) -non_zero_cols <- unname(apply(data, 2, function(x) !all(x==0))) -data <- data[,non_zero_cols] -data_axes <- data_axes[non_zero_cols,,drop=FALSE] +non_zero_cols <- unname(apply(data, 2, function(x) !all(x == 0))) +data <- data[, non_zero_cols] +data_axes <- data_axes[non_zero_cols, , drop = FALSE] # make pairs plot -if(nrow(data)>0){ - # options(repr.plot.width=10, repr.plot.height=10) - n_dim <- min(10, ncol(data)) - - pairs_plot <- ggpairs( - data = data, - mapping = ggplot2::aes(color = metadata[[metadata_col]]), - columns = 1:n_dim, - title = paste0("PCA pairs plot colored by ",metadata_col), - upper = list(continuous = wrap("density", alpha = 0.5, size=0.25)), - lower = list(continuous = wrap("points", alpha = pairs_alpha, size = pairs_size)), - diag = list(continuous = wrap("densityDiag", alpha = 0.5, size=0.25)), - params = NULL, - xlab = NULL, - ylab = NULL, - axisLabels = c("show", "internal", "none"), - columnLabels = data_axes[1:n_dim,'label'], - labeller = "label_value", - switch = NULL, - showStrips = NULL, - legend = legend, - cardinality_threshold = 15, - progress = NULL, - proportions = NULL - )+ - theme(legend.position = "bottom") + +if (nrow(data) > 0) { + # options(repr.plot.width=10, repr.plot.height=10) + n_dim <- min(10, ncol(data)) + + pairs_plot <- ggpairs( + data = data, + mapping = ggplot2::aes(color = metadata[[metadata_col]]), + columns = 1:n_dim, + title = paste0("PCA pairs plot colored by ", metadata_col), + upper = list(continuous = wrap("density", alpha = 0.5, size = 0.25)), + lower = list(continuous = wrap("points", alpha = pairs_alpha, size = pairs_size)), + diag = list(continuous = wrap("densityDiag", alpha = 0.5, size = 0.25)), + params = NULL, + xlab = NULL, + ylab = NULL, + axisLabels = c("show", "internal", "none"), + columnLabels = data_axes[1:n_dim, "label"], + labeller = "label_value", + switch = NULL, + showStrips = NULL, + legend = legend, + cardinality_threshold = 15, + progress = NULL, + proportions = NULL + ) + + theme(legend.position = "bottom") + labs(fill = metadata_col) - if (is.numeric(metadata[[metadata_col]])){ - pairs_plot <- pairs_plot + scale_color_gradient2(midpoint=0, low="royalblue4", mid="grey80", high="firebrick2", space ="Lab") - } -}else{ - pairs_plot <- ggplot() + annotate("text", x = 0.5, y = 0.5, label = "No group with more than 2 members in the data.") + theme_void() - n_dim <- 5 + if (is.numeric(metadata[[metadata_col]])) { + pairs_plot <- pairs_plot + scale_color_gradient2(midpoint = 0, low = "royalblue4", mid = "grey80", high = "firebrick2", space = "Lab") + } +} else { + pairs_plot <- ggplot() + + annotate("text", x = 0.5, y = 0.5, label = "No group with more than 2 members in the data.") + + theme_void() + n_dim <- 5 } # save pairs plot ggsave(basename(pairs_path), - plot = pairs_plot, - device = 'png', - path = dirname(pairs_path), - scale = 1, - dpi = 300, - width = n_dim, - height = n_dim, - limitsize = FALSE, - ) + plot = pairs_plot, + device = "png", + path = dirname(pairs_path), + scale = 1, + dpi = 300, + width = n_dim, + height = n_dim, + limitsize = FALSE, +) ### loadings plot print("Loadings plot") # plot specifications -n_col <- min(5,n_dim) -height_panel <- ceiling(n_dim/n_col)*4 +n_col <- min(5, n_dim) +height_panel <- ceiling(n_dim / n_col) * 4 width_panel <- n_col * 4 loading_plots <- list() -for(i in 1:(n_dim-1)){ - tmp_x <- i - tmp_y <- i+1 - - # determine top 10 features per PC combination - top_features <- rownames(data_loadings)[order(-(data_loadings[paste0("PC_",tmp_x)]^2 + data_loadings[paste0("PC_",tmp_y)]^2))][1:10] - - # subset data - tmp_loadings <- data_loadings[top_features, c(paste0("PC_",tmp_x), paste0("PC_",tmp_y))] - tmp_loadings$features <- rownames(tmp_loadings) - text_var <- "features" - - # plot data - loading_plots[[i]] <- ggplot(data=tmp_loadings, aes_string(x=paste0("PC_",tmp_x), y=paste0("PC_",tmp_y), label=text_var))+ - geom_segment(data=tmp_loadings, aes_string(x=0, y=0, xend=paste0("PC_",tmp_x), yend=paste0("PC_",tmp_y)), arrow=arrow(length=unit(0.2,"cm")), alpha=0.75, color="black") + - geom_label_repel(size = 2)+ - xlab(paste0("Principal Component ",tmp_x)) + - ylab(paste0("Principal Component ",tmp_y)) + +for (i in 1:(n_dim - 1)) { + tmp_x <- i + tmp_y <- i + 1 + + # determine top 10 features per PC combination + top_features <- rownames(data_loadings)[order(-(data_loadings[paste0("PC_", tmp_x)]^2 + data_loadings[paste0("PC_", tmp_y)]^2))][1:10] + + # subset data + tmp_loadings <- data_loadings[top_features, c(paste0("PC_", tmp_x), paste0("PC_", tmp_y))] + tmp_loadings$features <- rownames(tmp_loadings) + text_var <- "features" + + # plot data + loading_plots[[i]] <- ggplot(data = tmp_loadings, aes_string(x = paste0("PC_", tmp_x), y = paste0("PC_", tmp_y), label = text_var)) + + geom_segment(data = tmp_loadings, aes_string(x = 0, y = 0, xend = paste0("PC_", tmp_x), yend = paste0("PC_", tmp_y)), arrow = arrow(length = unit(0.2, "cm")), alpha = 0.75, color = "black") + + geom_label_repel(size = 2) + + xlab(paste0("Principal Component ", tmp_x)) + + ylab(paste0("Principal Component ", tmp_y)) + theme_linedraw() - -# print(tmp_plot) + + # print(tmp_plot) } loadings_plot_panel <- wrap_plots(loading_plots, ncol = n_col) @@ -260,61 +265,58 @@ loadings_plot_panel <- wrap_plots(loading_plots, ncol = n_col) # loadings_plot_panel ggsave(basename(loadingsplot_path), - plot = loadings_plot_panel, - device = 'png', - path = dirname(loadingsplot_path), - scale = 1, - dpi = 300, - width = width_panel, - height = height_panel, - limitsize = FALSE, - ) + plot = loadings_plot_panel, + device = "png", + path = dirname(loadingsplot_path), + scale = 1, + dpi = 300, + width = width_panel, + height = height_panel, + limitsize = FALSE, +) ### loadings lollipop plot print("Loadings Lollipop plot") - + # plot specifications -n_col <- min(5,n_dim) -height_panel <- ceiling(n_dim/n_col)*2 +n_col <- min(5, n_dim) +height_panel <- ceiling(n_dim / n_col) * 2 width_panel <- n_col * 3 - + lollipops <- list() for (i in 1:n_dim) { - - # determine top 10 features per PC -# top_features <- rownames(data_loadings)[order(-(data_loadings[paste0("PC_",i)]^2))][1:10] - top_features <- rownames(data_loadings)[order(-abs(data_loadings[paste0("PC_", i)]))][1:10] - - # subset data - tmp_loadings <- data_loadings[top_features, paste0("PC_",i), drop=FALSE] - colnames(tmp_loadings) <- c("Loadings") - tmp_loadings$Features <- factor(rownames(tmp_loadings), levels=rev(rownames(tmp_loadings))) - - # make plot - lollipops[[i]] <- ggplot(tmp_loadings, aes(x=Loadings, y=Features)) + - geom_point(color="blue") + - geom_segment(aes(xend=0, yend=Features), color="black") + + # determine top 10 features per PC + # top_features <- rownames(data_loadings)[order(-(data_loadings[paste0("PC_",i)]^2))][1:10] + top_features <- rownames(data_loadings)[order(-abs(data_loadings[paste0("PC_", i)]))][1:10] + + # subset data + tmp_loadings <- data_loadings[top_features, paste0("PC_", i), drop = FALSE] + colnames(tmp_loadings) <- c("Loadings") + tmp_loadings$Features <- factor(rownames(tmp_loadings), levels = rev(rownames(tmp_loadings))) + + # make plot + lollipops[[i]] <- ggplot(tmp_loadings, aes(x = Loadings, y = Features)) + + geom_point(color = "blue") + + geom_segment(aes(xend = 0, yend = Features), color = "black") + theme_bw() + theme(plot.title = element_text(size = 10)) + - labs(title=paste("Principal Component ", i), x="Loading", y="Feature") + labs(title = paste("Principal Component ", i), x = "Loading", y = "Feature") } -lollipop_plot_panel <- wrap_plots(lollipops, ncol=n_col) +lollipop_plot_panel <- wrap_plots(lollipops, ncol = n_col) # save diangostics plot # options(repr.plot.width=width_panel, repr.plot.height=height_panel) # lollipop_plot_panel ggsave(basename(loadings_lollipop_plot_path), - plot = lollipop_plot_panel, - device = 'png', - path = dirname(loadings_lollipop_plot_path), - scale = 1, - dpi = 300, - width = width_panel, - height = height_panel, - limitsize = FALSE, - ) - - + plot = lollipop_plot_panel, + device = "png", + path = dirname(loadings_lollipop_plot_path), + scale = 1, + dpi = 300, + width = width_panel, + height = height_panel, + limitsize = FALSE, +) diff --git a/workflow/scripts/plot_umap_connectivity.py b/workflow/scripts/plot_umap_connectivity.py index 3d9da01..52f3445 100644 --- a/workflow/scripts/plot_umap_connectivity.py +++ b/workflow/scripts/plot_umap_connectivity.py @@ -2,10 +2,8 @@ #### libraries # general -import os import pickle -import pandas as pd -import matplotlib.pyplot as plt + # dimensionality reduction import umap import umap.plot @@ -19,10 +17,12 @@ plot_connectivity_path = snakemake.output["connectivity_plot"] ### load data -with open(object_path, 'rb') as f: +with open(object_path, "rb") as f: umap_obj = pickle.load(f) - + ### generate & save UMAP connectivity plot # umap.plot.connectivity(umap_obj, show_points=True) -umap.plot.connectivity(umap_obj, edge_bundling='hammer').figure.savefig(plot_connectivity_path) +umap.plot.connectivity(umap_obj, edge_bundling="hammer").figure.savefig( + plot_connectivity_path +) diff --git a/workflow/scripts/plot_umap_diagnostics.py b/workflow/scripts/plot_umap_diagnostics.py index 552a8de..225afba 100644 --- a/workflow/scripts/plot_umap_diagnostics.py +++ b/workflow/scripts/plot_umap_diagnostics.py @@ -2,10 +2,8 @@ #### libraries # general -import os import pickle -import pandas as pd -import matplotlib.pyplot as plt + # dimensionality reduction import umap import umap.plot @@ -19,14 +17,19 @@ plot_diagnostics_path = snakemake.output["diagnostics_plot"] ### load data -with open(object_path, 'rb') as f: +with open(object_path, "rb") as f: umap_obj = pickle.load(f) - + ### generate & save UMAP specific diagnostic plots -fig_diag, ax_diag = umap.plot.plt.subplots(2, 2, figsize=(12,12)) -umap.plot.diagnostic(umap_obj, diagnostic_type='pca', ax=ax_diag[0,0]) -umap.plot.diagnostic(umap_obj, diagnostic_type='vq', ax=ax_diag[0,1]) -umap.plot.diagnostic(umap_obj, diagnostic_type='local_dim', ax=ax_diag[1,0]) -umap.plot.diagnostic(umap_obj, diagnostic_type='neighborhood', nhood_size=min(umap_obj.n_neighbors, 15), ax=ax_diag[1,1]) +fig_diag, ax_diag = umap.plot.plt.subplots(2, 2, figsize=(12, 12)) +umap.plot.diagnostic(umap_obj, diagnostic_type="pca", ax=ax_diag[0, 0]) +umap.plot.diagnostic(umap_obj, diagnostic_type="vq", ax=ax_diag[0, 1]) +umap.plot.diagnostic(umap_obj, diagnostic_type="local_dim", ax=ax_diag[1, 0]) +umap.plot.diagnostic( + umap_obj, + diagnostic_type="neighborhood", + nhood_size=min(umap_obj.n_neighbors, 15), + ax=ax_diag[1, 1], +) fig_diag.savefig(plot_diagnostics_path) diff --git a/workflow/scripts/subset_data.py b/workflow/scripts/subset_data.py index 05fa4eb..6a8d48a 100644 --- a/workflow/scripts/subset_data.py +++ b/workflow/scripts/subset_data.py @@ -2,7 +2,6 @@ #### libraries # general -import os import pandas as pd #### configurations @@ -13,7 +12,7 @@ metadata_features_path = snakemake.output["metadata_features"] # parameters -samples_by_features = int(snakemake.params['samples_by_features']) +samples_by_features = int(snakemake.params["samples_by_features"]) features_to_plot = set(snakemake.params["features_to_plot"]) ### load data @@ -23,7 +22,7 @@ data = pd.read_csv(data_path, index_col=0) else: data = pd.read_csv(data_path, index_col=0).T - + ### check if "ALL" features should be plotted and overlap with columns & subset data if features_to_plot == {"ALL"}: features_to_plot = list(data.columns) @@ -31,11 +30,13 @@ features_to_plot = list(features_to_plot.intersection(set(data.columns))) # subset data -if len(features_to_plot)!=0: - data = data.loc[:,features_to_plot] +if len(features_to_plot) != 0: + data = data.loc[:, features_to_plot] else: - print("requested features to plot are not in the provided data, first 10 features will be plotted instead") - data = data.iloc[:,:10] + print( + "requested features to plot are not in the provided data, first 10 features will be plotted instead" + ) + data = data.iloc[:, :10] # save data -data.to_csv(metadata_features_path) \ No newline at end of file +data.to_csv(metadata_features_path) diff --git a/workflow/scripts/umap_embed.py b/workflow/scripts/umap_embed.py index f1754a7..d46f96b 100644 --- a/workflow/scripts/umap_embed.py +++ b/workflow/scripts/umap_embed.py @@ -5,28 +5,39 @@ import os import pickle import pandas as pd + # dimensionality reduction import umap #### configurations # ipnuts -data_path = snakemake.input[0] #"/nobackup/lab_bock/projects/macroIC/results/AKsmall/KOcall_NonTargeting/counts/CORRECTED_RNA.csv" -graph_object_path = snakemake.input[1] #snakemake.input["knn_object"] #"/nobackup/lab_bock/projects/macroIC/results/AKsmall/unsupervised_analysis/AKsmall_KOcall_NonTargeting_CORRECTED/UMAP_correlation_100_graph.pickle" +data_path = snakemake.input[ + 0 # "/nobackup/lab_bock/projects/macroIC/results/AKsmall/KOcall_NonTargeting/counts/CORRECTED_RNA.csv" +] +graph_object_path = snakemake.input[ + 1 # snakemake.input["knn_object"] #"/nobackup/lab_bock/projects/macroIC/results/AKsmall/unsupervised_analysis/AKsmall_KOcall_NonTargeting_CORRECTED/UMAP_correlation_100_graph.pickle" +] # outputs -result_object_path = snakemake.output["result_object"] #"/nobackup/lab_bock/projects/macroIC/results/AKsmall/unsupervised_analysis/AKsmall_KOcall_NonTargeting_CORRECTED/UMAP_correlation_100_0.1_2_object.pickle" -result_data_path = snakemake.output["result_data"] #"/nobackup/lab_bock/projects/macroIC/results/AKsmall/unsupervised_analysis/AKsmall_KOcall_NonTargeting_CORRECTED/UMAP_correlation_100_0.1_2_data.csv" -result_axes_path = snakemake.output["result_axes"] #"/nobackup/lab_bock/projects/macroIC/results/AKsmall/unsupervised_analysis/AKsmall_KOcall_NonTargeting_CORRECTED/UMAP_correlation_100_0.1_2_axes.csv" +result_object_path = snakemake.output[ + "result_object" # "/nobackup/lab_bock/projects/macroIC/results/AKsmall/unsupervised_analysis/AKsmall_KOcall_NonTargeting_CORRECTED/UMAP_correlation_100_0.1_2_object.pickle" +] +result_data_path = snakemake.output[ + "result_data" # "/nobackup/lab_bock/projects/macroIC/results/AKsmall/unsupervised_analysis/AKsmall_KOcall_NonTargeting_CORRECTED/UMAP_correlation_100_0.1_2_data.csv" +] +result_axes_path = snakemake.output[ + "result_axes" # "/nobackup/lab_bock/projects/macroIC/results/AKsmall/unsupervised_analysis/AKsmall_KOcall_NonTargeting_CORRECTED/UMAP_correlation_100_0.1_2_axes.csv" +] result_dir = os.path.dirname(result_object_path) # parameters -samples_by_features = int(snakemake.params['samples_by_features']) #0 -metric = snakemake.params['metric'] #"correlation" -n_neighbors = int(snakemake.params['n_neighbors']) #100 -min_dist = float(snakemake.params['min_dist']) #0.1 -n_components = int(snakemake.params['n_components']) #2 -densmap = True if int(snakemake.params['densmap'])==1 else False #False +samples_by_features = int(snakemake.params["samples_by_features"]) # 0 +metric = snakemake.params["metric"] # "correlation" +n_neighbors = int(snakemake.params["n_neighbors"]) # 100 +min_dist = float(snakemake.params["min_dist"]) # 0.1 +n_components = int(snakemake.params["n_components"]) # 2 +densmap = True if int(snakemake.params["densmap"]) == 1 else False # False # make directory if not existing if not os.path.exists(result_dir): @@ -41,80 +52,83 @@ data = pd.read_csv(data_path, index_col=0).T # if less than 11 datapoints there is no pre-computed KNN graph -if data.shape[0]<11: +if data.shape[0] < 11: print("no pre-computed KNN graph will be used") knn = (None, None, None) else: # load pre-computed KNN graph - with open(graph_object_path, 'rb') as f: + with open(graph_object_path, "rb") as f: # The protocol version used is detected automatically, so we do not # have to specify it. knn = pickle.load(f) ### embed data in low dimensions -umap_obj=umap.umap_.UMAP(n_neighbors=n_neighbors, - n_components=n_components, - metric=metric, - metric_kwds=None, - output_metric='euclidean', - output_metric_kwds=None, - n_epochs=None, - learning_rate=1.0, - init='spectral', - min_dist=min_dist, - spread=1.0, - low_memory=True, - n_jobs=-1, - set_op_mix_ratio=1.0, - local_connectivity=1.0, - repulsion_strength=1.0, - negative_sample_rate=5, - transform_queue_size=4.0, - a=None, - b=None, - random_state=42, - angular_rp_forest=False, - target_n_neighbors=-1, - target_metric='categorical', - target_metric_kwds=None, - target_weight=0.5, - transform_seed=42, - transform_mode='embedding', - force_approximation_algorithm=False, - verbose=False, - tqdm_kwds=None, - unique=False, - densmap=densmap, - dens_lambda=2.0, - dens_frac=0.3, - dens_var_shift=0.1, - output_dens=False, - disconnection_distance=None, - precomputed_knn=knn - ).fit(data) +umap_obj = umap.umap_.UMAP( + n_neighbors=n_neighbors, + n_components=n_components, + metric=metric, + metric_kwds=None, + output_metric="euclidean", + output_metric_kwds=None, + n_epochs=None, + learning_rate=1.0, + init="spectral", + min_dist=min_dist, + spread=1.0, + low_memory=True, + n_jobs=-1, + set_op_mix_ratio=1.0, + local_connectivity=1.0, + repulsion_strength=1.0, + negative_sample_rate=5, + transform_queue_size=4.0, + a=None, + b=None, + random_state=42, + angular_rp_forest=False, + target_n_neighbors=-1, + target_metric="categorical", + target_metric_kwds=None, + target_weight=0.5, + transform_seed=42, + transform_mode="embedding", + force_approximation_algorithm=False, + verbose=False, + tqdm_kwds=None, + unique=False, + densmap=densmap, + dens_lambda=2.0, + dens_frac=0.3, + dens_var_shift=0.1, + output_dens=False, + disconnection_distance=None, + precomputed_knn=knn, +).fit(data) # data_embedding = umap_obj.fit_transform(data) -data_df = pd.DataFrame(umap_obj.embedding_, index=data.index,) +data_df = pd.DataFrame( + umap_obj.embedding_, + index=data.index, +) data_df = data_df.rename_axis(("sample_name")) if densmap: - data_df.columns = ["densMAP_{}".format(str(idx+1)) for idx in data_df.columns] + data_df.columns = ["densMAP_{}".format(str(idx + 1)) for idx in data_df.columns] else: - data_df.columns = ["UMAP_{}".format(str(idx+1)) for idx in data_df.columns] + data_df.columns = ["UMAP_{}".format(str(idx + 1)) for idx in data_df.columns] ### save data # save umap object -with open(result_object_path, 'wb') as f: +with open(result_object_path, "wb") as f: pickle.dump(umap_obj, f, pickle.HIGHEST_PROTOCOL) - + # save transformed data data_df.to_csv(result_data_path) # save axes information for visualization axes_info_df = pd.DataFrame(data_df.columns) -axes_info_df.columns = ['label'] -axes_info_df['label'] = [label.replace("_", "")for label in axes_info_df['label']] +axes_info_df.columns = ["label"] +axes_info_df["label"] = [label.replace("_", "") for label in axes_info_df["label"]] axes_info_df.to_csv(result_axes_path) - diff --git a/workflow/scripts/umap_graph.py b/workflow/scripts/umap_graph.py index a58ea70..aaa58af 100644 --- a/workflow/scripts/umap_graph.py +++ b/workflow/scripts/umap_graph.py @@ -5,23 +5,27 @@ import os import pickle import pandas as pd + # dimensionality reduction -import umap from umap.umap_ import nearest_neighbors #### configurations # ipnuts -data_path = snakemake.input["data"] #"/nobackup/lab_bock/projects/macroIC/results/CC001/counts/mdm_normCQN.csv" +data_path = snakemake.input[ + "data" # "/nobackup/lab_bock/projects/macroIC/results/CC001/counts/mdm_normCQN.csv" +] # outputs -result_object_path = snakemake.output["result_object"] # "/nobackup/lab_bock/projects/macroIC/results/CC001/unsupervised_analysis/mdm_normCQN/UMAP/UMAP_correlation_5_0.1_2_object.pickle" +result_object_path = snakemake.output[ + "result_object" # "/nobackup/lab_bock/projects/macroIC/results/CC001/unsupervised_analysis/mdm_normCQN/UMAP/UMAP_correlation_5_0.1_2_object.pickle" +] result_dir = os.path.dirname(result_object_path) # parameters -samples_by_features = int(snakemake.params['samples_by_features']) #0 -metric = snakemake.params['metric'] # "correlation" -n_neighbors = int(snakemake.params['n_neighbors']) #100 +samples_by_features = int(snakemake.params["samples_by_features"]) # 0 +metric = snakemake.params["metric"] # "correlation" +n_neighbors = int(snakemake.params["n_neighbors"]) # 100 # make directory if not existing if not os.path.exists(result_dir): @@ -36,25 +40,28 @@ data = pd.read_csv(data_path, index_col=0).T # if less than 11 datapoints the KNN graph object can not be serialized (PyNNdescent issue: https://github.com/Teichlab/bbknn/issues/48) -if data.shape[0]<11: +if data.shape[0] < 11: from pathlib import Path + Path(result_object_path).touch() import sys + sys.exit() - + ### get knn-graph -knn = nearest_neighbors(data, - n_neighbors=n_neighbors, - metric=metric, - metric_kwds=None, - angular=False, - random_state=42, - low_memory=True, - use_pynndescent=True, - n_jobs=-1, - verbose=False - ) +knn = nearest_neighbors( + data, + n_neighbors=n_neighbors, + metric=metric, + metric_kwds=None, + angular=False, + random_state=42, + low_memory=True, + use_pynndescent=True, + n_jobs=-1, + verbose=False, +) ## save knn graph object -with open(result_object_path, 'wb') as f: - pickle.dump(knn, f, pickle.HIGHEST_PROTOCOL) \ No newline at end of file +with open(result_object_path, "wb") as f: + pickle.dump(knn, f, pickle.HIGHEST_PROTOCOL) diff --git a/workflow/scripts/validation_external.py b/workflow/scripts/validation_external.py index ca3b7dd..a1f620b 100644 --- a/workflow/scripts/validation_external.py +++ b/workflow/scripts/validation_external.py @@ -17,36 +17,41 @@ # load the clustering results and categorical metadata metadata = pd.read_csv(metadata_path, index_col=0) -metadata = metadata.fillna('') +metadata = metadata.fillna("") clustering_results = pd.read_csv(clusterings_path, index_col=0) -indices = [s.split('external_index_')[1].split('.csv')[0] for s in result_paths] +indices = [s.split("external_index_")[1].split(".csv")[0] for s in result_paths] # identify categorical metadata meta_cat = [] for variable in metadata.columns: unique_vals = list(metadata[variable].unique()) - + # check if integer AND less than 25 unique values -> categorical metadata - if all([isinstance(i, (int, np.int64)) for i in unique_vals]) and len(unique_vals)<25: + if ( + all([isinstance(i, (int, np.int64)) for i in unique_vals]) + and len(unique_vals) < 25 + ): meta_cat.append(variable) continue - + if all([isinstance(i, (str, bool, np.bool_)) for i in unique_vals]): meta_cat.append(variable) - + else: print("variable {} not categorical".format(variable)) # subset for categorical data -categorical_metadata = metadata.loc[:,meta_cat] +categorical_metadata = metadata.loc[:, meta_cat] # Ensure that the clustering results and categorical metadata have the same indices # fix metadata indices if they do not agree with data as they come from outside the workflow (e.g., R) -if not(set(clustering_results.index) == set(categorical_metadata.index)): - if categorical_metadata.index.inferred_type=='string': - categorical_metadata.index = [idx.replace('-','.') for idx in categorical_metadata.index] +if not (set(clustering_results.index) == set(categorical_metadata.index)): + if categorical_metadata.index.inferred_type == "string": + categorical_metadata.index = [ + idx.replace("-", ".") for idx in categorical_metadata.index + ] # Reorder the rows of the clustering results DataFrame clustering_results = clustering_results.reindex(categorical_metadata.index) @@ -55,23 +60,35 @@ # Create a DataFrame for each index for idx in indices: - idx_dfs[idx] = pd.DataFrame(index=clustering_results.columns, columns=categorical_metadata.columns) - + idx_dfs[idx] = pd.DataFrame( + index=clustering_results.columns, columns=categorical_metadata.columns + ) + # For each clustering result for clustering in clustering_results.columns: # For each categorical metadata for metadata in categorical_metadata.columns: # Calculate & store the scores - idx_dfs["AMI"].loc[clustering, metadata] = metrics.adjusted_mutual_info_score(categorical_metadata[metadata], clustering_results[clustering]) - idx_dfs["ARI"].loc[clustering, metadata] = metrics.adjusted_rand_score(categorical_metadata[metadata], clustering_results[clustering]) - idx_dfs["FMI"].loc[clustering, metadata] = metrics.fowlkes_mallows_score(categorical_metadata[metadata], clustering_results[clustering]) - - homogeneity, completeness, v_measure = metrics.homogeneity_completeness_v_measure(categorical_metadata[metadata], clustering_results[clustering]) + idx_dfs["AMI"].loc[clustering, metadata] = metrics.adjusted_mutual_info_score( + categorical_metadata[metadata], clustering_results[clustering] + ) + idx_dfs["ARI"].loc[clustering, metadata] = metrics.adjusted_rand_score( + categorical_metadata[metadata], clustering_results[clustering] + ) + idx_dfs["FMI"].loc[clustering, metadata] = metrics.fowlkes_mallows_score( + categorical_metadata[metadata], clustering_results[clustering] + ) + + homogeneity, completeness, v_measure = ( + metrics.homogeneity_completeness_v_measure( + categorical_metadata[metadata], clustering_results[clustering] + ) + ) idx_dfs["Homogeneity"].loc[clustering, metadata] = homogeneity idx_dfs["Completeness"].loc[clustering, metadata] = completeness idx_dfs["V"].loc[clustering, metadata] = v_measure - + # # Save the DataFrames as CSV files for i, idx in enumerate(indices): idx_dfs[idx].to_csv(result_paths[i]) diff --git a/workflow/scripts/validation_internal.R b/workflow/scripts/validation_internal.R index a91a57e..798f0a3 100644 --- a/workflow/scripts/validation_internal.R +++ b/workflow/scripts/validation_internal.R @@ -8,14 +8,14 @@ set.seed(42) # helper function for BIC calculation do_BIC <- function(x) { - # Check if there are at least two unique values in the column - if(length(unique(x)) < 2){ - # Return Inf if there's only one unique value - return(Inf) - } else { - # Perform BIC calculation if there are two or more unique values - return(BIC(lm(data_mtx[,i] ~ as.factor(x)))) - } + # Check if there are at least two unique values in the column + if (length(unique(x)) < 2) { + # Return Inf if there's only one unique value + return(Inf) + } else { + # Perform BIC calculation if there are two or more unique values + return(BIC(lm(data_mtx[, i] ~ as.factor(x)))) + } } ### configurations @@ -30,15 +30,15 @@ pca_var_path <- snakemake@input[["pca_var"]] result_path <- snakemake@output[["internal_indices"]] # parameters -internal_index <- as.character(snakemake@params['internal_index']) #"Silhouette" -sample_proportion <- as.numeric(snakemake@params['sample_proportion']) #0.1 -metadata_of_interest <- unlist(c(snakemake@params['metadata_of_interest'])) +internal_index <- as.character(snakemake@params["internal_index"]) # "Silhouette" +sample_proportion <- as.numeric(snakemake@params["sample_proportion"]) # 0.1 +metadata_of_interest <- unlist(c(snakemake@params["metadata_of_interest"])) ### load data -metadata <- data.frame(fread(file.path(metadata_path), header=TRUE), row.names=1) -clusterings <- data.frame(fread(file.path(clusterings_path), header=TRUE), row.names=1) -pca_var <- data.frame(fread(file.path(pca_var_path), header=TRUE), row.names=1) -pca <- data.frame(fread(file.path(pca_path), header=TRUE), row.names=1) +metadata <- data.frame(fread(file.path(metadata_path), header = TRUE), row.names = 1) +clusterings <- data.frame(fread(file.path(clusterings_path), header = TRUE), row.names = 1) +pca_var <- data.frame(fread(file.path(pca_var_path), header = TRUE), row.names = 1) +pca <- data.frame(fread(file.path(pca_path), header = TRUE), row.names = 1) ### USED BEFORE when PCA was not configurable # # load PCs that explain >90% of the variance in the data @@ -53,37 +53,37 @@ pca <- data.frame(fread(file.path(pca_path), header=TRUE), row.names=1) # pca <- data.frame(fread(file.path(pca_path), colClasses = classes, header=TRUE), row.names=1) # subset metadata to metadata_of_interest -if(length(metadata_of_interest)==0){ - metadata <- metadata[,1,drop=FALSE] -}else{ - metadata <- metadata[,metadata_of_interest,drop=FALSE] +if (length(metadata_of_interest) == 0) { + metadata <- metadata[, 1, drop = FALSE] +} else { + metadata <- metadata[, metadata_of_interest, drop = FALSE] } # transform metadata na_cols <- c() -for (col in colnames(metadata)){ - # if NA or less than 2 unique values -> remove column and move on - if (any(is.na(metadata[[col]])) | length(unique(metadata[[col]]))<2){ - na_cols <- c(na_cols, col) - next - } - # if a metadata class is empty ("") fill with "unknown" - if (!any(is.na(metadata[[col]]))){ - if (any(metadata[[col]]=="")){ - metadata[metadata[[col]]=="", col] <- "unknown" - } +for (col in colnames(metadata)) { + # if NA or less than 2 unique values -> remove column and move on + if (any(is.na(metadata[[col]])) | length(unique(metadata[[col]])) < 2) { + na_cols <- c(na_cols, col) + next + } + # if a metadata class is empty ("") fill with "unknown" + if (!any(is.na(metadata[[col]]))) { + if (any(metadata[[col]] == "")) { + metadata[metadata[[col]] == "", col] <- "unknown" } - # convert metadata to categorical if less than 25 unique integer values - if (is.numeric(metadata[[col]]) & length(unique(metadata[[col]]))<=25){ - if(all(metadata[[col]] == round(metadata[[col]]))){ - metadata[col] <- as.factor(metadata[[col]]) - } + } + # convert metadata to categorical if less than 25 unique integer values + if (is.numeric(metadata[[col]]) & length(unique(metadata[[col]])) <= 25) { + if (all(metadata[[col]] == round(metadata[[col]]))) { + metadata[col] <- as.factor(metadata[[col]]) } + } } # remove columns with NA -metadata <- metadata[, !(colnames(metadata) %in% na_cols),drop=FALSE] +metadata <- metadata[, !(colnames(metadata) %in% na_cols), drop = FALSE] # add categorical metadata to clustering results with prefix "metadata_" -metadata_cat <- metadata[,sapply(metadata, function(x) !is.numeric(x)), drop=FALSE] +metadata_cat <- metadata[, sapply(metadata, function(x) !is.numeric(x)), drop = FALSE] # Convert all categorical columns to integer metadata_cat[colnames(metadata_cat)] <- lapply(metadata_cat[colnames(metadata_cat)], function(x) as.integer(factor(x))) colnames(metadata_cat) <- paste0("metadata_", colnames(metadata_cat)) @@ -97,32 +97,32 @@ indices_df <- data.frame(matrix(ncol = 1, nrow = ncol(clusterings), dimnames = l clusterings[colnames(clusterings)] <- lapply(clusterings[colnames(clusterings)], function(x) as.integer(x)) data_mtx <- as.matrix(pca) data_mtx <- data_mtx[sample(nrow(data_mtx), ceiling(sample_proportion * nrow(data_mtx))), ] -clusterings <- clusterings[rownames(data_mtx),,drop=FALSE] +clusterings <- clusterings[rownames(data_mtx), , drop = FALSE] # calculate internal cluster index -if(internal_index %in% c("Silhouette", "Calinski_Harabasz", "C_index", "Davies_Bouldin", "Dunn")){ - for(clust in colnames(clusterings)){ - indices_df[clust,internal_index] <- intCriteria(traj=data_mtx, part=clusterings[[clust]], crit=c(internal_index)) - } -} else if(internal_index=="AIC"){ # -> NOT USED - ### determine indices using AIC on top PC of PCA - AIC_sum <- rep(0L, ncol(clusterings)) - - for(i in 1:ncol(data_mtx)){ - AIC_sum <- AIC_sum + unlist(apply(clusterings,2,function(x) AIC(lm(data_mtx[,i]~as.factor(x)))))*pca_var[i,1] - } - indices_df$AIC <- AIC_sum -} else if(internal_index=="BIC"){ - ### determine indices using BIC on top PC of PCA - BIC_sum <- rep(0L, ncol(clusterings)) - - for(i in 1:ncol(pca)){ -# BIC_sum <- BIC_sum + unlist(apply(clusterings,2,function(x) BIC(lm(data_mtx[,i]~as.factor(x)))))*pca_var[i,1] # crashed in case of only 1 cluster - BIC_sum <- BIC_sum + unlist(apply(clusterings, 2, do_BIC)) * pca_var[i, 1] - } - indices_df$BIC <- BIC_sum +if (internal_index %in% c("Silhouette", "Calinski_Harabasz", "C_index", "Davies_Bouldin", "Dunn")) { + for (clust in colnames(clusterings)) { + indices_df[clust, internal_index] <- intCriteria(traj = data_mtx, part = clusterings[[clust]], crit = c(internal_index)) + } +} else if (internal_index == "AIC") { # -> NOT USED + ### determine indices using AIC on top PC of PCA + AIC_sum <- rep(0L, ncol(clusterings)) + + for (i in 1:ncol(data_mtx)) { + AIC_sum <- AIC_sum + unlist(apply(clusterings, 2, function(x) AIC(lm(data_mtx[, i] ~ as.factor(x))))) * pca_var[i, 1] + } + indices_df$AIC <- AIC_sum +} else if (internal_index == "BIC") { + ### determine indices using BIC on top PC of PCA + BIC_sum <- rep(0L, ncol(clusterings)) + + for (i in 1:ncol(pca)) { + # BIC_sum <- BIC_sum + unlist(apply(clusterings,2,function(x) BIC(lm(data_mtx[,i]~as.factor(x)))))*pca_var[i,1] # crashed in case of only 1 cluster + BIC_sum <- BIC_sum + unlist(apply(clusterings, 2, do_BIC)) * pca_var[i, 1] + } + indices_df$BIC <- BIC_sum } ### save results # write.csv(indices_df, file=result_path, row.names=TRUE) -fwrite(as.data.frame(indices_df), file=file.path(result_path), row.names=TRUE) \ No newline at end of file +fwrite(as.data.frame(indices_df), file = file.path(result_path), row.names = TRUE)