Full analysis (#23)

* Save analysis

* More cleanup

* Fixed

* Linting

---------

Co-authored-by: Aaron Meyer <git@asmlab.org>

Author: aarmey

pf2barcode/figures/figure4.py

@@ -1,19 +1,20 @@
 from pathlib import Path
 
-import anndata
-import hdf5plugin  # noqa: F401
 import numpy as np
 import pandas as pd
 import scanpy as sc
+from anndata import AnnData, concat
+from anndata.io import read_text
 from scipy.sparse import csr_array, csr_matrix
 from scipy.special import xlogy
+from sklearn.preprocessing import scale
 from sklearn.utils.sparsefuncs import (
     inplace_column_scale,
     mean_variance_axis,
 )
 
 
-def prepare_dataset(X: anndata.AnnData, geneThreshold: float) -> anndata.AnnData:
+def prepare_dataset(X: AnnData, geneThreshold: float) -> AnnData:
     assert isinstance(X.X, csr_matrix)
     assert np.amin(X.X.data) >= 0.0
 
@@ -37,7 +38,7 @@ def prepare_dataset(X: anndata.AnnData, geneThreshold: float) -> anndata.AnnData
     return X
 
 
-def prepare_dataset_dev(X: anndata.AnnData) -> anndata.AnnData:
+def prepare_dataset_dev(X: AnnData) -> AnnData:
     X.X = csr_array(X.X)  # type: ignore
     assert np.amin(X.X.data) >= 0.0
 
@@ -88,7 +89,7 @@ def prepare_dataset_dev(X: anndata.AnnData) -> anndata.AnnData:
     return X
 
 
-def import_CCLE(pca_option="dev_pca") -> anndata.AnnData:
+def import_CCLE(pca_option="dev_pca", n_comp=10) -> AnnData:
     # pca option should be passed as either pca or glm_pca
     """Imports barcoded cell data."""
     adatas = {}
@@ -103,7 +104,7 @@ def import_CCLE(pca_option="dev_pca") -> anndata.AnnData:
         # "T1_MDAMB231",
         "T2_MDAMB231",
     ):
-        data = anndata.read_text(current_dir / "data" / f"{name}_count_mtx.tsv.bz2").T
+        data = read_text(current_dir / "data" / f"{name}_count_mtx.tsv.bz2").T
         barcodes = pd.read_csv(
             current_dir / "data" / f"{name}_SW.txt", sep="\t", index_col=0, header=0
         )
@@ -115,7 +116,7 @@ def import_CCLE(pca_option="dev_pca") -> anndata.AnnData:
         barcode_dfs.append(barcodes)
         adatas[name] = data
 
-    X = anndata.concat(adatas, label="sample", index_unique="-")
+    X = concat(adatas, label="sample", index_unique="-")
     X.X = csr_matrix(X.X)
 
     counts = X.obs["SW"].value_counts()
@@ -132,9 +133,11 @@ def import_CCLE(pca_option="dev_pca") -> anndata.AnnData:
     # conditional statement for either dev_pca or pca
     if pca_option == "dev_pca":
         X = prepare_dataset_dev(X)
-        sc.pp.pca(X, n_comps=20, svd_solver="arpack")
+        X.X = scale(X.X)
+        sc.pp.pca(X, n_comps=n_comp, svd_solver="arpack")
     else:
         X = prepare_dataset(X, geneThreshold=0.001)
-        sc.pp.pca(X, n_comps=20, svd_solver="arpack")
+        X.X = scale(X.X)
+        sc.pp.pca(X, n_comps=n_comp, svd_solver="arpack")
 
     return X

pf2barcode/figures/figure6.py

@@ -8,19 +8,15 @@ authors = [
 license = { text = "MIT" }
 requires-python = ">=3.12,<3.13"
 dependencies = [
-    "numpy>=2.0",
-    "scipy>=1.14",
-    "scikit-learn>=1.5",
+    "numpy>=2.2",
+    "scipy==1.15.3",
+    "scikit-learn>=1.7",
     "seaborn>=0.13",
-    "pandas>=2.2",
+    "pandas>=2.3",
     "gseapy>=1.1",
-    "scanpy>=1.10",
-    "dask[dataframe]>=2024.3.1",
-    "ipykernel>=6.29.4",
-    "setuptools>=74.0.0",
-    "ipywidgets>=8.1.3",
+    "scanpy>=1.11.2",
+    "ipykernel>=6.29",
     "anndata>=0.11.4",
-    "hdf5plugin>=5.1.0",
 ]
 
 [tool.rye]

pf2barcode/imports.py

@@ -0,0 +1,51 @@
+---
+title: "Analysis"
+format: html
+jupyter: python3
+---
+
+Generates a bar plot visualizing the relationship of PCs and computed
+negative log10 p-values from the Kruskal-Wallis H-test
+
+Computed p-values determines if distributions of PCs are statistically
+significantly across different groups, and the negative log10 transformation
+of the p-values allows for easier identification and interpretation of signficant PCs
+
+```{python}
+import scanpy as sc
+from pf2barcode.imports import import_CCLE
+from pf2barcode.analysis import anova_pvalues, kruskal_pvalues
+
+import matplotlib.pyplot as plt
+import numpy as np
+import seaborn as sns
+
+adata = import_CCLE("dev_pca")
+
+# Implement anova_pvalues function
+pvalues = anova_pvalues(adata)
+
+# Implement kruskal_pvalues function
+# pvalues = kruskal_pvalues(adata)
+```
+
+```{python}
+# Barplot setup
+sns.barplot(x=np.arange(pvalues.shape[0]), y=-np.log10(pvalues))
+plt.xlabel("PC")
+plt.ylabel("-log10(p-value)")
+```
+
+```{python}
+sc.pl.pca_variance_ratio(adata)
+```
+
+
+```{python}
+sc.pl.pca(adata, color="SW", components="2,3")
+```
+
+
+```{python}
+sc.pl.pca_loadings(adata, components="1,2")
+```