diff --git a/docker/docker-compose.yml.template b/docker/docker-compose.yml.template index ae2f63fe5..c170dcaa8 100644 --- a/docker/docker-compose.yml.template +++ b/docker/docker-compose.yml.template @@ -86,6 +86,23 @@ services: - ../lib:/gEAR/lib - ../www/uploads/files:/gEAR/www/uploads/files + # AnnData object Upload RabbitMQ Consumer + # i.e. h5ad, mex, tar + anndata_upload_consumer: + build: + context: .. + dockerfile: listeners/Dockerfile.anndata_upload + depends_on: + - queue + image: adkinsrs/gear_anndata_upload_consumer:latest + networks: + - gear + restart: on-failure + volumes: + - ../listeners:/gEAR/listeners + - ../lib:/gEAR/lib + - ../www/uploads/files:/gEAR/www/uploads/files + volumes: db: driver: local diff --git a/docs/developer/services/rabbitmq_consumers.md b/docs/developer/services/rabbitmq_consumers.md index 960ac056b..8d00809dd 100644 --- a/docs/developer/services/rabbitmq_consumers.md +++ b/docs/developer/services/rabbitmq_consumers.md @@ -25,6 +25,16 @@ Result Storage (Database, File System) ## Available Consumers +### Anndata Upload Consume + +Facilitates uploading of various file formats into Anndata (H5AD) format. + +- **Listener**: `listeners/anndata_upload_consumer.py` +- **Queue**: `anndata_upload_jobs` +- **Service Template**: `systemd/anndata-upload-consumer@.service` +- **Service Group**: `systemd/anndata-upload-consumer.target` + + ### Gosling Upload Consume Facilitates uploading of track files for the epigenome uploader. diff --git a/lib/gear/anndata_processor.py b/lib/gear/anndata_processor.py new file mode 100644 index 000000000..6f5dde0eb --- /dev/null +++ b/lib/gear/anndata_processor.py @@ -0,0 +1,581 @@ +""" +AnndataProcessor - Process expression datasets in various formats to H5AD. + +Handles conversion, sanitization, and primary analysis of expression datasets +in formats: H5AD, 3-tab, Excel, MEX. +""" + +import gc +import json +import os +from pathlib import Path +import tarfile +import zipfile + +import anndata +import pandas as pd +from scipy import sparse + +import geardb +from gear.primary_analysis import add_primary_analysis_to_dataset, PrimaryAnalysisProcessingError +from gear.utils import update_var_with_ensembl_ids + + + +def write_status(status_file, status): + """Write status dictionary to JSON file.""" + with open(status_file, 'w') as f: + f.write(json.dumps(status, indent=2)) + +def process_anndata_synchronously(job_id, share_uid, staging_area, status_file, dataset_uid, dataset_format, perform_primary_analysis): + """Process anndata upload synchronously.""" + processor = AnndataProcessor( + job_id=job_id, + share_uid=share_uid, + staging_area=staging_area, + status_file=status_file, + dataset_uid=dataset_uid, + ) + return processor.process(dataset_format, perform_primary_analysis) + +def clean_chunk(chunk: pd.DataFrame) -> pd.DataFrame: + """ + Clean a dataframe chunk by removing whitespace and converting to numeric. + + Args: + chunk: DataFrame chunk to clean + + Returns: + Cleaned DataFrame + """ + chunk = chunk.replace(r'^\s+|\s+$', '', regex=True) + chunk = chunk.apply( + lambda col: col.map(lambda x: x.strip() if isinstance(x, str) else x) + ) + chunk = chunk.apply(pd.to_numeric, errors='coerce').fillna(0) + return chunk + +def sanitize_obs_for_h5ad(obs_df: pd.DataFrame) -> pd.DataFrame: + """Sanitize observation dataframe for H5AD storage.""" + for col in obs_df.columns: + if obs_df[col].dtype == 'object': + obs_df[col] = obs_df[col].fillna('').astype(str) + return obs_df + +def categorize_observation_columns(obs: pd.DataFrame) -> None: + """Categorize and convert specific observation columns.""" + for str_type in ['cell_type', 'condition', 'time_point', 'time_unit']: + if str_type in obs.columns: + obs[str_type] = pd.Categorical(obs[str_type]) + + for num_type in ['replicate', 'time_point_order']: + if num_type in obs.columns: + obs[num_type] = pd.to_numeric(obs[num_type]) + + +def package_content_type(filenames: list[str]) -> str | None: + #print("DEBUG: filenames", file=sys.stderr, flush=True) + #print(filenames, file=sys.stderr, flush=True) + """ + mex: + matrix.mtx + barcodes.tsv + genes.tsv + + threetab: + expression.tab + genes.tab + observations.tab + + None is returned if neither of these is true + + Added NEMO file format functionality. + DataMTX.tab -> expression.tab + COLmeta.tab -> observations.tab + ROWmeta.tab -> genes.tab + """ + if 'expression.tab' in filenames and 'genes.tab' in filenames and 'observations.tab' in filenames: + return 'threetab' + + if 'matrix.mtx' in filenames and 'barcodes.tsv' in filenames and 'genes.tsv' in filenames: + return 'mex' + + if 'DataMTX.tab' in filenames and 'COLmeta.tab' in filenames and 'ROWmeta.tab' in filenames: + return 'threetab' + + return None + +class ProcessingError(Exception): + """Raised when dataset processing fails.""" + pass + +class AnndataProcessor: + """Process and convert expression datasets to H5AD format.""" + + def __init__( + self, + job_id: str, + share_uid: str, + staging_area: Path, + status_file: Path, + dataset_uid: str, + ) -> None: + """ + Initialize the processor. + + Args: + job_id: Unique job identifier + share_uid: Share UID for the dataset + staging_area: Directory containing uploaded files + status_file: Path to status.json for progress updates + dataset_uid: Dataset UID for primary analysis + """ + self.job_id = job_id + self.share_uid = share_uid + self.staging_area = staging_area + self.status_file = status_file + self.dataset_uid = dataset_uid + self.status = { + "job_id": self.job_id, + "status": "processing", + "message": "", + "progress": 0, + } + + def process( + self, + dataset_format: str, + perform_primary_analysis: bool = False, + ) -> dict: + """ + Process the uploaded dataset. + + Args: + dataset_format: Format of the dataset (h5ad, threetab, excel, mex) + perform_primary_analysis: Whether to perform primary analysis after processing + + Returns: + Result dictionary with 'success' and 'message' keys + """ + try: + h5ad_path = self._process_by_format(dataset_format) + + if perform_primary_analysis: + self._update_progress(66, "Performing primary analysis...") + try: + add_primary_analysis_to_dataset( + self.dataset_uid, + self.share_uid, + self.staging_area, + dataset_format, + ) + except PrimaryAnalysisProcessingError as e: + raise ProcessingError(f"Primary analysis failed: {str(e)}") + + self._update_progress(100, "Dataset processed successfully.") + return {"success": 1, "message": "Dataset processed successfully."} + + except ProcessingError as e: + self._update_status("error", str(e)) + return {"success": 0, "message": str(e)} + except Exception as e: + import traceback + traceback.print_exc() + self._update_status("error", f"Unexpected error: {str(e)}") + return {"success": 0, "message": f"Unexpected error: {str(e)}"} + finally: + gc.collect() + + def _process_by_format(self, dataset_format: str) -> Path: + """Route to appropriate processor based on dataset format.""" + if dataset_format == "h5ad": + return self._process_h5ad() + elif dataset_format == "mex_3tab": + return self._process_mex_3tab() + elif dataset_format == "excel": + return self._process_excel() + elif dataset_format == "rdata": + raise NotImplementedError("RData format processing not yet implemented.") + return self._process_rdata() + else: + raise ProcessingError(f"Unsupported dataset format: {dataset_format}") + + def _process_rdata(self): + pass + + def _process_h5ad(self) -> Path: + """Process .h5ad file with backed mode for memory efficiency.""" + self._update_progress(5, "Reading H5AD file...") + + filepath = self.staging_area / f"{self.share_uid}.h5ad" + + # Read in backed mode to avoid loading full file into memory + adata = anndata.read_h5ad(filepath, backed='r') + + self._update_progress(15, "Sanitizing observation metadata...") + + # Convert obs to memory (typically small) + obs = adata.obs + categorize_observation_columns(obs) + obs = sanitize_obs_for_h5ad(obs) + + # Update var if gene_symbol missing + var = adata.var if "gene_symbol" not in adata.var.columns else None + if var is not None: + self._update_progress(25, "Mapping gene symbols via Ensembl...") + var = self._update_var_with_ensembl_ids(var) + + # Close backed file before final write + adata.file.close() + + self._update_progress(50, "Writing H5AD file...") + + # Re-read for full write (necessary for consistency) + adata = anndata.read_h5ad(filepath) + adata.obs = obs + if var is not None: + adata.var = var + + h5ad_temp = self.staging_area / f"{self.share_uid}.new.h5ad" + adata.write(h5ad_temp, compression='gzip') + + filepath.unlink() + h5ad_temp.rename(filepath) + + self._update_progress(65, "H5AD processing complete.") + return filepath + + def _process_mex_3tab(self) -> Path: + # Extract the file + compression_format = None + filename = self.staging_area / f"{self.share_uid}.tar.gz" + + if filename.exists(): + compression_format = 'tarball' + else: + filename = self.staging_area / f"{self.share_uid}.zip" + + if filename.exists(): + compression_format = 'zip' + else: + raise ProcessingError("No tarball or zip file found for MEX/3-tab dataset.") + + files_extracted = [] + + if compression_format == 'tarball': + try: + with tarfile.open(filename) as tf: + for entry in tf: + tf.extract(entry, path=self.staging_area) + + # Nemo suffixes + nemo_suffixes = ['DataMTX.tab', 'COLmeta.tab', 'ROWmeta.tab'] + suffix_found = None + + for suffix in nemo_suffixes: + if entry.name.endswith(suffix): + suffix_found = suffix + # Rename the file to the appropriate name + old_name = self.staging_area / entry.name + new_name = self.staging_area / suffix + old_name.replace(new_name) + + if suffix_found is not None: + files_extracted.append(suffix_found) + else: + files_extracted.append(entry.name) + except tarfile.ReadError: + raise ProcessingError("Bad tarball file. Couldn't extract the tarball.") + + if compression_format == 'zip': + try: + with zipfile.ZipFile(filename) as zf: + for entry in zf.infolist(): + zf.extract(entry, path=self.staging_area) + + # Nemo suffixes + nemo_suffixes = ['DataMTX.tab', 'COLmeta.tab', 'ROWmeta.tab'] + suffix_found = None + + for suffix in nemo_suffixes: + if entry.filename.endswith(suffix): + suffix_found = suffix + # Rename the file to the appropriate name + old_name = self.staging_area / entry.filename + new_name = self.staging_area / suffix + old_name.replace(new_name) + + if suffix_found is not None: + files_extracted.append(suffix_found) + else: + files_extracted.append(entry.filename) + except zipfile.BadZipFile: + raise ProcessingError("Bad zip file. Couldn't extract the zip file.") + + # Determine the dataset type + dataset_type = package_content_type(files_extracted) + + if dataset_type is None: + raise ProcessingError("Unsupported dataset format. Couldn't tell type from file names within the archive.") + + # Call the appropriate function + if dataset_type == 'threetab': + return self._process_threetab() + elif dataset_type == 'mex': + return self._process_mex() + + # If no valid dataset type is found, raise an error + raise ProcessingError("Failed to process MEX/3-tab dataset. No valid dataset type identified.") + + + def _process_threetab(self) -> Path: + """Process 3-tab format (expression, genes, observations) with chunking.""" + self._update_progress(5, "Extracting and validating files...") + + # Extract/find the three required files + expression_matrix_path, obs, var = self._extract_threetab_files() + + self._update_progress(15, "Categorizing observations...") + categorize_observation_columns(obs) + + self._update_progress(25, "Processing expression matrix in chunks...") + + # Process expression matrix in chunks + chunk_size = 500 + expression_matrix = self._read_expression_matrix_chunks( + expression_matrix_path, chunk_size, total_rows=len(var) + ) + + # Create AnnData object + adata = anndata.AnnData(X=expression_matrix, obs=var, var=obs) + adata = adata.transpose() + adata.obs = sanitize_obs_for_h5ad(adata.obs) # type: ignore + + self._update_progress(50, "Writing H5AD file...") + + h5ad_path = self.staging_area / f"{self.share_uid}.h5ad" + adata.write(h5ad_path, compression='gzip') + + self._update_progress(65, "3-tab processing complete.") + return h5ad_path + + def _extract_threetab_files(self) -> tuple[Path, pd.DataFrame, pd.DataFrame]: + """Extract and read the three required 3-tab files.""" + expression_matrix_path = None + obs = None + var = None + + for infile in self.staging_area.iterdir(): + if infile.name.startswith('.'): + continue + + filepath = infile + + if infile.name in ['expression.tab', 'DataMTX.tab']: + expression_matrix_path = filepath + elif infile.name in ['observations.tab', 'COLmeta.tab']: + obs = pd.read_table(filepath, sep='\t', index_col=0, header=0) + elif infile.name in ['genes.tab', 'ROWmeta.tab']: + var = pd.read_table(filepath, sep='\t', index_col=0, header=0) + + if obs is None: + raise ProcessingError("No observations file found (expected observations.tab or COLmeta.tab).") + if var is None: + raise ProcessingError("No genes file found (expected genes.tab or ROWmeta.tab).") + if expression_matrix_path is None: + raise ProcessingError("No expression file found (expected expression.tab or DataMTX.tab).") + + return expression_matrix_path, obs, var + + def _process_excel(self) -> Path: + """Process Excel file with expression, observations, and genes sheets.""" + self._update_progress(5, "Reading Excel file...") + + filepath = self.staging_area / f"{self.share_uid}.xlsx" + exp_df = pd.read_excel(filepath, sheet_name='expression', index_col=0).transpose() + + try: + X = exp_df.to_numpy()[:, 0:].astype(float) + except ValueError: + raise ProcessingError("Encountered unexpected value type. Expected float type in expression matrix.") + + number_obs_from_exp, number_genes_from_exp = X.shape + + self._update_progress(20, "Reading observations sheet...") + + try: + obs_df = pd.read_excel(filepath, sheet_name='observations', index_col='observations') + except ValueError: + raise ProcessingError("No observations sheet found. Expected spreadsheet sheet named 'observations'.") + + # Validate observations + number_obs = len(obs_df) + if number_obs != number_obs_from_exp: + raise ProcessingError( + f"Observations sheet error. Row count ({number_obs}) must match " + f"expression sheet row count ({number_obs_from_exp})." + ) + + if not obs_df.index.equals(exp_df.index): + raise ProcessingError( + "Observations sheet error. Index names and order must match expression sheet rows." + ) + + self._update_progress(35, "Reading genes sheet...") + + genes_df = self._extract_genes_sheet(filepath, exp_df) + + # Validate gene_symbol column + if 'gene_symbol' not in genes_df.columns: + raise ProcessingError("Failed to find gene_symbol column in genes tab.") + + digit_count = genes_df['gene_symbol'].str.isnumeric().sum() + if digit_count > 0: + raise ProcessingError(f"{digit_count} gene symbols are numbers, not gene symbols.") + + categorize_observation_columns(obs_df) + + # Validate gene count + number_genes = len(genes_df) + if number_genes != number_genes_from_exp: + raise ProcessingError( + f"Genes sheet error. Row count ({number_genes}) must match " + f"expression sheet row count ({number_genes_from_exp})." + ) + + if not genes_df.index.equals(exp_df.columns): + raise ProcessingError("Genes sheet error. Index names and order must match expression sheet columns.") + + self._update_progress(50, "Creating AnnData object...") + + adata = anndata.AnnData(X=X, obs=obs_df, var=genes_df) + adata.obs = sanitize_obs_for_h5ad(adata.obs) # type: ignore + + self._update_progress(60, "Writing H5AD file...") + + h5ad_path = self.staging_area / f"{self.share_uid}.h5ad" + adata.write(h5ad_path, compression='gzip') + + self._update_progress(65, "Excel processing complete.") + return h5ad_path + + def _extract_genes_sheet(self, filepath: Path, exp_df: pd.DataFrame) -> pd.DataFrame: + """Extract genes sheet, falling back to expression sheet if needed.""" + try: + genes_df = pd.read_excel( + filepath, sheet_name='genes', index_col=0, + converters={'gene_symbol': str} + ) + except ValueError: + try: + # Try getting genes from expression sheet + genes_df = pd.read_excel( + filepath, sheet_name='expression', index_col=0, + usecols=[0, 1] + ) + genes_df = genes_df.drop(genes_df.columns[0], axis=1) + except Exception as err: + raise ProcessingError(f"No 'genes' sheet found. {str(err)}") + + return genes_df + + def _process_mex(self) -> Path: + """Process MEX format (matrix.mtx, barcodes.tsv, genes.tsv).""" + raise NotImplementedError("MEX format processing not yet implemented.") + + def _read_expression_matrix_chunks( + self, filepath: Path, chunk_size: int, total_rows: int + ) -> sparse.csr_matrix: + """Read expression matrix in chunks with progress updates.""" + total_chunks = (total_rows + chunk_size - 1) // chunk_size + + # Try reading without cleanup first + expression_matrix = self._read_chunks_with_cleanup( + filepath, chunk_size, total_rows, total_chunks, cleanup=False + ) + + if expression_matrix: + return sparse.csr_matrix(sparse.vstack(expression_matrix)) + + # Retry with data cleanup if first attempt fails + expression_matrix = self._read_chunks_with_cleanup( + filepath, chunk_size, total_rows, total_chunks, cleanup=True + ) + if expression_matrix: + return sparse.csr_matrix(sparse.vstack(expression_matrix)) + + raise ProcessingError("Failed to read expression matrix in chunks, even after cleanup.") + + def _read_chunks_with_cleanup( + self, + filepath: Path, + chunk_size: int, + total_rows: int, + total_chunks: int, + cleanup: bool = False, + ) -> list: + """ + Read expression matrix in chunks, optionally with data cleanup. + + Args: + filepath: Path to expression matrix file + chunk_size: Number of rows per chunk + total_rows: Total number of rows + total_chunks: Total number of chunks + cleanup: Whether to apply data cleanup (whitespace, type conversion) + + Returns: + List of sparse matrices, empty list if reading fails + """ + expression_matrix = [] + rows_read = 0 + reader = pd.read_csv(filepath, sep='\t', index_col=0, chunksize=chunk_size) + + try: + for chunk_idx, chunk in enumerate(reader, 1): + if cleanup: + chunk = clean_chunk(chunk) + + rows_read += len(chunk) + pct = int((rows_read / total_rows) * 100) + self._update_progress( + 25 + round(pct * 0.4), + f"Processing chunk {chunk_idx}/{total_chunks}" + ) + expression_matrix.append(sparse.csr_matrix(chunk.values)) + + return expression_matrix + + except Exception: + return [] + + def _update_var_with_ensembl_ids(self, var_df: pd.DataFrame) -> pd.DataFrame: + """Update var dataframe with Ensembl IDs.""" + metadata_file = self.staging_area / 'metadata.json' + with open(metadata_file, 'r') as f: + metadata = json.load(f) + + sample_taxid = metadata.get("sample_taxid", None) + organism_id = geardb.get_organism_id_by_taxon_id(sample_taxid) + + if not organism_id: + raise ProcessingError("Could not determine organism ID from sample taxonomic ID.") + + return update_var_with_ensembl_ids(var_df, organism_id, "UNMAPPED_") + + def _update_progress(self, progress: int, message: str) -> None: + """Update progress and write status file.""" + progress = max(0, min(100, progress)) # Clamp to 0-100 + self.status['progress'] = progress + self.status['message'] = message + self._write_status_file() + + def _update_status(self, status_name: str, message: str) -> None: + """Update status and write status file.""" + self.status['status'] = status_name + self.status['message'] = message + self._write_status_file() + + def _write_status_file(self) -> None: + """Write current status to status.json.""" + write_status(self.status_file, self.status) \ No newline at end of file diff --git a/lib/gear/dataarchive.py b/lib/gear/dataarchive.py deleted file mode 100644 index b4f73bb8b..000000000 --- a/lib/gear/dataarchive.py +++ /dev/null @@ -1,178 +0,0 @@ -import os,sys -import anndata -import pandas as pd -import scanpy as sc - - -class DataArchive: - """ - This class handles: - 1. Get dataset type (mex or 3tab) - 2. Create h5ad and return - 3. Write h5ad to specfic directory - """ - - def __init__(self, base_dir=None, extracted_file_paths=None, format=None, tarball_path=None): - # If the archive is a tarball, this is the source path for it - self.tarball_path = tarball_path - - # If the archive is extracted into a directory, this should be - # that directory base, directly within which all files are found. - self.base_dir = base_dir - - # The full path to each file in the archive - self.extracted_file_paths = list() - - # Can be explicit or derived (mex, 3tab, etc.) - self.format = format - - - def get_archive_type(self, data_path = None): - """ - Determines type of dataset based on file names and extensions. - Input - ----- - data_path - path/to/extracted/data - Output - ----- - return "mex" or "3tab" - """ - archive_filenames = os.listdir(data_path) - archive_files=''.join(archive_filenames) - data_type = None - if 'expression.tab' in archive_files and 'genes.tab' in archive_files and 'observations.tab' in archive_files: - data_type = "3tab" - if 'matrix.mtx' in archive_files and 'barcodes.tsv' in archive_files and 'genes.tsv' in archive_files: - data_type = "mex" - if 'DataMTX.tab' in archive_files and 'COLmeta.tab' in archive_files and 'ROWmeta.tab' in archive_files: - data_type = "3tab" - return data_type - - def read_mex_files(self, data_path = None): - """ - Reads all files in the direcory and returns an Anndata object. - - Input - ----- - data_path - path/to/extracted/archive - - Output - ------ - Output is anndata object created using the files in directory - 'adata' is an AnnData object where data of each MEX file is assigned to the AnnData object: - AnnData.X = matrix.mtx - AnnData.obs = barcodes.tsv - AnnData.var = genes.tsv - - """ - is_en=False - archive_files_list = os.listdir(data_path) - archive_files = [os.path.join(data_path, s) for s in archive_files_list] - for entry in archive_files: - if 'matrix.mtx' in entry: - adata = sc.read(entry, cache=False).transpose() - elif 'barcodes.tsv' in entry: - obs = pd.read_csv(entry, sep='\t', header=None,index_col = 0, names=['observations']) - elif 'genes.tsv' in entry: - var = pd.read_csv(entry, sep='\t', header=None, index_col = 0, names=['genes', 'gene_symbol']) - if var.index.str.contains('ENS').sum() >1: - is_en=True - elif 'EXPmeta.json' in entry: - continue - else: - raise Exception("Unexpected file name: '{0}'. Excepted 'matrix.mtx', 'genes.tsv', or 'barcodes.tsv'.".format(entry)) - adata.var = var - adata.obs = obs - self.adata = adata - self.originalPath = data_path - return self, is_en - - def read_3tab_files(self, data_path= None): - """ - Reads files in direcoty and returns an Anndata object - - Input - ----- - data_path - /path/to/extracted/archive - - Output - ------ - Output is anndata object created using the files in directory - 'adata' is an AnnData object where data of each TAB file is assigned to the AnnData object: - AnnData.X = expression file - AnnData.obs = observations file - AnnData.var = genes file - """ - is_en = False - archive_files_list = os.listdir(data_path) - archive_files = [os.path.join(data_path, s) for s in archive_files_list] - for entry in archive_files: - if 'expression.tab' in entry or 'DataMTX.tab' in entry: - # Get columns and rows of expression data in list form. - exp = pd.read_csv(entry, sep='\t', index_col=0, header=0) - exp_obs = list(exp.columns) - exp_genes= list(exp.index) - adata = sc.read(entry, first_column_names=True, cache=False).transpose() - elif 'observations.tab' in entry or 'COLmeta.tab' in entry: - obs = pd.read_csv(entry, sep='\t', index_col=0, header=0) - elif 'genes.tab' in entry or 'ROWmeta.tab' in entry: - var = pd.read_csv(entry, sep='\t', index_col=0, header=0) - if var.index.str.contains('ENS').sum() >1: - is_en=True - elif 'EXPmeta' in entry: - continue - ### Needs to be changed to account for other file types PCA etc - else: - raise Exception("Unexpected file name: '{0}'.".format(entry)) - - for str_type in ['cell_type', 'condition', 'time_point', 'time_unit']: - if str_type in obs.columns: - obs[str_type] = pd.Categorical(obs[str_type]) - for num_type in ['replicate', 'time_point_order']: - if num_type in obs.columns: - obs[num_type] = pd.to_numeric(obs[num_type]) - - # Ensure observations and genes are sorted the same as found in the expressions file - obs_index = list(obs.index) - if set(obs_index) != set(exp_obs): - raise Exception("Observation IDs from 'expressions' and 'observations' files are not the same.") - obs = obs.reindex(exp_obs) - - genes_index = list(var.index) - if set(genes_index) != set(exp_genes): - raise Exception("Gene IDs from 'expressions' and 'genes' files are not the same.") - var = var.reindex(exp_genes) - - adata.var = var - adata.obs = obs - self.adata = adata - self.originalPath = data_path - - return self, is_en - - def write_h5ad(self, output_path, gear_identifier): - """ - Writes anndata object as h5ad. Returns path to file. - - Input - ----- - output_path - path/to/write/converted/file - gear_identifier - unique name for output file - - Output - ------ - path to successfully generated h5ad format file. - - """ - if self.adata is None: - raise Exception("No AnnData object present to write to file.") - if output_path is None: - raise Exception("No destination file path given. Provide one to write file.") - try: - self.adata.write(filename = output_path) - except Exception as err: - raise Exception("Error occurred while writing to file: ", err) - return "" - - - diff --git a/lib/gear/datasetstats.py b/lib/gear/datasetstats.py deleted file mode 100755 index 14bbd1b28..000000000 --- a/lib/gear/datasetstats.py +++ /dev/null @@ -1,442 +0,0 @@ -import numpy as np -import pandas as pd - -""" -Class and associated methods needed to calculate the following statistics at time of upload: - - replicate averages - - p-values - - FDR (corrected p-values) - - standard deviations - - standard errors (of the mean) - - NOTE: - NumPy stats: https://docs.scipy.org/doc/numpy-dev/reference/routines.statistics.html - SciPy stats: https://docs.scipy.org/doc/scipy/reference/stats.html - - Averages: np.nanmean - https://docs.scipy.org/doc/numpy-dev/reference/generated/numpy.nanmean.html - standard_deviations: np.nanstd - https://docs.scipy.org/doc/numpy-dev/reference/generated/numpy.nanstd.html - standard_errors (standard error of the mean): https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.sem.html - -""" - - -def check_all_nans(values, total_rep): - """ - Checks to see if numpy array contains only nan values. Returns True - - Input: numpy array of replicate values - - Output: True/False - True if all replicates are nan - False if at least 1 value is present - """ - - all_nans = False - if np.isnan(values).any(): - #Only do this if an nan is present. - count_nans = 0 - for value in values: - if np.isnan(value): - count_nans +=1 - - if count_nans == total_rep: - all_nans = True - - return all_nans - -def get_replicate_counts(obs=None): - """ - Input: AnnData.obs observations (pandas dataframe) - - Output: python list containing the number of replicates there are for each condition. - - Example: Gene 'y' has 3 conditions. If Condition 1 has 2 replicates, Condition 2 has 3 replicates, - and Condition 3 has 3 replicates, the list [2,3,3] is returned. - """ - - # This my 1st attempt. It worked, but it didnt utilize the power of python pandas - - # # Get replicate column - # replicate_col = 'replicate' - # replicates = obs.loc[:, replicate_col] - # - # # Build list of total number of replicates for each group - # replicate_count = list() - # prior_rep = None - # for r in replicates: - # if prior_rep is None or r > prior_rep: - # prior_rep = r - # else: - # replicate_count.append(prior_rep) - # prior_rep = None - # - # # Append the last one - # replicate_count.append(prior_rep) - # - # # Return the count list - # return replicate_count - - # Get replicate counts by grouping them with pandas built-ins - if 'time_point' in obs.columns: - groups = obs.groupby(['cell_type', 'condition', 'time_point']).count().reset_index() - else: - groups = obs.groupby(['cell_type', 'condition']).count().reset_index() - - #pandas series containing number of replicates for each condition - replicate_count = groups.loc[:, 'replicate'] - - return replicate_count - - -class DatasetStats: - # This class contains the methods needed to calculate various dataset statistics - - def get_replicate_averages(adata=None): - """ - Input is AnnData object - Output is AnnData.X np.ndarray named .Xmean - - .Xmean is the same shape as .X but contains the averaged replicate value in place of the individual replicate values - Example: - .X has 3 replicates for each condition (18 columns). 72, 92, 51 are 3 replicates for one condition. Their average is 71.6667. - .Xmean also has 18 columns. 71.6667 appears 3x. Once for each replicate - adata.X[:,0] = [ 72. 92. 51. 93. 1. 46. 0. 33. 46. 75. 56. 28. 90. 100. 7. 25. 40. 81.] - adata.Xmean[:,0] = [71.6667 71.6667 71.6667 46.6667 46.6667 46.6667 26.3333 26.3333 26.3333 53. 53. 53. 65.6667 65.6667 65.6667 48.6667 48.6667 48.6667] - """ - # numpy stats: https://docs.scipy.org/doc/numpy/reference/routines.statistics.html - - X = adata.X - col_count, row_count = X.shape - - # Get replicate column - replicate_count = get_replicate_counts(adata.obs) - - average_list = list() - for row in range(row_count): - replicate_vals = X[:,row] - averages = list() - start = 0 - end = 0 - # for total_rep in replicate_count: - for i, total_rep in replicate_count.iteritems(): - end += total_rep - - # Subsample the replicates from 1 condition - values = replicate_vals[start:end] - - # Is there an expression value? Skip calculation, if none found - all_nans = check_all_nans(values, total_rep) - if all_nans: - rep_avg = None - else: - # Calcuate average. Ignoring NaN values and round to 4 decimals - rep_avg = round( np.nanmean(values), 4 ) - - # Add the average for each replicate - for r in range(total_rep): - averages.append(rep_avg) - r += 1 - - # Move starting point to next set of replicates - start += total_rep - - average_list.append(averages) - row += 1 - - # Convert list of lists to numpy array - Xmean = np.array(average_list) - # Xmean = pd.DataFrame(np.asarray(average_list)) - - # Transpose back to match the original .X - return Xmean.T - - - def get_replicate_std(adata=None): - """ - Input is AnnData object - Output is AnnData.X np.ndarray named .Xstd - - .Xstd is the same shape as .X but contains the standard deviation value in place of the individual replicate values - Example: - .X has 3 replicates for each condition (18 columns). 72, 92, 51 are 3 replicates for one condition. Their standard deviation is 16.7398. - .Xstd also has 18 columns. 16.7398 appears 3x. Once for each replicate - adata.X[:,0] = [ 72. 92. 51. 93. 1. 46. 0. 33. 46. 75. 56. 28. 90. 100. 7. 25. 40. 81.] - adata.Xstd[:,0] = [16.7398 16.7398 16.7398 37.5618 37.5618 37.5618 19.362 19.362 19.362 19.3046 19.3046 19.3046 41.684 41.684 41.684 23.669 23.669 23.669 ] - https://docs.scipy.org/doc/numpy-dev/reference/generated/numpy.nanstd.html - """ - - X = adata.X - col_count, row_count = X.shape - - # Get replicate column - replicate_count = get_replicate_counts(adata.obs) - - std_list = list() - for row in range(row_count): - replicate_vals = X[:,row] - stds = list() - start = 0 - end = 0 - # for total_rep in replicate_count: - for i, total_rep in replicate_count.iteritems(): - end += total_rep - - # Subsample the replicates from 1 condition - values = replicate_vals[start:end] - - # Is there an expression value? Skip calculation, if none found - all_nans = check_all_nans(values, total_rep) - if all_nans: - rep_std = None - else: - # Calcuate standard deviations. Ignoring NaN values and round to 4 decimals - rep_std = round( np.nanstd(values), 4 ) - - # Add the std for each replicate - for r in range(total_rep): - stds.append(rep_std) - r += 1 - - # Move starting point to next set of replicates - start += total_rep - - std_list.append(stds) - row += 1 - - # Convert list of lists to numpy array - Xstd = np.array(std_list) - # Xstd = pd.DataFrame(np.asarray(std_list)) - - # Transpose back to match the original .X - return Xstd.T - - - def get_replicate_sem(adata=None): - """ - Input is AnnData object - Output is AnnData.X np.ndarray named .Xsem - - .Xsem is the same shape as .X but contains the standard error of the mean (SEM) in place of the individual replicate values - Example: - .X has 3 replicates for each condition (18 columns). 72, 92, 51 are 3 replicates for one condition. Their SEM is 11.8369. - .Xsem also has 18 columns. 11.8369 appears 3x. Once for each replicate - adata.X[:,0] = [ 72. 92. 51. 93. 1. 46. 0. 33. 46. 75. 56. 28. 90. 100. 7. 25. 40. 81.] - adata.Xsem[:,0] = [11.8369 11.8369 11.8369 26.5602 26.5602 26.5602 13.691 13.691 13.691 13.6504 13.6504 13.6504 29.475 29.475 29.475 16.7365 16.7365 16.7365] - https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.sem.html - """ - import scipy.stats as stats - - X = adata.X - col_count, row_count = X.shape - - # Get replicate column - replicate_count = get_replicate_counts(adata.obs) - - sem_list = list() - for row in range(row_count): - replicate_vals = X[:,row] - sems = list() - start = 0 - end = 0 - # for total_rep in replicate_count: - for i, total_rep in replicate_count.iteritems(): - end += total_rep - - # Subsample the replicates from 1 condition - values = replicate_vals[start:end] - - # Is there an expression value? Skip calculation, if none found - all_nans = check_all_nans(values, total_rep) - if all_nans: - rep_sem = None - else: - # Calcuate SEM. Ignoring NaN values and round to 4 decimals - rep_sem = round( stats.sem(values, nan_policy="omit"), 4 ) - - # Add the SEM for each replicate - for r in range(total_rep): - sems.append(rep_sem) - r += 1 - - # Move starting point to next set of replicates - start += total_rep - - sem_list.append(sems) - row += 1 - - # Convert list of lists to numpy array - Xsem = np.array(sem_list) - # Xsem = pd.DataFrame(np.asarray(sem_list)) - - # Transpose back to match the original .X - return Xsem.T - - - def get_replicate_pvalue(adata=None): - """ - Input is AnnData object - Output is AnnData.X np.ndarray named .Xpval - - .Xpval is the same shape as .X but contains the p-value in place of the individual replicate values - Example: - .X has 3 replicates for each condition (18 columns). - Control group are 3 replicates = 72, 92, 51 - Treated group are 3 replicates = 1, 46, 0 - .Xpval also has 18 columns. The pvalue to the above example 0.6836206 appears 3x. Once for each replicate. - adata.X[:,0] = [ 72. 92. 51. 93. 1. 46. 0. 33. 46. 75. 56. 28. 90. 100. 7. 25. 40. 81.] - adata.Xpval[:,0] = [None None None 0.6836206 0.6836206 0.6836206 None None None 0.5876048 0.5876048 0.5876048 None None None 0.4909728 0.4909728 0.4909728] - - Display for Treated when comparing Control vs treated - One value for a set of replicates - - P-value is calculated from Pearson calculation using python scipy.stats: - https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html#scipy.stats.pearsonr - """ - import scipy.stats as stats - - X = adata.X - obs = adata.obs - - # Group conditions to get total number of replicates for each - grouped = obs.groupby(['cell_type', 'condition', 'time_point']).count() - grouped.reset_index().T - - # Standardized keywords that denote the control group. - control_string = ['control', 'ctrl', 'wildtype', 'wt', 'input'] - - p_values_list = list() - for g, gene in enumerate(X.T): - p_values = [np.nan] * len(gene) - cell_type = None - time_point = None - control_list = list() - treated_list = list() - treated_indices = list() #track index position so p-value can later be inserted - count_replicates = 0 - count_conditions = 1 # adds control group here - - #Go through each condition for each expression row - for col, replicate in obs.iterrows(): - cell_type = replicate['cell_type'] - condition = replicate['condition'] - time_point = replicate['time_point'] - - # Replicate is control group - if any(condition.lower() in c for c in control_string): - # Add replicate to list - if len(control_list) < grouped.loc[cell_type, condition, time_point]['replicate']: - # Still adding replicates of current condition - control_list.append(gene[count_replicates]) - - # Replicate is non-control - else: - # Do we already get all the replicates for this condition? - if len(treated_list) < grouped.loc[cell_type, condition, time_point]['replicate']: - # Still adding replicates of current condition - treated_list.append(gene[count_replicates]) - treated_indices.append(count_replicates) - - # All replicates accounted for. Get p-value then reset - if len(control_list) == grouped.loc[cell_type, condition, time_point]['replicate'] and \ - len(treated_list) == grouped.loc[cell_type, condition, time_point]['replicate']: - - # Calculate pvalue - p_coef, p_value = stats.pearsonr(control_list, treated_list) - p_value = round(p_value, 7) - - #Control replicates remain None, treated replicates get the pvalue - for t in range(len(treated_list)): - p_values[treated_indices[t]] = p_value - - # Reset. Allows for additional non-control groups to be compared to the current control - treated_list = list() - treated_indices = list() - if count_conditions == grouped.index.levels[1].size: - # P-values for this set of conditions are calculated - # Reset control group for next condition - control_list = list() - - count_replicates += 1 - - p_values_list.append(p_values) - - # Convert list of lists to numpy array - Xpval = np.array(p_values_list) - - # Transpose back to match the original .X - return pd.DataFrame(Xpval.T) - - - def get_replicate_fdr(adata): - """ - Input is AnnData object - Output is AnnData.X np.ndarray named .Xfdr - - .Xfdr is the same shape as .X but contains the FDR (false discovery rate) in place of the individual replicate values - Example: - .X has 3 replicates for each condition (18 columns). - .Xpval also has 18 columns. The pvalue to the above example 0.6836206 appears 3x. Once for each replicate. - adata.Xpval[:,0] = [None None None 0.6836206 0.6836206 0.6836206 None None None 0.5876048 0.5876048 0.5876048 None None None 0.4909728 0.4909728 0.4909728] - adata.Xfdr[:,0] = [None None None 0.6836206 0.6836206 0.6836206 None None None 0.6836206 0.6836206 0.6836206 None None None 0.6836206 0.6836206 0.6836206] - - FDR calculation is based on the following: - https://stackoverflow.com/a/21739593/2900840 - """ - - # Get pvalue matrix - # Xpval = adata.Xpval.T - Xpval = adata.uns['Xpval'].T - - # Get replicate count - replicates = get_replicate_counts(adata.obs) - - fdrs_list = list() - for index, row in enumerate(Xpval): - current_fdrs = [np.nan] * len(row) - start = 0 - end = 0 - pvalue_list = list() - pvalue_indices = list() # tuple of int values (index_of_first_replicate, number_of_replicates) - for i, rep_count in enumerate(replicates): - end += rep_count - - # Take 1st replicate of group (skip None groups aka controls) - # and store starting index and number of replicates - if row[start:end][0] is not None: - pvalue_list.append(row[start:end][0]) - pvalue_indices.append((start, rep_count)) - - start += rep_count - - # The following block is heavily based on the FDR calculation provided - # here: https://stackoverflow.com/a/21739593/2900840 - # orignally written by user: emre (https://stackoverflow.com/users/955278/emre) - pvalue_list = [ (pvalue, i) for i, pvalue in enumerate(pvalue_list) ] - pvalue_list.sort() - pvalue_list.reverse() - n = len(pvalue_list) - fdr_list = list() - for i, vals in enumerate(pvalue_list): - rank = n - i - pvalue, index = vals - fdr_list.append( round((n/rank) * pvalue, 7) ) - for i in range(0, int(n)-1): - if fdr_list[i] < fdr_list[i+1]: - fdr_list[i+1] = fdr_list[i] - for i, vals in enumerate(pvalue_list): - pvalue, index = vals - fdr_list[index] = fdr_list[i] - - - # Add the FDR to fdr_list for each replicate. Control replicates remain None - for (i, rep_count), fdr in zip(pvalue_indices, fdr_list): - for t in range(rep_count): - current_fdrs[i+t] = fdr - - fdrs_list.append(current_fdrs) - - Xfdr = np.array(fdrs_list) - # Xfdr = pd.DataFrame(np.asarray(fdrs_list)) - - # Return the transposed np.array - return Xfdr.T diff --git a/lib/gear/datasetuploader.py b/lib/gear/datasetuploader.py deleted file mode 100755 index c89d5db3a..000000000 --- a/lib/gear/datasetuploader.py +++ /dev/null @@ -1,98 +0,0 @@ -import os, sys -import tarfile - -class FileType(object): - # filetypes = [] - def __init__(self, filetypes=None, adata=None): - self.filetypes = filetypes - self.adata = adata - - if self.filetypes is None: - self.filetypes = list() - - #TODO: .adata added via ExcelUploader child-class - if self.adata is None: - pass - -class DatasetUploader: - - def get_by_filetype(self, filetype=None, filepath=None): - """ - This factory nests the dataset filetype classes. Preventing them from being directly called. - - dataset = DatasetUploader.get_by_filetype('xlsx') - dataset.read_file(filepath) - - Creates an Excel class object which now can be used for process and upload the dataset. - - Follows example: - http://python-3-patterns-idioms-test.readthedocs.io/en/latest/Factory.html#preventing-direct-creation - """ - - if filetype == "xlsx" or filetype == "xls": - import gear.exceluploader as exceluploader - return exceluploader.ExcelUploader() - - if filetype == "tar" or filetype == "gz": - if filepath is None: - raise Exception("filepath is a required argument for tarball uploading") - else: - tartype = self.tarball_type(filepath) - if tartype == 'mex': - import gear.mexuploader as mexuploader - return mexuploader.MexUploader() - elif tartype == 'threetab': - import gear.threetabuploader as threetabuploader - return threetabuploader.ThreeTabUploader() - if filetype == "h5ad": - import gear.h5aduploader as h5aduploader - return h5aduploader.H5adUploader() - assert 0, "Do not recognize file type given: " + filetype - - def tarball_type(self, file=None): - """ - Since we allow users to upload tarballs (.tar or .tar.gz) we don't know - from the filename alone if this is MEX or ThreeTab. This method reads the contents - of the tarball and returns either 'mex' or 'threetab' based on what it finds. The basenames variable - contains the filenames in case the files are within a forlder inside the tarball. - - mex: - matrix.mtx - barcodes.tsv - genes.tsv - - threetab: - expression.tab - genes.tab - observations.tab - - None is returned if neither of these is true - - Added NEMO file format functionality. - DataMTX.tab -> expression.tab - COLmeta.tab -> observations.tab - ROWmeta.tab -> genes.tab - """ - t = tarfile.open(file, 'r') - filenames = t.getnames() - filestr=''.join(filenames) - basenames=[] - for f_path in filenames: - fname=os.path.basename(f_path) - basenames.append(fname) - - if 'expression.tab' in basenames and 'genes.tab' in basenames and 'observations.tab' in basenames: - return 'threetab' - - if 'matrix.mtx' in basenames and 'barcodes.tsv' in basenames and 'genes.tsv' in basenames: - return 'mex' - - if 'DataMTX.tab' in filestr and 'COLmeta.tab' in filestr and 'ROWmeta.tab' in filestr: - return 'threetab' - - # ? subject to change - if "spatial" in filestr and "tissue_positions_list.csv" in filestr and "filtered_feature_bc_matrix.h5" in filestr: - return 'visium' - if "spatial" in filestr and "tissue_positions.parquet" in filestr and "filtered_feature_bc_matrix.h5" in filestr: - return 'visiumhd' - return None diff --git a/lib/gear/exceluploader.py b/lib/gear/exceluploader.py deleted file mode 100755 index 3df9bcada..000000000 --- a/lib/gear/exceluploader.py +++ /dev/null @@ -1,276 +0,0 @@ -import anndata -import pandas as pd -import os, sys -import scanpy as sc -sc.settings.verbosity = 0 - -from gear.datasetuploader import FileType - - -class ExcelUploader(FileType): - """ - Called by datasetuploader.py (factory) when excel file is going uploaded - - Standardized names for different sheets: - 'expression' - contains the expression values matrix. - It is assumed the 1st column is the row index and the 1st row is the column index - 'observations' - contains details of the column index from 'expression'. - 'genes' - contains details of the row index from 'expression' sheet. This sheet is optional. - - Standardized names for columns in each sheet: - 'expression' - 'genes' - Must be the first column on the left. This can contain either: - 1) Ensembl ids, or; - 2) Gene symbols - 'observations' - - 'observations' - This column must contain the identical names of the first row from the 'expression' sheet. - 'cell_type' - This column is the anatomical region or cell_type. Examples: 'utricle', 'inner_hair_cell' - 'condition' - This column is the experimental condition of the observation. Examples: 'control', 'treated' - 'replicate' - This column is the number of replicate the observation is. If this is replicate 1 of 3, put '1'. - 'time_point' - This column is the time point, or age, at which the observation was taken. Examples: 'P0', '24' - 'time_point_order' - This column is a numeric value allowing for sorting of the time_point column values. - 'time_unit' - What is the unit of time used? hours? minutes? - """ - - def _read_file(self, filepath): - """ - Currently reads excel file 3x (once for each sheet) for validation purposes. In doing so, this allows for closer - validation and error handling. - - Input - ----- - filepath - /path/to/your/excel.xlsx - - Output - ------ - Output is assigned to the DatasetUploader object: - DatasetUploader.adata = adata - DatasetUploader.originalFile = filepath - - 'adata' is an AnnData object where data of each XLSX sheet is assigned to - the AnnData object: - AnnData.X = expression sheet - AnnData.obs = observations sheet - AnnData.var = genes sheet - (if absent, this only contains a pandas dataframe index of gene ids) - - 'filepath' is the file path of the original file - """ - - # Get the expression matrix - try: - exp_df = pd.read_excel(filepath, sheet_name='expression', index_col=0).transpose() - except: - raise Exception("No expression sheet found. Expected spreadsheet sheet named 'expression'.") - - try: - X = exp_df.values[:, 0:].astype(float) - except: - raise Exception("Encountered unexpected value type. Expected float type in expression matrix.") - - # Get counts of genes and observations - number_obs_from_exp, number_genes_from_exp = X.shape - - # Get the observations - try: - obs_df = pd.read_excel(filepath, sheet_name='observations', index_col='observations') - except: - raise Exception("No observations sheet found. Expected spreadsheet sheet named 'observations'.") - - # Verify number observations equal those found in expression sheet - number_obs, number_cond = obs_df.shape - if number_obs != number_obs_from_exp: - raise Exception("Observations sheet error. Row count ({0}) in 'observations' sheet must match column count of 'expression' sheet({1}).".format(number_obs, number_obs_from_exp)) - - # Verify observations index matches expression sheet col index - if not obs_df.index.equals(exp_df.index): - raise Exception("Observations sheet error. The names and order of the index column in 'observations' sheet must match the columns of 'expression' sheet.") - - # Get the genes (if present), else use the .var from exp_df - try: - genes_df = pd.read_excel(filepath, sheet_name='genes', index_col=0, converters={'gene_symbol': str}) - except: - # With 'genes' sheet absent try to get the genes from 'expression' - try: - # read expression sheet. Set genes as the index and only parse 1 column - genes_df = pd.read_excel(filepath, sheet_name='expression', index_col=0, usecols=[0,1]) - - # remove the 1 parsed column. This leaves an empty dataframe - # with an index of gene ids to match other datasets. - genes_df = genes_df.drop(genes_df.columns[0], axis=1) - except Exception as err: - error = "No 'genes' sheet found. Tried using genes column from 'expression' sheet as .var, but " + str(err) - raise Exception(error) - - # Check for numeric gene symbols - if 'gene_symbol' in genes_df.columns: - digit_count = genes_df['gene_symbol'].str.isnumeric().sum() - if digit_count > 0: - raise Exception("Genes sheet error. {0} gene symbols are listed as numbers, not gene symbols".format(str(digit_count))) - else: - raise Exception("Failed to find gene_symbol column in genes tab") - - for str_type in ['cell_type', 'condition', 'time_point', 'time_unit']: - if str_type in obs_df.columns: - obs_df[str_type] = pd.Categorical(obs_df[str_type]) - - for num_type in ['replicate', 'time_point_order']: - if num_type in obs_df.columns: - obs_df[num_type] = pd.to_numeric(obs_df[num_type]) - - # Verify number genes equal those found in expression sheet - number_genes, number_columns = genes_df.shape - if number_genes != number_genes_from_exp: - raise Exception("Genes sheet error. Row count ({0}) in 'genes' sheet must match row count of 'expression' sheet({1}).".format(number_genes, number_genes_from_exp)) - - # Verify genes index matches expression sheet columns - if not genes_df.index.equals(exp_df.columns): - raise Exception("Genes sheet error. The names and order of the index column in 'genes' sheet must match the rows of 'expression' sheet.") - - # Create AnnData object and return it - self.adata = anndata.AnnData(X=X, obs=obs_df, var=genes_df) - self.originalFile = filepath - return self - - def get_plot_preview_expression(self): - """ - Searches the adata.X values for the first set of values that are all present and - non-zero. - - Returns tuple containing (gene_symbol, expression): - 1) Gene symbol (if present) or Ensembl id of gene. - 2) Pandas dataframe containing the expression values of the first - gene containing values for all observations and are non-zero. - The dataframe contains all columns of observation characteristics - from adata.obs. - """ - expression = None - index = None - gene = None - X = pd.DataFrame(self.adata.X) - var = self.adata.var - obs = self.adata.obs - - index_length = len(X.index) - for i in range(0, len(X.columns)): - # Find a gene whose values are all populated and non-zero - if pd.notnull(X).all(1).sum() == index_length: - index = i - - # Get the gene's expression and add the observation characteristics - expression = pd.DataFrame(X[i]) - expression.columns = ['raw_value'] - - for col in obs.columns: - expression[col] = obs[col].values - break - - #Get gene symbol (if present), otherwise use index label (Ensembl id) - if 'gene_symbol' in var.columns: - gene = var.iloc[index]['gene_symbol'] - else: - gene = var.index.values[index] - - # Return gene and it's expression dataframe - return gene, expression - - - def _add_calculated_values(self, adata): - """ - TODO: calculations in DatasetStats make inaccurate assumption when determining - the number of replicates within a condition (observation). - - **Correcting this is now hold until handling replicates observational and technical is addressed and finalized. ** - - Runs the statistical calculations on AnnData.X and stores them to AnnData.uns - - Input - ----- - AnnData object - - Output - ------ - AnnData object with statistics. statistics are stored as numpy arrays in - the structured AnnData dictionary: - AnnData.uns['Xmean'] contains a numpy array martix of means - - """ - from gear.datasetstats import DatasetStats - - # Calculate replicate averages. Return them as numpy array called .Xmean - Xmean = DatasetStats.get_replicate_averages(adata) - adata.uns['Xmean'] = Xmean - - # Calculate replicate standard deviations - Xstd = DatasetStats.get_replicate_std(adata) - adata.uns['Xstd'] = Xstd - - # Calculate replicate standard errors (of the mean) - Xsem = DatasetStats.get_replicate_sem(adata) - adata.uns['Xsem'] = Xsem - - # Calculate replicate p-value (uses Pearson from scipy.stats) - Xpval = DatasetStats.get_replicate_pvalue(adata) - adata.uns['Xpval'] = Xpval - - # TODO Needs refactored now that Xpval is stored in .uns['Xpval'] - # Calculate replicate FDR - # Xfdr = DatasetStats.get_replicate_fdr(adata) - # adata.uns['Xfdr'] = Xfdr - - self.adata = adata - return self - - - def _add_color_values(self, adata): - #calculate raw and absolute coloring - """ - TODO: calculations in DatasetColoring make inaccurate assumption when determining - the number of replicates within a condition (observation). - - Dataset scoped coloring needs completed - - **Correcting this is now hold until handling replicates observational and technical is addressed and finalized. ** - - """ - from gear.datasetcoloring import DatasetColoring - - # Calcuate gene level coloring - raw - XColGenRaw = DatasetColoring.get_gene_coloring(adata, color_mode='raw') - adata.uns['XColGenRaw'] = XColGenRaw - - # Calcuate gene level coloring - absolute - XColGenAbs = DatasetColoring.get_gene_coloring(adata, color_mode='absolute') - adata.uns['XColGenAbs'] = XColGenAbs - - # Calcuate tissue level coloring - raw - XColTisRaw = DatasetColoring.get_tissue_coloring(adata, color_mode='raw') - adata.uns['XColTisRaw'] = XColTisRaw - - # Calcuate tissue level coloring - absolute - XColTisAbs = DatasetColoring.get_tissue_coloring(adata, color_mode='absolute') - adata.uns['XColTisAbs'] = XColTisAbs - - # TODO - # get_color_dataset() - # get_color_abs_dataset() - - - def _write_to_h5ad(self, filepath=None): - #TODO this might not be used. It's only a anndata.write command - """ - Write AnnData object to file in h5ad format. - """ - - if self.adata is None: - raise Exception("No AnnData object present to write to file.") - if filepath is None: - raise Exception("No destination file path given. Provide one to write file.") - - - # write AnnData object to file - try: - self.adata.write(filename=filepath) - except Exception as err: - raise Exception("Error occurred while writing to file:\n", err) - - return self diff --git a/lib/gear/gearemail.py b/lib/gear/gearemail.py deleted file mode 100755 index e69de29bb..000000000 diff --git a/lib/gear/mexuploader.py b/lib/gear/mexuploader.py deleted file mode 100755 index 662fa8a61..000000000 --- a/lib/gear/mexuploader.py +++ /dev/null @@ -1,103 +0,0 @@ -import anndata -import pandas as pd -import numpy as np -import os, sys -import scanpy as sc -sc.settings.verbosity = 0 - -from gear.datasetuploader import FileType - - -class MexUploader(FileType): - """ - Called by datasetuploader.py (factory) when a TAR archive (mexfiles_are_inside.tar) is being uploaded. - - """ - - def _read_file(self, filepath): - """ - Reads in TAR archive that contains the 3 MEX files: matrix.mtx, barcodes.tsv, and genes.tsv. - Output is directed to the FileUploader object. - - Input - ----- - filepath - /path/to/your/upload.tar - - TAR archive must end with file extension '.tar' - - The must contain 3 files by the exact names as below: - matrix.mtx - barcodes.tsv - genes.tsv - - Output - ------ - Output is assigned to the DatasetUploader object: - DatasetUploader.adata = adata - DatasetUploader.originalFile = filepath - - 'adata' is an AnnData object where data of each MEX file is assigned to - the AnnData object: - AnnData.X = matrix.mtx - AnnData.obs = barcodes.tsv - AnnData.var = genes.tsv - - 'filepath' is the file path of the original file - - """ - import tarfile - - #Get tar filename so tmp directory can be assigned - tar_filename = filepath.rsplit('/', 1)[1].rsplit('.')[0] - tmp_dir = '/tmp/' + tar_filename - - # Open the archive and extract each file into the tmp directory - with tarfile.open(filepath) as tf: - for entry in tf: - - # Extract file into tmp dir - filepath = "{0}/{1}".format(tmp_dir, entry.name) - tf.extract(entry, path=tmp_dir) - - # Read each file (.mtx as adata.X and .tsv as pandas dataframes) - if entry.name == 'matrix.mtx' or os.path.basename(filepath)== 'matrix.mtx': - adata = sc.read(filepath, first_column_names=True, cache=False).transpose() - elif entry.name == 'barcodes.tsv' or os.path.basename(filepath)== 'barcodes.tsv': - obs = pd.read_csv(filepath, sep='\t', index_col=0, header=None, names=['observations']) - elif entry.name == 'genes.tsv' or os.path.basename(filepath)=='genes.tsv': - var = pd.read_csv(filepath, sep='\t', index_col=0, header=None, names=['genes', 'gene_symbol']) - else: - raise Exception("Unexpected file name: '{0}'. Excepted 'matrix.mtx', 'genes.tsv', or 'barcodes.tsv'.".format(entry.name)) - - # Assign genes and observations to AnnData object - adata.var = var - adata.obs = obs - - # Apply AnnData obj and filepath to uploader obj - self.adata = adata - self.originalFile = filepath - return self - - - def get_plot_preview_expression(self): - #TODO: Currently no plot types support this data type. - pass - - - def _write_to_h5ad(self, filepath=None): - #TODO this might not be used. It's only a anndata.write command - """ - Write AnnData object to file in h5ad format. - """ - - if self.adata is None: - raise Exception("No AnnData object present to write to file.") - if filepath is None: - raise Exception("No destination file path given. Provide one to write file.") - - - # write AnnData object to file - try: - self.adata.write(filename=filepath) - except Exception as err: - raise Exception("Error occurred while writing to file:\n", err) - - return self diff --git a/lib/gear/threetabuploader.py b/lib/gear/threetabuploader.py deleted file mode 100755 index 22e80584a..000000000 --- a/lib/gear/threetabuploader.py +++ /dev/null @@ -1,121 +0,0 @@ -import anndata -import pandas as pd -import os, sys -import scanpy as sc -import tarfile -sc.settings.verbosity = 0 - -from gear.datasetuploader import FileType - - -class ThreeTabUploader(FileType): - """ - Called by datasetuploader.py (factory) when tarball containing three tab-delimited files is going uploaded - - Standardized names for different files: - 'expression.tab' - contains the expression values matrix. - It is assumed the 1st column is the row index and the 1st row is the column index - 'observations.tab' - contains details of the column index from 'expression'. - 'genes.tab' - contains details of the row index from 'expression' sheet. Should at least have gene symbol - - Standardized names for columns in each file: - 'expression.tab' - 'genes' - Must be the first column on the left. This should be the Ensembl ID - 'gene_symbol' - Such as 'Rfx7' or 'TonB' - 'observations' - - 'observations' - This column must contain the identical names of the first row from the 'expression' sheet. - 'cell_type' - This column is the anatomical region or cell_type. Examples: 'utricle', 'inner_hair_cell' - 'condition' - This column is the experimental condition of the observation. Examples: 'control', 'treated' - 'replicate' - This column is the number of replicate the observation is. If this is replicate 1 of 3, put '1'. - 'time_point' - This column is the time point, or age, at which the observation was taken. Examples: 'P0', '24' - 'time_point_order' - This column is a numeric value allowing for sorting of the time_point column values. - 'time_unit' - What is the unit of time used? hours? minutes? - """ - - def _read_file(self, filepath): - """ - Input - ----- - filepath - /path/to/your/some.tar[.gz] - - Output - ------ - Output is assigned to the DatasetUploader object: - DatasetUploader.adata = adata - DatasetUploader.originalFile = filepath - - 'adata' is an AnnData object where data of each TAB file is assigned to - the AnnData object: - AnnData.X = expression file - AnnData.obs = observations file - AnnData.var = genes file - - 'filepath' is the file path of the original tarball - """ - - # Get tar filename so tmp directory can be assigned - tar_filename = filepath.rsplit('/', 1)[1].rsplit('.')[0] - tmp_dir = '/tmp/' + tar_filename - - with tarfile.open(filepath) as tf: - for entry in tf: - # Extract file into tmp dir - filepath = "{0}/{1}".format(tmp_dir, entry.name) - tf.extract(entry, path=tmp_dir) - - # Read each file as pandas dataframes - # Added functionality to read files inside folder within tarball and for NEMO three tab files - if entry.name == 'expression.tab' or os.path.basename(filepath)== 'expression.tab' or 'DataMTX.tab' in entry.name: - # Get columns and rows of expression data in list form. - exp = pd.read_table(filepath, sep='\t', index_col=0, header=0) - exp_obs = list(exp.columns) - exp_genes= list(exp.index) - - # Read in expressions as AnnData object - adata = sc.read(filepath, first_column_names=True, cache=False).transpose() - elif entry.name == 'observations.tab' or os.path.basename(filepath)== 'observations.tab' or 'COLmeta.tab' in entry.name: - obs = pd.read_table(filepath, sep='\t', index_col=0, header=0) - elif entry.name == 'genes.tab' or os.path.basename(filepath)== 'genes.tab' or 'ROWmeta.tab' in entry.name: - var = pd.read_table(filepath, sep='\t', index_col=0, header=0) - - for str_type in ['cell_type', 'condition', 'time_point', 'time_unit']: - if str_type in obs.columns: - obs[str_type] = pd.Categorical(obs[str_type]) - - for num_type in ['replicate', 'time_point_order']: - if num_type in obs.columns: - obs[num_type] = pd.to_numeric(obs[num_type]) - - - # Ensure observations and genes are sorted the same as found in the expressions file - obs_index = list(obs.index) - if set(obs_index) != set(exp_obs): - raise Exception("Observation IDs from 'expressions' and 'observations' files are not the same.") - obs = obs.reindex(exp_obs) - - genes_index = list(var.index) - if set(genes_index) != set(exp_genes): - raise Exception("Gene IDs from 'expressions' and 'genes' files are not the same.") - var = var.reindex(exp_genes) - - # Assign genes and observations to AnnData object - adata.var = var - adata.obs = obs - - # Apply AnnData obj and filepath to uploader obj - self.adata = adata - self.originalFile = filepath - return self - - - def _write_to_h5ad(self, filepath=None): - ###Added this subroutine from mexuploader as it was not writing to h5ad without it. - if self.adata is None: - raise Exception("No AnnData object present to write to file.") - if filepath is None: - raise Exception("No destination file path given. Provide one to write file.") - try: - self.adata.write(filename=filepath) - except Exception as err: - raise Exception("Error occurred while writing to file: ", err) - return self diff --git a/listeners/Dockerfile.anndata_upload b/listeners/Dockerfile.anndata_upload new file mode 100644 index 000000000..0e53d4f8b --- /dev/null +++ b/listeners/Dockerfile.anndata_upload @@ -0,0 +1,37 @@ +FROM python:3.14-slim + +# NOTE: Call in context of the gEAR root directory + +# For pybigwig dependency +ENV PYO3_USE_ABI3_FORWARD_COMPATIBILITY=1 + +WORKDIR /gEAR + +# Install system dependencies +RUN apt update && apt install -y \ + build-essential \ + curl \ + zlib1g-dev \ + && rm -rf /var/lib/apt/lists/* + +# Copy requirements and install Python dependencies (context is gEAR root) +COPY listeners/requirements.txt . +RUN pip install --no-cache-dir -r requirements.txt --prefer-binary \ + && rm -rf ~/.cache/pip + +# Install lib/gear as editable package +COPY lib/gear ./lib/gear +COPY lib/setup.py ./lib/ +#RUN pip install -e ./lib/ + +# Copy consumer script +RUN mkdir -p listeners +COPY listeners/anndata_upload_consumer.py ./listeners/ + +# Ensure log directory exists +RUN mkdir -p /var/log/gEAR_queue + +# Copy gear.ini +COPY gear.ini ./gear.ini + +CMD ["python3", "listeners/anndata_upload_consumer.py"] \ No newline at end of file diff --git a/listeners/Dockerfile.gosling_upload b/listeners/Dockerfile.gosling_upload index c6e0a857c..16455a68f 100644 --- a/listeners/Dockerfile.gosling_upload +++ b/listeners/Dockerfile.gosling_upload @@ -1,4 +1,7 @@ FROM python:3.14-slim + +# NOTE: Call in context of the gEAR root directory + # For pybigwig dependency ENV PYO3_USE_ABI3_FORWARD_COMPATIBILITY=1 diff --git a/listeners/Dockerfile.projectr b/listeners/Dockerfile.projectr index c0fd8de58..4cd798597 100644 --- a/listeners/Dockerfile.projectr +++ b/listeners/Dockerfile.projectr @@ -1,5 +1,7 @@ FROM python:3.14-slim +# NOTE: Call in context of the gEAR root directory + # For pybigwig dependency # TODO: Have specific requirements files for each listener ENV PYO3_USE_ABI3_FORWARD_COMPATIBILITY=1 diff --git a/listeners/anndata_upload_consumer.py b/listeners/anndata_upload_consumer.py new file mode 100644 index 000000000..209b7fa80 --- /dev/null +++ b/listeners/anndata_upload_consumer.py @@ -0,0 +1,172 @@ +#!/opt/bin/python3 + +""" +anndata_upload_consumer.py - RabbitMQ consumer for expression dataset upload jobs. + +Processes H5AD, 3-tab, Excel, and MEX format datasets. +""" + +import gc +import json +import os +import sys +import traceback +from pathlib import Path + +# Gear root for file operations +gear_root = Path(__file__).resolve().parents[1] +gear_lib = gear_root / "lib" +sys.path.insert(0, str(gear_lib)) + +from gear.serverconfig import ServerConfig # noqa: I001 + +servercfg = ServerConfig().parse() + +queue_name = "anndata_upload_jobs" +os.makedirs("/var/log/gEAR_queue", exist_ok=True) +logfile = f"/var/log/gEAR_queue/{queue_name}.log" +pid = os.getpid() + +user_upload_base = gear_root / 'www' / 'uploads' / 'files' + + +def _on_request(channel, method_frame, properties, body) -> None: + """Callback to handle new anndata upload job message.""" + from gear.anndata_processor import AnndataProcessor # noqa: E402 + + delivery_tag = method_frame.delivery_tag + deserialized_body = json.loads(body) + + job_id = deserialized_body["job_id"] + share_uid = deserialized_body["share_uid"] + dataset_uid = deserialized_body["dataset_uid"] + dataset_format = deserialized_body["dataset_format"] + perform_primary_analysis = deserialized_body.get("perform_primary_analysis", False) + + with open(logfile, "a") as fh: + print( + f"{pid} - [x] Received request for anndata job {job_id}", + flush=True, + file=fh, + ) + + if not user_upload_base.is_dir(): + print( + f"{pid} - ERROR: User upload base directory {user_upload_base} does not exist", + flush=True, + file=fh, + ) + channel.basic_nack(delivery_tag=delivery_tag, requeue=False) + return + + try: + # Infer staging_area from share_uid directory structure + staging_area = None + for session_dir in user_upload_base.iterdir(): + candidate = session_dir / share_uid + if candidate.is_dir(): + staging_area = candidate + break + + if not staging_area: + raise FileNotFoundError(f"Could not find staging area for {share_uid}") + + status_file = staging_area / "status.json" + + # Process the job + processor = AnndataProcessor( + job_id=job_id, + share_uid=share_uid, + staging_area=staging_area, + status_file=status_file, + dataset_uid=dataset_uid, + ) + + result = processor.process( + dataset_format=dataset_format, + perform_primary_analysis=perform_primary_analysis, + ) + + print(f"{pid} - Job {job_id}: {result['message']}", flush=True, file=fh) + channel.basic_ack(delivery_tag=delivery_tag) + + except Exception as e: + traceback.print_exc() + print(f"{pid} - Caught error '{str(e)}'", flush=True, file=fh) + channel.basic_nack(delivery_tag=delivery_tag, requeue=False) + finally: + gc.collect() + + +class Consumer: + """RabbitMQ consumer with automatic reconnection for anndata uploads.""" + + def __init__(self, host: str) -> None: + self._reconnect_delay = 0 + self.host = host + + import gearqueue # noqa: F401 + + self._consumer = gearqueue.AsyncConnection( + host=self.host, + publisher_or_consumer="consumer", + queue_name=queue_name, + on_message_callback=_on_request, + pid=pid, + logfile=logfile, + purge_queue=False, + ) + + def run(self) -> None: + """Run the consumer with automatic reconnection.""" + while True: + try: + self._consumer.run() + except KeyboardInterrupt: + self._consumer.stop() + break + self._maybe_reconnect() + + def _maybe_reconnect(self) -> None: + """Attempt reconnection with exponential backoff.""" + import time + + import gearqueue # noqa: F401 + + if self._consumer.should_reconnect: + self._consumer.stop() + reconnect_delay = self._get_reconnect_delay() + print( + f"{pid} - Reconnecting after {reconnect_delay} seconds", + flush=True, + ) + time.sleep(reconnect_delay) + self._consumer = gearqueue.AsyncConnection( + host=self.host, + publisher_or_consumer="consumer", + queue_name=queue_name, + on_message_callback=_on_request, + pid=pid, + logfile=logfile, + ) + + def _get_reconnect_delay(self) -> int: + """Calculate reconnect delay with exponential backoff.""" + if self._consumer.was_consuming: + self._reconnect_delay = 0 + else: + self._reconnect_delay += 1 + if self._reconnect_delay > 30: + self._reconnect_delay = 30 + return self._reconnect_delay + + +def main() -> None: + """Start the anndata processing consumer.""" + host = servercfg["dataset_uploader"]["queue_host"] + consumer = Consumer(host=host) + consumer.run() + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/www/api/resources/dataset_processing.py b/www/api/resources/dataset_processing.py index 44ee0fc40..80fd92323 100644 --- a/www/api/resources/dataset_processing.py +++ b/www/api/resources/dataset_processing.py @@ -1,5 +1,7 @@ import json import os +import subprocess +import sys from pathlib import Path import geardb @@ -16,11 +18,28 @@ def _check_process_running(status_data: dict) -> None: status_data: The status dictionary to update if process is not running """ process_id = status_data.get('process_id', -1) - if process_id > 0 and os.system(f'ps -p {process_id} > /dev/null') != 0: - status_data['status'] = 'error' - status_data['message'] = 'The processing step failed. Please contact the gEAR team.' - status_data['progress'] = 0 + if process_id == -1: + raise Exception("No process ID found in status data to check if process is running") + result = subprocess.run( + ['ps', '-p', str(process_id)], + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL + ) + if result.returncode != 0: + raise Exception(f"Process with ID {process_id} is not running") + + +def _check_job_running(status_data: dict) -> None: + """ + Verify that the job is still running; mark as error if not. + + Args: + status_data: The status dictionary to update if job is not running + """ + job_id = status_data.get('job_id', -1) + if job_id == -1: + raise Exception("No job ID found in status data to check if job is running") def _error_response(message: str, status_code: int) -> tuple: """ @@ -78,6 +97,20 @@ def post(self, share_uid): if current_status in 'complete': status_data['progress'] = 100 elif current_status == 'processing' and dataset_format not in ["gosling"]: - _check_process_running(status_data) + # ? Remove the "gosling" check + try: + _check_job_running(status_data) + except Exception as e: + print(f"Error checking job status: {e}", file=sys.stderr) + # if it doesn't work, check for a process ID (legacy implementation) + try: + _check_process_running(status_data) + except Exception as e: + print(f"Error checking process status: {e}", file=sys.stderr) + # If both checks fail, we assume the process has failed + status_data['status'] = 'error' + status_data['message'] = 'The processing step failed. Please contact the gEAR team.' + status_data['progress'] = 0 + return status_data, 400 return status_data, 200 diff --git a/www/cgi/check_dataset_processing_status.cgi b/www/cgi/check_dataset_processing_status.cgi index 15cbe0455..b987a12a7 100755 --- a/www/cgi/check_dataset_processing_status.cgi +++ b/www/cgi/check_dataset_processing_status.cgi @@ -9,7 +9,7 @@ making sure that process is in fact still running. Structure returned: { - "process_id": 1234, + "job_id": 1234, "status": "processing", "message": "Processing the dataset. This may take a while.", "progress": 0 @@ -20,7 +20,7 @@ Where status can be 'extracting', 'processing', 'error', or 'complete'. import cgi import json -import os +import subprocess import sys from pathlib import Path @@ -29,8 +29,28 @@ def main(): print('Content-Type: application/json\n\n') form = cgi.FieldStorage() - session_id = form.getvalue('session_id') - share_uid = form.getvalue('share_uid') + session_id = form.getfirst('session_id') + share_uid = form.getfirst('share_uid') + + if session_id is None: + status = { + "job_id": -1, + "status": "error", + "message": "No session_id provided.", + "progress": 0 + } + print(json.dumps(status)) + return status + + if share_uid is None: + status = { + "job_id": -1, + "status": "error", + "message": "No share_uid provided.", + "progress": 0 + } + print(json.dumps(status)) + return status user_upload_file_root = Path(__file__).resolve().parents[1] / 'uploads' / 'files' user_upload_file_base = user_upload_file_root / session_id / share_uid @@ -38,7 +58,7 @@ def main(): if not status_file.is_file(): status = { - "process_id": -1, + "job_id": -1, "status": "error", "message": "No status file found. Please upload a dataset first.", "progress": 0 @@ -57,13 +77,29 @@ def main(): return status if state == 'processing': - # Check if the process is still running - process_id = status.get('process_id', -1) - if process_id > 0: - # TODO: check that the process is the correct name too - if os.system(f'ps -p {process_id} > /dev/null') != 0: + job_id = status.get("job_id", -1) + if job_id > 0: + pass + + else: + # Check if the process is still running. + # Fallback for legacy uploads + try: + process_id = status.get('process_id', -1) + if process_id == -1: + raise Exception("No process ID found in status data to check if process is running") + + result = subprocess.run( + ['ps', '-p', str(process_id)], + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL + ) + if result.returncode != 0: + raise Exception(f"Process with ID {process_id} is not running") + + except Exception: status['status'] = 'error' - status['message'] = 'The processing step failed. Please contact the gEAR team.' + status['message'] = 'No process nor job ID found for a job in processing state. Please contact the gEAR team.' status['progress'] = 0 return status diff --git a/www/cgi/get_uploads_in_progress.cgi b/www/cgi/get_uploads_in_progress.cgi index 3b93ed2fe..afe68c204 100755 --- a/www/cgi/get_uploads_in_progress.cgi +++ b/www/cgi/get_uploads_in_progress.cgi @@ -15,8 +15,11 @@ Data structure returned: """ -import cgi, json -import os, sys +import cgi +import json +import os +import subprocess +import sys lib_path = os.path.abspath(os.path.join('..', '..', 'lib')) sys.path.append(lib_path) @@ -109,14 +112,24 @@ def main(): result['uploads'][-1]['status'] = 'processing' result['uploads'][-1]['load_step'] = 'process-dataset' - # Check if the process is still running - process_id = status_json.get('process_id', -1) + # Check job IDs + job_id = status_json.get('job_id', -1) + if job_id == -1: + # Check if the process is still running + # Legacy implementation + process_id = status_json.get('process_id', -1) - if process_id > 0: - # TODO: check that the process is the correct name too - if os.system(f'ps -p {process_id} > /dev/null') != 0: + if process_id == -1: result['uploads'][-1]['status'] = 'error' + sp_result = subprocess.run( + ['ps', '-p', str(process_id)], + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL + ) + if sp_result.returncode != 0: + result['uploads'][-1]['status'] = 'error' + elif processing_status == 'complete': result['uploads'][-1]['status'] = 'processed' result['uploads'][-1]['load_step'] = 'finalize-dataset' diff --git a/www/cgi/process_uploaded_expression_dataset.cgi b/www/cgi/process_uploaded_expression_dataset.cgi index 4fcb35bf7..3c4138d32 100755 --- a/www/cgi/process_uploaded_expression_dataset.cgi +++ b/www/cgi/process_uploaded_expression_dataset.cgi @@ -1,551 +1,167 @@ #!/opt/bin/python3 """ -Process the uploaded expression dataset, regardless of type. As the data are uploaded, -a process which can take hours, the following data structure is periodically -written to the same directory as the dataset: +Process the uploaded expression dataset by publishing to RabbitMQ queue. -status.json - -{ - "process_id": 1234, - "status": "processing", - "message": "Processing the dataset. This may take a while.", - "progress": 0 -} - -Where status can be 'uploaded', 'extracting', 'processing', 'error', or 'complete'. +This CGI endpoint queues dataset processing jobs for asynchronous handling +by the anndata_upload_consumer worker. """ import cgi +import configparser import json import os import sys -import zipfile +import uuid from pathlib import Path -import anndata -import pandas as pd -import scanpy as sc -from scipy import sparse - -# This has a huge dependency stack of libraries. Occasionally, one of them has methods -# which prints debugging information on STDOUT, killing this CGI. So here we redirect -# STDOUT until we need it. +# Redirect stdout to suppress debug output from dependencies original_stdout = sys.stdout sys.stdout = open(os.devnull, 'w') -lib_path = Path(__file__).resolve().parents[2] / 'lib' +gear_root = Path(__file__).resolve().parents[2] +lib_path = gear_root / 'lib' sys.path.append(str(lib_path)) + +_config = configparser.ConfigParser() +_config.read(gear_root / 'gear.ini') + import geardb -from gear.primary_analysis import add_primary_analysis_to_dataset, PrimaryAnalysisProcessingError +from gear.anndata_processor import process_anndata_synchronously, write_status from gear.spatialhandler import SPATIALTYPE2CLASS -from gear.utils import update_adata_with_ensembl_ids -share_uid = None -session_id = None user_upload_file_base = '../uploads/files' -# Initial status (will be updated as we progress) -status = { - "process_id": None, - "status": "extracting", - "message": "", - "progress": 0 -} +class QueueDisabledError(Exception): + """Custom exception to indicate that the queue is disabled in configuration.""" + pass -def main(): - result = {'success':0, "message": ""} - global share_uid - global session_id +def main() -> tuple: + """Queue expression dataset processing job.""" + result = {'success': 0, "message": ""} form = cgi.FieldStorage() - share_uid = form.getvalue('share_uid') - session_id = form.getvalue('session_id') - dataset_format = form.getvalue('dataset_format') - spatial_format = form.getvalue('spatial_format') # may be None + share_uid = form.getfirst('share_uid') + session_id = form.getfirst('session_id') + dataset_format = form.getfirst('dataset_format') + spatial_format = form.getfirst('spatial_format') # May be None + + # Validate required parameters + if share_uid is None: + result['message'] = 'Share UID is required.' + return result, 400 - if share_uid is None or session_id is None or dataset_format is None: - result['message'] = 'Missing one or more required parameters.' - return result + if session_id is None: + result['message'] = 'Session ID is required.' + return result, 400 + + if dataset_format is None: + result['message'] = 'Dataset format is required.' + return result, 400 user = geardb.get_user_from_session_id(session_id) if user is None: result['message'] = 'User ID not found. Please log in to continue.' - return result + return result, 401 - # values are mex_3tab, excel, rdata, h5ad - dataset_formats = ['mex_3tab', 'excel', 'rdata', 'h5ad', 'spatial'] - dataset_upload_dir = Path(user_upload_file_base) / session_id / share_uid + # Validate dataset format + supported_adata_formats = ['h5ad', 'mex_3tab', 'excel', 'rdata'] + spatial_formats = list(SPATIALTYPE2CLASS.keys()) - # quickly write the status so the page doesn't error out + if dataset_format not in supported_adata_formats and dataset_format != 'spatial': + result['message'] = f'Unsupported dataset format: {dataset_format}' + return result, 400 - status_file = Path(dataset_upload_dir) / 'status.json' - with open(status_file, 'w') as f: - f.write(json.dumps(status)) + if dataset_format == "spatial": + if spatial_format not in spatial_formats: + result['message'] = f'Invalid spatial format: {spatial_format}' + return result, 400 + + # Locate dataset upload directory + dataset_upload_dir = Path(user_upload_file_base) / session_id / share_uid - # if the upload directory doesn't exist, we can't process the dataset if not dataset_upload_dir.is_dir(): result['message'] = 'Dataset/directory not found.' - write_status(dataset_upload_dir, 'error', result['message']) - return result + return result, 404 - if dataset_format not in dataset_formats: - result['message'] = 'Unsupported dataset format.' - write_status(dataset_upload_dir, 'error', result['message']) - return result + job_id = str(uuid.uuid4()) - if dataset_format == "spatial": - if spatial_format not in SPATIALTYPE2CLASS: - result['message'] = 'Invalid spatial format specified.' - write_status(dataset_upload_dir, 'error', result['message']) - return result + # Initialize status file + status = { + "job_id": job_id, + "status": "queued", + "message": "Job queued for processing.", + "progress": 0, + } + status_file = dataset_upload_dir / 'status.json' + write_status(status_file, status) + # Load and update metadata metadata_file = dataset_upload_dir / 'metadata.json' if not metadata_file.is_file(): - write_status(dataset_upload_dir, 'error', "No metadata JSON file found.") - return - with open(metadata_file, 'r') as f: - metadata = json.load(f) - dataset_uid = metadata.get('dataset_uid', '') - dataset_type = metadata.get('dataset_type', '') - - # Update metadata for downstream uses - metadata["dataset_format"] = dataset_format - metadata["perform_primary_analysis"] = True if dataset_type in ['single-cell-rnaseq', 'spatial'] else False - - with open(metadata_file, 'w') as f: - json.dump(metadata, f, indent=4) - - # Since this process can take a while, we want to fork off of apache and continue - # processing in the background. - - ############################################### - # This is the fork off of apache - # https://stackoverflow.com/a/22181041/1368079 - # https://stackoverflow.com/questions/6024472/start-background-process-daemon-from-cgi-script - # https://groups.google.com/g/comp.lang.python/c/gSRnd0RoVKY?pli=1 - do_fork = False - if do_fork: - sys.stdout = original_stdout - result['success'] = 1 - print(json.dumps(result)) - - sys.stdout.flush() - sys.stdout.close() - - sys.stderr.flush() - sys.stderr.close() - - f = os.fork() - - if f != 0: - # Terminate the parent - sys.exit(0) - # PARENT DEAD - - # CHILD CONTINUES FROM HERE - - status['process_id'] = os.getpid() - - # new child command - if dataset_format == 'mex_3tab': - process_mex_3tab(dataset_upload_dir, metadata["perform_primary_analysis"]) - elif dataset_format == 'excel': - process_excel(dataset_upload_dir, metadata["perform_primary_analysis"]) - elif dataset_format == "h5ad": - process_h5ad(dataset_upload_dir, metadata["perform_primary_analysis"]) - elif dataset_format == "spatial": - process_spatial(dataset_upload_dir, spatial_format, metadata["perform_primary_analysis"]) - else: - result["success"] = 0 - result["message"] = f"Unsupported dataset format: {dataset_format}" - return result - - if metadata["perform_primary_analysis"]: - try: - result["success"] = add_primary_analysis_to_dataset(dataset_uid, share_uid, dataset_upload_dir, dataset_format) - except PrimaryAnalysisProcessingError as e: - write_status(dataset_upload_dir, 'error', f"Error during primary analysis: {str(e)}") - return result - - status["progress"] = 100 - write_status(dataset_upload_dir, 'complete', "Dataset processed successfully.") - - result["success"] = 1 - result["message"] = "Dataset processed successfully." - return result - -def process_h5ad(upload_dir: Path, perform_primary_analysis: bool) -> None: - """ - Processes an uploaded .h5ad (AnnData) file in the specified upload directory by performing the following steps: - 1. Reads the .h5ad file as an AnnData object. - 2. Sanitizes and categorizes the observation (obs) dataframe columns. - 3. Writes the sanitized data to a new .h5ad file. - 4. Replaces the original file with the sanitized version. - 5. Updates the processing status at each stage. - - Args: - upload_dir (Path): The directory containing the uploaded .h5ad file. - - Returns: - None - """ - # If the file is an h5ad, it should be formatted as an AnnData object already. - # But we still want to do some sanitization of the obs dataframe. - - # TODO: Read in chunks to save memory - - write_status(upload_dir, 'processing', 'Initializing dataset processing.') - - filepath = upload_dir / f"{share_uid}.h5ad" - adata = anndata.read_h5ad(filepath) - obs = adata.obs - - total_steps = 4 if perform_primary_analysis else 3 - step_counter = 1 - status["progress"] = int((step_counter / total_steps) * 100) - write_status(upload_dir, 'processing', 'Sanitizing AnnData object') - - categorize_observation_columns(obs) - adata.obs = sanitize_obs_for_h5ad(obs) - - if "gene_symbol" not in adata.var.columns: - # get organism_id by converting sample_taxid - metadata_file = upload_dir / 'metadata.json' - if not metadata_file.is_file(): - write_status(upload_dir, 'error', "No metadata JSON file found.") - return + result['message'] = "No metadata JSON file found." + return result, 400 + try: with open(metadata_file, 'r') as f: metadata = json.load(f) - sample_taxid = metadata.get("sample_taxid", None) - organism_id=geardb.get_organism_id_by_taxon_id(sample_taxid) - if not organism_id: - write_status(upload_dir, 'error', "Could not determine organism ID from sample taxonomic ID.") - return - - adata = update_adata_with_ensembl_ids(adata, organism_id, "UNMAPPED_") - - step_counter += 1 - status["progress"] = int((step_counter / total_steps) * 100) - write_status(upload_dir, 'processing', 'Writing sanitized data to new H5AD.') - - h5ad_path = upload_dir / f"{share_uid}.new.h5ad" - adata.write(h5ad_path) - - # Replace the original file with the sanitized one - filepath.unlink() # remove original - h5ad_path.rename(filepath) # rename new to original name - - step_counter += 1 - status["progress"] = int((step_counter / total_steps) * 100) - write_status(upload_dir, 'processing', f"Finished processing dataset. {'Performing primary analysis...' if perform_primary_analysis else ''}") - -def process_3tab(upload_dir: Path, perform_primary_analysis: bool) -> None: - import subprocess - - chunk_size = 500 - expression_matrix_path = None - obs = None - var = None - - write_status(upload_dir, 'processing', 'Initializing dataset processing.') - - for infile in upload_dir.iterdir(): - # skip any files beginning with a dot - if infile.name.startswith('.'): - continue - - filepath = "{0}/{1}".format(upload_dir, infile.name) - - # Read each file as pandas dataframes - if infile.name == 'expression.tab' or 'DataMTX.tab' in infile.name: - expression_matrix_path = filepath - elif infile.name == 'observations.tab' or 'COLmeta.tab' in infile.name: - obs = pd.read_table(filepath, sep='\t', index_col=0, header=0) - elif infile.name == 'genes.tab' or 'ROWmeta.tab' in infile.name: - var = pd.read_table(filepath, sep='\t', index_col=0, header=0) - - if obs is None: - write_status(upload_dir, 'error', "No observations file found. Expected observations.tab or COLmeta.tab.") - return - - if var is None: - write_status(upload_dir, 'error', "No genes file found. Expected genes.tab or ROWmeta.tab.") - return - if expression_matrix_path is None: - write_status(upload_dir, 'error', "No expression file found. Expected expression.tab or DataMTX.tab.") - return - - categorize_observation_columns(obs) - - # Read in expressions as AnnData object in a memory-efficient manner - adata = sc.AnnData(obs=var, var=obs) - reader = pd.read_csv(expression_matrix_path, sep='\t', index_col=0, chunksize=chunk_size) + dataset_uid = metadata.get('dataset_uid', '') + dataset_type = metadata.get('dataset_type', '') - # Count rows safely without shell execution (https://github.com/IGS/gEAR/security/code-scanning/229) - try: - result = subprocess.run( - ['/usr/bin/wc', '-l', str(expression_matrix_path)], - capture_output=True, - text=True, - check=True + # Determine if primary analysis should be performed + perform_primary_analysis = ( + dataset_type in ['single-cell-rnaseq', 'spatial'] and + dataset_format != 'spatial' # Exclude spatial from anndata processor ) - total_rows = int(result.stdout.split()[0]) - except (subprocess.CalledProcessError, ValueError, FileNotFoundError): - # Fallback to Python if wc fails - total_rows = sum(1 for _ in open(expression_matrix_path)) - 1 - - if perform_primary_analysis: - total_rows += 1 # account for the additional primary analysis step that will be performed after this - - expression_matrix = [] - rows_read = 0 - - ## Try to process the file the quickest way first, assuming things are peachy. Then, if not, - # do some checks and conversions (slower) as a backup. - try: - for chunk in reader: - rows_read += chunk_size - percentage = int((rows_read / total_rows) * 100) - expression_matrix.append(sparse.csr_matrix(chunk.values)) - - status['progress'] = percentage - message = f"Processed {rows_read}/{total_rows} expression matrix chunks ..." - write_status(upload_dir, 'processing', message) - - adata.X = sparse.vstack(expression_matrix) # type: ignore - except Exception: - #print(f"\nOriginal vstack failed: {e}") - #print("Retrying with per-chunk cleanup...") - - expression_matrix.clear() - chunk_shapes = [] - rows_read = 0 - status['progress'] = 0 - - # Re-open reader here - reader = pd.read_csv(expression_matrix_path, sep='\t', index_col=0, chunksize=chunk_size) - - for chunk_index, chunk in enumerate(reader, start=1): - try: - # Clean each cell: strip string values - chunk = chunk.replace(r'^\s+|\s+$', '', regex=True) - chunk = chunk.apply(lambda col: col.map(lambda x: x.strip() if isinstance(x, str) else x)) - - # Convert to numeric (non-numeric → NaN → fill with 0) - chunk_numeric = chunk.apply(pd.to_numeric, errors='coerce').fillna(0) - - matrix = sparse.csr_matrix(chunk_numeric.values) - expression_matrix.append(matrix) - chunk_shapes.append(matrix.shape) - - rows_read += chunk_size - percentage = int((rows_read / total_rows) * 100) - - message = f"Processed {rows_read}/{total_rows} expression matrix chunks ..." - write_status(upload_dir, 'processing', message) - - except Exception: - #print(f"\nError in chunk {chunk_index}: {inner_e}") - #print("Chunk head:") - #print(chunk.head()) - raise - - # Try stacking the cleaned chunks - try: - adata.X = sparse.vstack(expression_matrix) # type: ignore - except Exception: - #print(f"\nFinal vstack still failed: {final_e}") - - #print("Collected chunk shapes:") - #for i, shape in enumerate(chunk_shapes): - # print(f" Chunk {i+1}: {shape}") - raise - - adata = adata.transpose() - adata.obs = sanitize_obs_for_h5ad(adata.obs) - h5ad_path = upload_dir / f"{share_uid}.h5ad" - adata.write(h5ad_path) + metadata["dataset_format"] = dataset_format + metadata["perform_primary_analysis"] = perform_primary_analysis - # Progress is accounted for in chunk processing - write_status(upload_dir, 'processing', f"Finished processing dataset. {'Performing primary analysis...' if perform_primary_analysis else ''}") + with open(metadata_file, 'w') as f: + json.dump(metadata, f, indent=4) -def process_excel(upload_dir: Path, perform_primary_analysis: bool) -> None: - filepath = upload_dir / f"{share_uid}.xlsx" + except (json.JSONDecodeError, IOError) as e: + result['message'] = f"Error reading metadata: {str(e)}" + return result, 500 - write_status(upload_dir, 'processing', 'Initializing dataset processing.') + # Queue the job (skips spatial format) + if dataset_format == 'spatial': + if spatial_format is None: + result['message'] = 'Spatial format is required for spatial datasets.' + return result, 400 + process_spatial(job_id, dataset_upload_dir, share_uid, spatial_format, metadata["perform_primary_analysis"]) + return result, 200 - exp_df = pd.read_excel(filepath, sheet_name='expression', index_col=0).transpose() - - try: - X = exp_df.to_numpy()[:, 0:].astype(float) - except ValueError: - write_status(upload_dir, 'error', "Encountered unexpected value type. Expected float type in expression matrix.") - return - - # Get counts of genes and observations - number_obs_from_exp, number_genes_from_exp = X.shape - - # Get the observations - try: - obs_df = pd.read_excel(filepath, sheet_name='observations', index_col='observations') - except ValueError: - write_status(upload_dir, 'error', "No observations sheet found. Expected spreadsheet sheet named 'observations'.") - return - - # Verify number observations equal those found in expression sheet - number_obs, number_cond = obs_df.shape - if number_obs != number_obs_from_exp: - write_status(upload_dir, 'error', "Observations sheet error. Row count ({0}) in 'observations' sheet must match row count of 'expression' sheet({1}).".format(number_obs, number_obs_from_exp)) - return - - # Verify observations index matches expression sheet col index - if not obs_df.index.equals(exp_df.index): - write_status(upload_dir, 'error', "Observations sheet error. The names and order of the index column in 'observations' sheet must match the rows of 'expression' sheet.") - return - - # Get the genes (if present), else use the .var from exp_df + # Queue the job try: - genes_df = pd.read_excel(filepath, sheet_name='genes', index_col=0, converters={'gene_symbol': str}) - except ValueError: - # With 'genes' sheet absent try to get the genes from 'expression' - try: - # read expression sheet. Set genes as the index and only parse 1 column - genes_df = pd.read_excel(filepath, sheet_name='expression', index_col=0, usecols=[0,1]) - - # remove the 1 parsed column. This leaves an empty dataframe - # with an index of gene ids to match other datasets. - genes_df = genes_df.drop(genes_df.columns[0], axis=1) - except Exception as err: - write_status(upload_dir, 'error', "No 'genes' sheet found. Tried using genes column from 'expression' sheet as .var, but " + str(err)) - return - - # Check for numeric gene symbols - if 'gene_symbol' in genes_df.columns: - digit_count = genes_df['gene_symbol'].str.isnumeric().sum() - if digit_count > 0: - write_status(upload_dir, 'error', "Genes sheet error. {0} gene symbols are listed as numbers, not gene symbols".format(str(digit_count))) - return - else: - write_status(upload_dir, 'error', "Failed to find gene_symbol column in genes tab") - return + queue_anndata_job(job_id, share_uid, dataset_uid, dataset_format, metadata["perform_primary_analysis"]) + result["success"] = True + result["message"] = "Dataset upload processing job queued" + return result, 202 # Accepted + except QueueDisabledError: + + result_sync = process_anndata_synchronously( + job_id=job_id, + share_uid=share_uid, + staging_area=dataset_upload_dir, + status_file=status_file, + dataset_uid=dataset_uid, + dataset_format=dataset_format, + perform_primary_analysis=metadata["perform_primary_analysis"], - categorize_observation_columns(obs_df) - - # Verify number genes equal those found in expression sheet - number_genes, number_columns = genes_df.shape - if number_genes != number_genes_from_exp: - write_status(upload_dir, 'error', "Genes sheet error. Row count ({0}) in 'genes' sheet must match row count of 'expression' sheet({1}).".format(number_genes, number_genes_from_exp)) - return - - # Verify genes index matches expression sheet columns - if not genes_df.index.equals(exp_df.columns): - write_status(upload_dir, 'error', "Genes sheet error. The names and order of the index column in 'genes' sheet must match the rows of 'expression' sheet.") - return - - # Create AnnData object and return it - adata = anndata.AnnData(X=X, obs=obs_df, var=genes_df) - adata.obs = sanitize_obs_for_h5ad(adata.obs) - - h5ad_path = upload_dir / f"{share_uid}.h5ad" - adata.write(h5ad_path) - - total_steps = 2 if perform_primary_analysis else 1 - status["progress"] = int((1 / total_steps) * 100) - write_status(upload_dir, 'processing', f"Finished processing dataset. {'Performing primary analysis...' if perform_primary_analysis else ''}") - - -def process_mex(upload_dir: Path, perform_primary_analysis: bool) -> None: - pass - -def process_mex_3tab(upload_dir: Path, perform_primary_analysis: bool) -> None: - # Extract the file - import tarfile - compression_format = None - filename = upload_dir / f"{share_uid}.tar.gz" - - if filename.exists(): - compression_format = 'tarball' - else: - filename = upload_dir / f"{share_uid}.zip" + ) - if filename.exists(): - compression_format = 'zip' - else: - write_status(upload_dir, 'error', "No tarball or zip file found.") - return + result["success"] = result_sync["success"] + result["message"] = result_sync["message"] + return result, 200 if result["success"] else 500 - files_extracted = [] + except Exception as e: + result["message"] = f"Error processing track hub: {str(e)}" + print(f"AnndataUpload error: {str(e)}", file=sys.stderr) + return result, 500 - if compression_format == 'tarball': - try: - with tarfile.open(filename) as tf: - for entry in tf: - tf.extract(entry, path=upload_dir) - - # Nemo suffixes - nemo_suffixes = ['DataMTX.tab', 'COLmeta.tab', 'ROWmeta.tab'] - suffix_found = None - - for suffix in nemo_suffixes: - if entry.name.endswith(suffix): - suffix_found = suffix - # Rename the file to the appropriate name - old_name = upload_dir / entry.name - new_name = upload_dir / suffix - old_name.replace(new_name) - - if suffix_found is not None: - files_extracted.append(suffix_found) - else: - files_extracted.append(entry.name) - except tarfile.ReadError: - write_status(upload_dir, 'error', "Bad tarball file. Couldn't extract the tarball.") - return - - if compression_format == 'zip': - try: - with zipfile.ZipFile(filename) as zf: - for entry in zf.infolist(): - zf.extract(entry, path=upload_dir) - - # Nemo suffixes - nemo_suffixes = ['DataMTX.tab', 'COLmeta.tab', 'ROWmeta.tab'] - suffix_found = None - - for suffix in nemo_suffixes: - if entry.filename.endswith(suffix): - suffix_found = suffix - # Rename the file to the appropriate name - old_name = upload_dir / entry.filename - new_name = upload_dir / suffix - old_name.replace(new_name) - - if suffix_found is not None: - files_extracted.append(suffix_found) - else: - files_extracted.append(entry.filename) - except zipfile.BadZipFile: - write_status(upload_dir, 'error', "Bad zip file. Couldn't extract the zip file.") - return - - # Determine the dataset type - dataset_type = package_content_type(files_extracted) - - if dataset_type is None: - write_status(upload_dir, 'error', "Unsupported dataset format. Couldn't tell type from file names within the tarball") - - # Call the appropriate function - if dataset_type == 'threetab': - process_3tab(upload_dir, perform_primary_analysis) - elif dataset_type == 'mex': - process_mex(upload_dir, perform_primary_analysis) - -def process_spatial(upload_dir: Path, spatial_format: str, perform_primary_analysis: bool) -> None: +def process_spatial(job_id, upload_dir: Path, share_uid: str, spatial_format: str, perform_primary_analysis: bool) -> None: """ Processes a spatial transcriptomics dataset uploaded to a specified directory. @@ -563,12 +179,22 @@ def process_spatial(upload_dir: Path, spatial_format: str, perform_primary_analy Writes error status if the metadata file is missing or if reading/converting the spatial file fails. """ - write_status(upload_dir, 'processing', 'Initializing dataset processing.') + status = { + "job_id": job_id, + "status": "processing", + "message": "Initializing dataset processing.", + "progress": 0, + } + status_file = upload_dir / 'status.json' + write_status(status_file, status) spatial_obj = SPATIALTYPE2CLASS[spatial_format]() # instantiate the appropriate handler class metadata_file = upload_dir / 'metadata.json' if not metadata_file.is_file(): - write_status(upload_dir, 'error', "No metadata JSON file found.") + status["status"] = "error" + status["message"] = "No metadata JSON file found." + write_status(status_file, status) + return # get organism_id by converting sample_taxid(needed for some but not all spatial handlers) with open(metadata_file, 'r') as f: @@ -583,7 +209,9 @@ def process_spatial(upload_dir: Path, spatial_format: str, perform_primary_analy except Exception as e: import traceback traceback.print_exc() - write_status(upload_dir, 'error', f"Error in uploading spatial file: {e}") + status["status"] = "error" + status["message"] = f"Error in uploading spatial file: {e}" + write_status(status_file, status) return output_filename = f"{share_uid}.zarr" @@ -597,82 +225,61 @@ def process_spatial(upload_dir: Path, spatial_format: str, perform_primary_analy total_steps = 3 if perform_primary_analysis else 2 step_counter = 1 status["progress"] = int((step_counter / total_steps) * 100) - write_status(upload_dir, 'processing', 'Writing Zarr store') + write_status(status_file, status) spatial_obj.write_to_zarr(filepath=output_path) step_counter += 1 status["progress"] = int((step_counter / total_steps) * 100) - write_status(upload_dir, 'processing', f"Finished processing spatial dataset. {'Performing primary analysis...' if perform_primary_analysis else ''}") - -def sanitize_obs_for_h5ad(obs_df: pd.DataFrame) -> pd.DataFrame: - for col in obs_df.columns: - if obs_df[col].dtype == 'object': - obs_df[col] = obs_df[col].fillna('').astype(str) - return obs_df - -def write_status(upload_dir: Path, status_name: str, message: str) -> None: - status['status'] = status_name - status['message'] = message - with open(upload_dir / 'status.json', 'w') as f: - f.write(json.dumps(status)) - -def package_content_type(filenames: list[str]) -> str | None: - #print("DEBUG: filenames", file=sys.stderr, flush=True) - #print(filenames, file=sys.stderr, flush=True) - """ - mex: - matrix.mtx - barcodes.tsv - genes.tsv - - threetab: - expression.tab - genes.tab - observations.tab - - None is returned if neither of these is true - - Added NEMO file format functionality. - DataMTX.tab -> expression.tab - COLmeta.tab -> observations.tab - ROWmeta.tab -> genes.tab - """ - if 'expression.tab' in filenames and 'genes.tab' in filenames and 'observations.tab' in filenames: - return 'threetab' - - if 'matrix.mtx' in filenames and 'barcodes.tsv' in filenames and 'genes.tsv' in filenames: - return 'mex' - - if 'DataMTX.tab' in filenames and 'COLmeta.tab' in filenames and 'ROWmeta.tab' in filenames: - return 'threetab' - - return None - -def categorize_observation_columns(obs: pd.DataFrame) -> None: - """ - Categorizes and converts specific columns in the observation DataFrame. + write_status(status_file, status) - For each of the string-type columns ('cell_type', 'condition', 'time_point', 'time_unit') present in the DataFrame, - converts the column to a pandas Categorical type. For each of the numeric-type columns ('replicate', 'time_point_order') - present in the DataFrame, converts the column to a numeric type. +def queue_anndata_job(job_id: str, share_uid: str, dataset_uid: str, dataset_format: str, perform_primary_analysis: bool) -> None: + """Queue anndata processing job to RabbitMQ.""" - Parameters: - obs (pd.DataFrame): The observation DataFrame whose columns will be categorized or converted. + # If queue is not enabled, return False + if not _config.getboolean('dataset_uploader', 'queue_enabled', fallback=False): + print("Queue is disabled in configuration. Cannot queue trackhub job. Falling back to synchronous processing.", file=sys.stderr) + raise QueueDisabledError() - Returns: - None: The function modifies the DataFrame in place. - """ - for str_type in ['cell_type', 'condition', 'time_point', 'time_unit']: - if str_type in obs.columns: - obs[str_type] = pd.Categorical(obs[str_type]) + import gearqueue + host = _config["dataset_uploader"]["queue_host"] + + try: + # Connect as a blocking RabbitMQ publisher + connection = gearqueue.Connection( + host=host, publisher_or_consumer="publisher" + ) + except Exception as e: + print(f"Error connecting to RabbitMQ: {e}", file=sys.stderr) + raise Exception(f"Error connecting to RabbitMQ: {e}") + + with connection: + connection.open_channel() + + payload = { + "job_id": job_id, + "share_uid": share_uid, + "dataset_uid": dataset_uid, + "dataset_format": dataset_format, + "perform_primary_analysis": perform_primary_analysis, + } + + try: + connection.publish( + queue_name="anndata_upload_jobs", + message=payload, # method dumps JSON + ) + except Exception as e: + print(f"Error publishing message to RabbitMQ: {e}", file=sys.stderr) + raise + return - for num_type in ['replicate', 'time_point_order']: - if num_type in obs.columns: - obs[num_type] = pd.to_numeric(obs[num_type]) if __name__ == '__main__': - result = main() + result, code = main() sys.stdout = original_stdout - print('Content-Type: application/json\n\n', flush=True) - print(json.dumps(result), flush=True) + # Return JSON result and status code + print(f"Status: {code}", flush=True) + print('Content-Type: application/json', flush=True) + print('', flush=True) # blank line ends headers + print(json.dumps(result), flush=True) \ No newline at end of file diff --git a/www/cgi/store_expression_dataset.cgi b/www/cgi/store_expression_dataset.cgi index 54bb1f262..0ce29b2cb 100755 --- a/www/cgi/store_expression_dataset.cgi +++ b/www/cgi/store_expression_dataset.cgi @@ -19,10 +19,10 @@ import geardb def main(): print('Content-Type: application/json\n\n') form = cgi.FieldStorage() - session_id = form.getvalue('session_id') - share_uid = form.getvalue('share_uid') - dataset_format = form.getvalue('dataset_format') - spatial_format = form.getvalue('spatial_format') # may be None + session_id = form.getfirst('session_id') + share_uid = form.getfirst('share_uid') + dataset_format = form.getfirst('dataset_format') + spatial_format = form.getfirst('spatial_format') # may be None if not share_uid: # should never happen error_msg = f"Unexpected missing share_uid in store_expression_dataset.cgi. session_id={session_id!r}" @@ -84,7 +84,7 @@ def main(): result['message'] = 'Dataset file saved successfully.' status = { - "process_id": None, + "job_id": None, "status": "uploaded", "message": "The dataset has been uploaded and is pending processing", "progress": 0 diff --git a/www/cgi/upload_epigenetic_data.cgi b/www/cgi/upload_epigenetic_data.cgi deleted file mode 100755 index e0f68696a..000000000 --- a/www/cgi/upload_epigenetic_data.cgi +++ /dev/null @@ -1,61 +0,0 @@ -#!/opt/bin/python3 - -import cgi -import json -import os, sys - -# This has a huge dependency stack of libraries. Occasionally, one of them has methods -# which prints debugging information on STDERR, killing this CGI. So here we redirect -# STDOUT until we need it. -original_stdout = sys.stdout -sys.stdout = open(os.devnull, 'w') - -lib_path = os.path.abspath(os.path.join('..', '..', 'lib')) -sys.path.append(lib_path) -import geardb - -from gear.datasetuploader import FileType, DatasetUploader - -# This is another attempt to fix the STDOUT hack above in a cleaner way, but I couldn't -# get it to actually change the loglevel of another imported module. -#import logging -#logging.getLogger("matplotlib").setLevel(logging.WARNING) - -def main(): - user_upload_file_base = '../datasets_epigenetic/uploads/files' - - result = {'success':0 } - result['filename'] = filename - - form = cgi.FieldStorage() - filename = form.getvalue('filename') - - # This means the PHP upload failed - if filename is None: - result['message'] = "File upload failed. Try again and contact us if this continues." - print_and_go(json.dumps(result)) - - filepath = user_upload_file_base + "/" + filename - session_id = form.getvalue('session_id') - - user = geardb.get_user_from_session_id(session_id) - - if user is None: - result['message'] = 'User ID not found. Please log in to continue.' - else: - # gEAR User found. Read and Validate expression file - file_type = filename.rsplit('.', 1)[1] - - result['message'] = 'File uploaded - ' + str(filename) - result['success'] = 1 - - print_and_go(json.dumps(result)) - -def print_and_go(content): - sys.stdout = original_stdout - print('Content-Type: application/json\n\n', flush=True) - print(content) - sys.exit(0) - -if __name__ == '__main__': - main() diff --git a/www/cgi/validate_expression.cgi b/www/cgi/validate_expression.cgi deleted file mode 100755 index 161de6d23..000000000 --- a/www/cgi/validate_expression.cgi +++ /dev/null @@ -1,81 +0,0 @@ -#!/opt/bin/python3 - -import cgi -import json -import os, sys - -# This has a huge dependency stack of libraries. Occasionally, one of them has methods -# which prints debugging information on STDERR, killing this CGI. So here we redirect -# STDOUT until we need it. -original_stdout = sys.stdout -sys.stdout = open(os.devnull, 'w') - -lib_path = os.path.abspath(os.path.join('..', '..', 'lib')) -sys.path.append(lib_path) -import geardb - -from gear.datasetuploader import FileType, DatasetUploader - -# This is another attempt to fix the STDOUT hack above in a cleaner way, but I couldn't -# get it to actually change the loglevel of another imported module. -#import logging -#logging.getLogger("matplotlib").setLevel(logging.WARNING) - -def main(): - user_upload_file_base = '../uploads/files' - - result = {'success':0, 'shape': {}, 'preview_gene': '', 'plots': [] } - - form = cgi.FieldStorage() - filename = form.getvalue('filename') - - # This means the PHP upload failed - if filename is None: - result['message'] = "File upload failed. Try again and contact us if this continues." - print_and_go(json.dumps(result)) - - filepath = user_upload_file_base + "/" + filename - h5_out_path = user_upload_file_base + "/" + filename.replace('.', '_') + '.h5ad' - session_id = form.getvalue('session_id') - - user = geardb.get_user_from_session_id(session_id) - - if user is None: - result['message'] = 'User ID not found. Please log in to continue.' - else: - # gEAR User found. Read and Validate expression file - file_type = filename.rsplit('.', 1)[1] - - #Initize uploader for expression data (with factory) - dsu = DatasetUploader() - dataset_uploader = dsu.get_by_filetype(filetype=file_type, filepath=filepath) - try: - dataset_uploader._read_file(filepath) - - #Let's give the user some feedback on the shape of the dataset. - rows_X, cols_X = dataset_uploader.adata.X.shape - cols_obs, rows_obs = dataset_uploader.adata.obs.shape - cols_var, rows_var = dataset_uploader.adata.var.shape - - dataset_uploader.adata.write(h5_out_path) - - result['success'] = 1 - result['message'] = 'File successfully read.' - - result['shape']['X'] = {'rows': str(rows_X), 'cols': str(cols_X)} - result['shape']['obs'] = {'rows': str(rows_obs), 'cols': str(cols_obs)} - result['shape']['var'] = {'rows': str(rows_var), 'cols': str(cols_var)} - - except Exception as err: - result['message'] = str(err) - - print_and_go(json.dumps(result)) - -def print_and_go(content): - sys.stdout = original_stdout - print('Content-Type: application/json\n\n', flush=True) - print(content) - sys.exit(0) - -if __name__ == '__main__': - main() diff --git a/www/include/navigation_bar.html b/www/include/navigation_bar.html index 84dbc8007..60edfef81 100644 --- a/www/include/navigation_bar.html +++ b/www/include/navigation_bar.html @@ -68,7 +68,6 @@