Merge pull request #389 from broadinstitute/development

Release 1.41.0

Author: bistline

ingest/anndata_.py

@@ -36,6 +36,7 @@ def __init__(self, file_path, study_file_id, study_id, **kwargs):
         IngestFiles.__init__(
             self, file_path, allowed_file_types=self.ALLOWED_FILE_TYPES
         )
+        self.kwargs = kwargs
 
     def obtain_adata(self):
         try:
@@ -58,6 +59,36 @@ def basic_validation(self):
         except ValueError:
             return False
 
+    def validate_raw_location(self):
+        """
+        Confirm file has data at raw_location
+        """
+        adata = self.obtain_adata()
+        raw_location = self.kwargs.get("raw_location")
+        if raw_location is not None:
+            try:
+                if raw_location == ".raw":
+                    if adata.raw is None:
+                        msg = f'No data found in .raw slot'
+                        log_exception(
+                            IngestFiles.dev_logger, IngestFiles.user_logger, msg
+                        )
+                        raise ValueError(msg)
+                else:
+                    if raw_location not in adata.layers.keys():
+                        msg = f'No data found at adata.layers["{raw_location}"]'
+                        log_exception(
+                            IngestFiles.dev_logger, IngestFiles.user_logger, msg
+                        )
+                        raise ValueError(msg)
+                return True
+            except ValueError:
+                return False
+        else:
+            msg = 'Must specify location of raw counts in AnnData object'
+            log_exception(IngestFiles.dev_logger, IngestFiles.user_logger, msg)
+            return False
+
     def create_cell_data_arrays(self):
         """Extract cell name DataArray documents for raw data"""
         adata = self.obtain_adata()

ingest/cli_parser.py

@@ -1,5 +1,4 @@
-"""Helper functions for ingest_pipeline.py
-"""
+"""Helper functions for ingest_pipeline.py"""
 
 import argparse
 import ast
@@ -451,6 +450,12 @@ def create_parser():
         help="Array of obsm key(s) to extract as cluster files",
     )
 
+    parser_anndata.add_argument(
+        "--raw-location",
+        help="location of raw counts. '.raw' for raw slot, "
+        "else adata.layers key value or None if no raw counts",
+    )
+
     parser_anndata.add_argument(
         "--extract",
         type=ast.literal_eval,

ingest/de.py

@@ -367,7 +367,7 @@ def write_de_result(adata, group, annotation, rank_key, cluster_name, extra_para
             clean_group = DifferentialExpression.sanitize_string(group)
             out_file = f'{cluster_name}--{clean_annotation}--{clean_group}--{annot_scope}--{method}.tsv'
             DifferentialExpression.de_logger.info(
-                f"Writing DE output for {clean_group} vs rest"
+                f"Writing DE output for {clean_group} vs restq"
             )
         elif de_type == "pairwise":
             # rank_genes_groups accepts a list. For SCP pairwise, should be a list with one item
@@ -403,6 +403,7 @@ def run_scanpy_de(
     ):
         method = extra_params.get("method")
         de_type = extra_params.get("de_type")
+        raw_location = extra_params.get("raw_location")
 
         try:
             DifferentialExpression.assess_annotation(annotation, metadata, extra_params)
@@ -432,24 +433,50 @@ def run_scanpy_de(
             )
 
         if matrix_file_type == "h5ad":
-            if orig_adata.raw is not None:
-                adata = AnnData(
-                    # using .copy() for the AnnData components is good practice
-                    # but we won't be using orig_adata for analyses
-                    # choosing to avoid .copy() for memory efficiency
-                    X=orig_adata.raw.X,
-                    obs=orig_adata.obs,
-                    var=orig_adata.var,
-                )
+            if raw_location == ".raw":
+                if orig_adata.raw is not None:
+                    DifferentialExpression.de_logger.info(
+                        f"Performing DE on {raw_location} data"
+                    )
+                    adata = AnnData(
+                        # using .copy() for the AnnData components is good practice
+                        # but we won't be using orig_adata for analyses
+                        # choosing to avoid .copy() for memory efficiency
+                        X=orig_adata.raw.X,
+                        obs=orig_adata.obs,
+                        var=orig_adata.var,
+                    )
+                else:
+                    msg = f'{matrix_file_path} does not have a .raw attribute'
+                    print(msg)
+                    log_exception(
+                        DifferentialExpression.dev_logger,
+                        DifferentialExpression.de_logger,
+                        msg,
+                    )
+                    raise ValueError(msg)
             else:
-                msg = f'{matrix_file_path} does not have a .raw attribute'
-                print(msg)
-                log_exception(
-                    DifferentialExpression.dev_logger,
-                    DifferentialExpression.de_logger,
-                    msg,
-                )
-                raise ValueError(msg)
+                if raw_location in orig_adata.layers.keys():
+                    DifferentialExpression.de_logger.info(
+                        f"Performing DE on adata.layers['{raw_location}'] data"
+                    )
+                    adata = AnnData(
+                        # using .copy() for the AnnData components is good practice
+                        # but we won't be using orig_adata for analyses
+                        # choosing to avoid .copy() for memory efficiency
+                        X=orig_adata.layers[raw_location],
+                        obs=orig_adata.obs,
+                        var=orig_adata.var,
+                    )
+                else:
+                    msg = f'{matrix_file_path} does not have adata.layers["{raw_location}"]'
+                    print(msg)
+                    log_exception(
+                        DifferentialExpression.dev_logger,
+                        DifferentialExpression.de_logger,
+                        msg,
+                    )
+                    raise ValueError(msg)
         # AnnData expects gene x cell so dense and mtx matrices require transposition
         else:
             adata = adata.transpose()

ingest/ingest_pipeline.py

@@ -55,7 +55,7 @@
 python ingest_pipeline.py  --study-id 5d276a50421aa9117c982845 --study-file-id 5dd5ae25421aa910a723a337 ingest_anndata --ingest-anndata --anndata-file ../tests/data/anndata/trimmed_compliant_pbmc3K.h5ad  --extract "['processed_expression']"
 
 # Ingest AnnData - happy path raw count cell name only extraction
-python ingest_pipeline.py  --study-id 5d276a50421aa9117c982845 --study-file-id 5dd5ae25421aa910a723a337 ingest_anndata --ingest-anndata --anndata-file ../tests/data/anndata/trimmed_compliant_pbmc3K.h5ad  --extract "['raw_counts']"
+python ingest_pipeline.py  --study-id 5d276a50421aa9117c982845 --study-file-id 5dd5ae25421aa910a723a337 ingest_anndata --ingest-anndata --anndata-file ../tests/data/anndata/trimmed_compliant_pbmc3K.h5ad  --extract "['raw_counts']" --raw-location ".raw"
 
 # Ingest AnnData - happy path cluster and metadata extraction
 python ingest_pipeline.py  --study-id 5d276a50421aa9117c982845 --study-file-id 5dd5ae25421aa910a723a337 ingest_anndata --ingest-anndata --anndata-file ../tests/data/anndata/trimmed_compliant_pbmc3K.h5ad  --extract "['cluster', 'metadata']" --obsm-keys "['X_umap','X_tsne']"
@@ -66,17 +66,20 @@
 # Differential expression analysis (sparse matrix)
 python ingest_pipeline.py --study-id addedfeed000000000000000 --study-file-id dec0dedfeed1111111111111 differential_expression --annotation-name cell_type__ontology_label --annotation-type group --annotation-scope study --matrix-file-path ../tests/data/differential_expression/sparse/sparsemini_matrix.mtx --gene-file ../tests/data/differential_expression/sparse/sparsemini_features.tsv --barcode-file ../tests/data/differential_expression/sparse/sparsemini_barcodes.tsv --matrix-file-type mtx --annotation-file ../tests/data/differential_expression/sparse/sparsemini_metadata.txt --cluster-file ../tests/data/differential_expression/sparse/sparsemini_cluster.txt --cluster-name de_sparse_integration --study-accession SCPsparsemini --differential-expression
 
-# Differential expression analysis (h5ad matrix)
-python ingest_pipeline.py --study-id addedfeed000000000000000 --study-file-id dec0dedfeed1111111111111 differential_expression --annotation-name louvain --annotation-type group --annotation-scope study --matrix-file-path ../tests/data/anndata/trimmed_compliant_pbmc3K.h5ad --matrix-file-type h5ad --annotation-file ../tests/data/anndata/h5ad_frag.metadata.tsv --cluster-file ../tests/data/anndata/h5ad_frag.cluster.X_umap.tsv --cluster-name umap --study-accession SCPdev --differential-expression
+# Differential expression analysis (h5ad matrix, raw count in raw slot)
+python ingest_pipeline.py --study-id addedfeed000000000000000 --study-file-id dec0dedfeed1111111111111 differential_expression  --raw-location '.raw' --annotation-name cell_type__ontology_label --de-type rest  --annotation-type group --annotation-scope study --annotation-file ../tests/data/anndata/compliant_liver_h5ad_frag.metadata.tsv.gz --cluster-file ../tests/data/anndata/compliant_liver_h5ad_frag.cluster.X_umap.tsv.gz --cluster-name umap --matrix-file-path ../tests/data/anndata/compliant_liver.h5ad  --matrix-file-type h5ad --study-accession SCPdev --differential-expression
+
+# Differential expression analysis (h5ad matrix, raw count in adata.layers['counts'])
+python ingest_pipeline.py --study-id addedfeed000000000000000 --study-file-id dec0dedfeed1111111111111 differential_expression  --raw-location 'counts' --annotation-name cell_type__ontology_label --de-type rest  --annotation-type group --annotation-scope study --annotation-file ../tests/data/anndata/compliant_liver_h5ad_frag.metadata.tsv.gz --cluster-file ../tests/data/anndata/compliant_liver_h5ad_frag.cluster.X_umap.tsv.gz --cluster-name umap --matrix-file-path ../tests/data/anndata/compliant_liver_layers_counts.h5ad  --matrix-file-type h5ad --study-accession SCPdev --differential-expression
 
 # Pairwise differential expression analysis (dense matrix)
 python ingest_pipeline.py --study-id addedfeed000000000000000 --study-file-id dec0dedfeed1111111111111 differential_expression --annotation-name cell_type__ontology_label --de-type pairwise --group1 "['cholinergic neuron']" --group2 "cranial somatomotor neuron" --annotation-type group --annotation-scope study --matrix-file-path ../tests/data/differential_expression/de_dense_matrix.tsv --matrix-file-type dense --annotation-file ../tests/data/differential_expression/de_dense_metadata.tsv --cluster-file ../tests/data/differential_expression/de_dense_cluster.tsv --cluster-name de_integration --study-accession SCPdev --differential-expression
 
 # Pairwise differential expression analysis (sparse matrix)
 python ingest_pipeline.py --study-id addedfeed000000000000000 --study-file-id dec0dedfeed1111111111111 differential_expression --annotation-name cell_type__ontology_label --de-type pairwise --group1 "['endothelial cell']" --group2 "smooth muscle cell" --annotation-type group --annotation-scope study --matrix-file-path ../tests/data/differential_expression/sparse/sparsemini_matrix.mtx --gene-file ../tests/data/differential_expression/sparse/sparsemini_features.tsv --barcode-file ../tests/data/differential_expression/sparse/sparsemini_barcodes.tsv --matrix-file-type mtx --annotation-file ../tests/data/differential_expression/sparse/sparsemini_metadata.txt --cluster-file ../tests/data/differential_expression/sparse/sparsemini_cluster.txt --cluster-name de_sparse_integration --study-accession SCPsparsemini --differential-expression
 
-# Pairwise differential expression analysis (h5ad matrix)
-python ingest_pipeline.py --study-id addedfeed000000000000000 --study-file-id dec0dedfeed1111111111111 differential_expression --annotation-name cell_type__ontology_label --de-type pairwise --group1 "['mature B cell']" --group2 "plasma cell" --annotation-type group --annotation-scope study --annotation-file ../tests/data/anndata/compliant_liver_h5ad_frag.metadata.tsv.gz --cluster-file ../tests/data/anndata/compliant_liver_h5ad_frag.cluster.X_umap.tsv.gz --cluster-name umap --matrix-file-path ../tests/data/anndata/compliant_liver.h5ad  --matrix-file-type h5ad --study-accession SCPdev --differential-expression
+# Pairwise differential expression analysis (h5ad matrix, raw count in raw slot)
+python ingest_pipeline.py --study-id addedfeed000000000000000 --study-file-id dec0dedfeed1111111111111 differential_expression  --raw-location '.raw' --annotation-name cell_type__ontology_label --de-type pairwise --group1 "mature B cell" --group2 "plasma cell" --annotation-type group --annotation-scope study --annotation-file ../tests/data/anndata/compliant_liver_h5ad_frag.metadata.tsv.gz --cluster-file ../tests/data/anndata/compliant_liver_h5ad_frag.cluster.X_umap.tsv.gz --cluster-name umap --matrix-file-path ../tests/data/anndata/compliant_liver.h5ad  --matrix-file-type h5ad --study-accession SCPdev --differential-expression
 
 """
 
@@ -559,7 +562,11 @@ def extract_from_anndata(self):
             if self.kwargs.get('extract') and "raw_counts" in self.kwargs.get(
                 'extract'
             ):
-                self.anndata.ingest_raw_cells()
+                if self.anndata.validate_raw_location():
+                    self.anndata.ingest_raw_cells()
+                else:
+                    self.report_validation("failure")
+                    return 1
             self.report_validation("success")
             return 0
         # scanpy unable to open AnnData file

ingest/validation/ontologies/cl.min.tsv.gz

@@ -1,2 +1 @@
-1738072997 # validation cache key
-
+1742404288 # validation cache key