Feature/output names

docker/spatial_py/requirements.txt

@@ -1,13 +1,13 @@
 argparse==1.4.0
-anndata==0.10.9
-numpy==2.0.1
-pandas==2.2.2
-scipy==1.15.2
+anndata==0.11.4
+numpy==2.2.0
+pandas==2.3.1
+scipy==1.15.3
 squidpy==1.6.2
 matplotlib==3.10.0
 seaborn==0.13.2
 harmonypy==0.0.10
-scanpy==1.10.4
+scanpy==1.11.3
 dask==2024.11.2
 distributed==2024.11.2
 numcodecs==0.11.0

docker/spatial_py/scripts/geomx_export_final_artifacts

@@ -0,0 +1,47 @@
+#!/usr/bin/env python3
+
+import argparse
+import scanpy as sc
+
+
+def main(args):
+    ############################
+    ## EXPORT FINAL ARTIFACTS ##
+    ############################
+    adata = sc.read_h5ad(args.adata_input)
+
+    # Save outputs
+    metadata = adata.obs
+    metadata.to_csv(f"{args.cohort_id}.final_metadata.csv")
+    adata.write_h5ad(filename=args.adata_output, compression="gzip")
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description="Export final AnnData object and its metadata"
+    )
+    parser.add_argument(
+        "-c",
+        "--cohort-id",
+        type=str,
+        required=True,
+        help="Cohort ID"
+    )
+    parser.add_argument(
+        "-i",
+        "--adata-input",
+        type=str,
+        required=True,
+        help="Adata object to finalize"
+    )
+    parser.add_argument(
+        "-o",
+        "--adata-output",
+        type=str,
+        required=True,
+        help="Output file name for the final AnnData object"
+    )
+
+    args = parser.parse_args()
+
+    main(args)

docker/spatial_py/scripts/geomx_merge_and_prep

@@ -24,16 +24,31 @@ def main(args):
     ##################
     ## ANNOTATE HVG ##
     ##################
+    merged_adata.layers["counts"] = merged_adata.X.copy()
     sc.pp.log1p(merged_adata)
 
     sc.pp.highly_variable_genes(
         merged_adata,
-        flavor="seurat",
         n_top_genes=args.n_top_genes,
+        flavor="seurat",
         inplace=True,
+        batch_key=args.batch_key,
+    )
+    full_features = merged_adata.var.copy()
+    hvg_full = (
+        merged_adata.var[merged_adata.var["highly_variable"]]
+        .sort_values(
+            by=["highly_variable_nbatches", "dispersions_norm"],
+            ascending=[False, False]
+        )
+        .head(args.n_top_genes)
+        .index.tolist()
     )
+    hvg_adata = merged_adata[:, merged_adata.var.index.isin(hvg_full)]
+    hvg_features = hvg_adata.var.copy()
+
     sc.pl.highly_variable_genes(
-        merged_adata,
+        hvg_adata,
     )
     fig = plt.gcf()
     fig.suptitle(f"Highly variable genes dispersion plot - {args.output_prefix}", va="center", ha="center", fontsize=16)
@@ -44,20 +59,22 @@ def main(args):
     ## DIMENSIONALITY REDUCTION ##
     ##############################
     sc.pp.pca(
-        merged_adata,
+        hvg_adata,
         n_comps=args.n_comps,
         svd_solver="arpack",
     )
 
     # Save outputs
     metadata = merged_adata.obs
-    metadata.to_csv(f"{args.output_prefix}.merged_adata_metadata.csv")
-    merged_adata.write_h5ad(filename=args.adata_output, compression="gzip")
+    metadata.to_csv(f"{args.output_prefix}.merged_metadata.csv")
+    full_features.to_csv(f"{args.output_prefix}.all_genes.csv", index=True)
+    hvg_features.to_csv(f"{args.output_prefix}.hvg_genes.csv", index=True)
+    hvg_adata.write_h5ad(filename=args.adata_output, compression="gzip")
 
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(
-        description="Merge adata objects and process to prepare for downstream analysis"
+        description="Merge adata objects and process to prepare for downstream analysis by identifying highly variable genes (HVGs) and PCA"
     )
     parser.add_argument(
         "-i",
@@ -80,6 +97,13 @@ if __name__ == "__main__":
         required=True,
         help="Number of principal components to compute using `scanpy.pp.pca` [30]"
     )
+    parser.add_argument(
+        "-b",
+        "--batch-key",
+        type=str,
+        required=True,
+        help="Key in AnnData object for batch information ['batch_id']"
+    )
     parser.add_argument(
         "-p",
         "--output-prefix",

docker/spatial_py/scripts/visium_merge_and_plot_qc

@@ -41,39 +41,41 @@ def main(args):
         multi_panel=True,
     )
     fig = plt.gcf()
-    fig.suptitle(f"QC violin plot - {args.qc_plots_prefix}", va="center", ha="center", fontsize=16)
-    plt.savefig(f"{args.qc_plots_prefix}.qc_violin.png", dpi=300, bbox_inches="tight")
+    fig.suptitle(f"QC violin plot - {args.output_prefix}", va="center", ha="center", fontsize=16)
+    plt.savefig(f"{args.output_prefix}.qc_violin.png", dpi=300, bbox_inches="tight")
 
     # Total counts and n genes by counts distribution plots
     fig, axs = plt.subplots(1, 4, figsize=(15,4))
-    fig.suptitle(f"Covariates for filtering - {args.qc_plots_prefix}")
-    sns.distplot(
+    fig.suptitle(f"Covariates for filtering - {args.output_prefix}")
+    sns.histplot(
         merged_adata.obs["total_counts"],
         kde=False,
         ax=axs[0],
     )
-    sns.distplot(
+    sns.histplot(
         merged_adata.obs["total_counts"][merged_adata.obs["total_counts"]<10000],
         kde=False,
         bins=40,
         ax=axs[1],
     )
-    sns.distplot(
+    sns.histplot(
         merged_adata.obs["n_genes_by_counts"],
         kde=False,
         bins=60,
         ax=axs[2],
     )
-    sns.distplot(
+    sns.histplot(
         merged_adata.obs["n_genes_by_counts"][merged_adata.obs["n_genes_by_counts"]<4000],
         kde=False,
         bins=60,
         ax=axs[3],
     )
 
-    plt.savefig(f"{args.qc_plots_prefix}.qc_dist.png", dpi=300, bbox_inches="tight")
+    plt.savefig(f"{args.output_prefix}.qc_dist.png", dpi=300, bbox_inches="tight")
 
     # Save outputs
+    metadata = merged_adata.obs
+    metadata.to_csv(f"{args.output_prefix}.merged_metadata.csv")
     merged_adata.write_h5ad(filename=args.merged_adata_output, compression="gzip")
 
 
@@ -90,10 +92,10 @@ if __name__ == "__main__":
     )
     parser.add_argument(
         "-p",
-        "--qc-plots-prefix",
+        "--output-prefix",
         type=str,
         required=True,
-        help="Output file name prefix for the QC violin and distribution plots"
+        help="Output file name prefix for the QC violin, distribution plots, and metadata"
     )
     parser.add_argument(
         "-o",

docker/spatial_py/scripts/visium_process

@@ -43,29 +43,45 @@ def main(args):
     ##################
     sc.pp.highly_variable_genes(
         adata,
-        flavor="seurat",
         n_top_genes=args.n_top_genes,
+        flavor="seurat",
         inplace=True,
+        batch_key=args.batch_key,
+    )
+    full_features = adata.var.copy()
+    hvg_full = (
+        adata.var[adata.var["highly_variable"]]
+        .sort_values(
+            by=["highly_variable_nbatches", "dispersions_norm"],
+            ascending=[False, False]
+        )
+        .head(args.n_top_genes)
+        .index.tolist()
     )
+    hvg_adata = adata[:, adata.var.index.isin(hvg_full)]
+    hvg_features = hvg_adata.var.copy()
+
     sc.pl.highly_variable_genes(
-        adata,
+        hvg_adata,
     )
     fig = plt.gcf()
-    fig.suptitle(f"Highly variable genes dispersion plot - {args.plots_prefix}", va="center", ha="center", fontsize=16)
-    plt.savefig(f"{args.plots_prefix}.hvg_dispersion.png", dpi=300, bbox_inches="tight")
+    fig.suptitle(f"Highly variable genes dispersion plot - {args.output_prefix}", va="center", ha="center", fontsize=16)
+    plt.savefig(f"{args.output_prefix}.hvg_dispersion.png", dpi=300, bbox_inches="tight")
 
 
     ##############################
     ## DIMENSIONALITY REDUCTION ##
     ##############################
     sc.pp.pca(
-        adata,
+        hvg_adata,
         n_comps=args.n_comps,
         svd_solver="arpack",
     )
 
     # Save outputs
-    adata.write_h5ad(filename=args.adata_output, compression="gzip")
+    full_features.to_csv(f"{args.output_prefix}.all_genes.csv", index=True)
+    hvg_features.to_csv(f"{args.output_prefix}.hvg_genes.csv", index=True)
+    hvg_adata.write_h5ad(filename=args.adata_output, compression="gzip")
 
 
 if __name__ == "__main__":
@@ -128,12 +144,19 @@ if __name__ == "__main__":
         required=True,
         help="Number of principal components to compute using `scanpy.pp.pca` [30]"
     )
+    parser.add_argument(
+        "-b",
+        "--batch-key",
+        type=str,
+        required=True,
+        help="Key in AnnData object for batch information ['batch_id']"
+    )
     parser.add_argument(
         "-p",
-        "--plots-prefix",
+        "--output-prefix",
         type=str,
         required=True,
-        help="Output file name prefix for the HVGs plot"
+        help="Output file name prefix for the HVGs plot and CSVs"
     )
     parser.add_argument(
         "-o",

docker/spatial_py/scripts/visium_spatially_variable_genes

@@ -43,6 +43,8 @@ def main(args):
 
     # Save table and adata object
     top_10_variable_genes.to_csv(f"{args.cohort_id}.moran_top_10_variable_genes.csv")
+    metadata = adata.obs
+    metadata.to_csv(f"{args.cohort_id}.final_metadata.csv")
     adata.write_h5ad(filename=args.adata_output, compression="gzip")