Skip to content

Commit e7c5c11

Browse files
authored
Purge fast_bioservices (#246)
* feat: use bioservice's `MyGeneInfo` Signed-off-by: Josh Loecker <joshloecker@icloud.com> * refactor: convert from fast_bioservices to COMO's built-in pipeline Signed-off-by: Josh Loecker <joshloecker@icloud.com> * refactor: remove fast-bioservices from pyproject Signed-off-by: Josh Loecker <joshloecker@icloud.com> * chore: use bioservices over fast_bioservices Signed-off-by: Josh Loecker <joshloecker@icloud.com> * feat(test): added unit tests for identifier conversion pipeline Signed-off-by: Josh Loecker <joshloecker@icloud.com> --------- Signed-off-by: Josh Loecker <joshloecker@icloud.com>
1 parent eb9a551 commit e7c5c11

12 files changed

Lines changed: 882 additions & 180 deletions

main/como/combine_distributions.py

Lines changed: 34 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,8 @@
1818
_OutputCombinedSourceFilepath,
1919
_SourceWeights,
2020
)
21-
from como.utils import LogLevel, log_and_raise_error, _num_columns, get_missing_gene_data
21+
from como.pipelines.identifier import convert
22+
from como.utils import LogLevel, get_missing_gene_data, log_and_raise_error, num_columns
2223

2324

2425
def _combine_z_distribution_for_batch(
@@ -55,9 +56,13 @@ def _combine_z_distribution_for_batch(
5556
f"source: '{source.value}', "
5657
f"context: '{context_name}'"
5758
)
58-
if _num_columns(matrix) < 2:
59-
logger.trace(f"A single sample exists for batch '{batch.batch_num}'. Returning as-is because no additional combining can be done")
60-
return matrix
59+
if num_columns(matrix) < 2:
60+
logger.trace(
61+
f"A single sample exists for batch '{batch.batch_num}'. Returning as-is because no additional combining can be done"
62+
)
63+
with_batch_num = matrix.copy()
64+
with_batch_num.columns = [batch.batch_num]
65+
return with_batch_num
6166

6267
values = matrix.values
6368
weighted_matrix = np.clip(
@@ -116,7 +121,7 @@ def _combine_z_distribution_for_source(
116121
A pandas dataframe of the weighted z-distributions
117122
"""
118123
# If only one batch column exists, return as-is (rename like R path-through).
119-
if _num_columns(merged_source_data) <= 1:
124+
if num_columns(merged_source_data) <= 1:
120125
logger.warning("A single batch exists for this source; returning matrix as-is.")
121126
out_df = merged_source_data.copy()
122127
out_df.columns = ["combine_z"]
@@ -203,8 +208,10 @@ def _combine_z_distribution_for_context(
203208
z_matrices.append(matrix)
204209

205210
z_matrix = pd.concat(z_matrices, axis="columns", join="outer", ignore_index=False)
206-
if _num_columns(z_matrix) <= 1:
207-
logger.trace(f"Only 1 source exists for '{context}', returning dataframe as-is becuase no data exists to combine")
211+
if num_columns(z_matrix) <= 1:
212+
logger.trace(
213+
f"Only 1 source exists for '{context}', returning dataframe as-is becuase no data exists to combine"
214+
)
208215
z_matrix.columns = ["combine_z"]
209216
return z_matrix
210217

@@ -275,8 +282,26 @@ async def _begin_combining_distributions(
275282
batch_results: list[pd.DataFrame] = []
276283
for batch in batch_names[source.value]:
277284
batch: _BatchEntry
278-
matrix_subset = cast(pd.DataFrame, matrix[[GeneIdentifier.ensembl_gene_id.value, *batch.sample_names]])
279-
matrix_subset = matrix_subset.set_index(keys=[GeneIdentifier.ensembl_gene_id.value], drop=True)
285+
if isinstance(matrix, pd.DataFrame):
286+
matrix_subset = matrix[[GeneIdentifier.entrez_gene_id.value, *batch.sample_names]]
287+
matrix_subset = matrix_subset.set_index(keys=[GeneIdentifier.entrez_gene_id.value], drop=True)
288+
matrix_subset = matrix_subset.drop(columns=["gene_symbol", "ensembl_gene_id"], errors="ignore")
289+
elif isinstance(matrix, sc.AnnData):
290+
conversion = convert(ids=matrix.var_names.tolist(), taxon=taxon)
291+
conversion.reset_index(drop=False, inplace=True)
292+
matrix: pd.DataFrame = matrix.to_df().T
293+
matrix.reset_index(inplace=True, drop=False, names=["gene_symbol"])
294+
matrix: pd.DataFrame = matrix.merge(
295+
conversion, left_on="gene_symbol", right_on="gene_symbol", how="left"
296+
)
297+
matrix_subset: pd.DataFrame = matrix[[GeneIdentifier.entrez_gene_id.value, *batch.sample_names]]
298+
matrix_subset: pd.DataFrame = matrix_subset.set_index(
299+
keys=[GeneIdentifier.entrez_gene_id.value], drop=True
300+
)
301+
else:
302+
raise TypeError(
303+
f"Unexpected matrix type for source '{source.value}'; expected 'pandas.DataFrame' or 'scanpy.AnnData': {type(matrix)}"
304+
)
280305

281306
output_fp = output_filepaths[source.value].parent / f"{source.value}_batch{batch.batch_num}_combined_z_distrobution_{context_name}.csv"
282307
batch_results.append(

main/como/merge_xomics.py

Lines changed: 10 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,12 @@
11
from __future__ import annotations
22

3-
import asyncio
43
import re
54
import sys
65
from pathlib import Path
7-
from typing import TextIO, cast
6+
from typing import TextIO
87

98
import numpy as np
109
import pandas as pd
11-
from fast_bioservices.biothings.mygene import MyGene
1210
from loguru import logger
1311

1412
from como.combine_distributions import _begin_combining_distributions
@@ -24,7 +22,7 @@
2422
_SourceWeights,
2523
)
2624
from como.project import Config
27-
from como.utils import log_and_raise_error, read_file, set_up_logging, get_missing_gene_data, return_placeholder_data
25+
from como.utils import get_missing_gene_data, log_and_raise_error, read_file, return_placeholder_data, set_up_logging
2826

2927

3028
class _MergedHeaderNames:
@@ -348,17 +346,16 @@ async def _update_missing_data(input_matrices: _InputMatrices, taxon_id: int) ->
348346
"proteomics": [input_matrices.proteomics],
349347
}
350348
logger.trace(f"Gathering missing data for data sources: {','.join(key for key in matrix_keys if key is not None)}")
349+
# ruff: disable[E501]
351350
# fmt: off
352-
results = await asyncio.gather(
353-
*[
354-
# Using 'is not None' is required because the truth value of a Dataframe is ambiguous
355-
get_missing_gene_data(values=input_matrices.trna, taxon_id=taxon_id) if input_matrices.trna is not None else asyncio.sleep(0),
356-
get_missing_gene_data(values=input_matrices.mrna, taxon_id=taxon_id) if input_matrices.mrna is not None else asyncio.sleep(0),
357-
get_missing_gene_data(values=input_matrices.scrna, taxon_id=taxon_id) if input_matrices.scrna is not None else asyncio.sleep(0),
358-
get_missing_gene_data(values=input_matrices.proteomics, taxon_id=taxon_id) if input_matrices.proteomics is not None else asyncio.sleep(0),
359-
]
351+
results: tuple[pd.DataFrame | None, ...] = (
352+
get_missing_gene_data(values=input_matrices.trna, taxon_id=taxon_id) if input_matrices.trna is not None else None,
353+
get_missing_gene_data(values=input_matrices.mrna, taxon_id=taxon_id) if input_matrices.mrna is not None else None,
354+
get_missing_gene_data(values=input_matrices.scrna, taxon_id=taxon_id) if input_matrices.scrna is not None else None,
355+
get_missing_gene_data(values=input_matrices.proteomics, taxon_id=taxon_id) if input_matrices.proteomics is not None else None,
360356
)
361357
# fmt: on
358+
# ruff: enable[E501]
362359
for i, key in enumerate(matrix_keys):
363360
matrix_keys[key].append(results[i])
364361

@@ -591,7 +588,7 @@ async def merge_xomics( # noqa: C901
591588
log_location: str | TextIO = sys.stderr,
592589
):
593590
"""Merge expression tables of multiple sources (RNA-seq, proteomics) into one."""
594-
_set_up_logging(level=log_level, location=log_location)
591+
set_up_logging(level=log_level, location=log_location)
595592
logger.info(f"Starting to merge all omics data for context: '{context_name}'")
596593

597594
# fmt: off

main/como/pipelines/identifier.py

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
from collections.abc import Sequence
2+
from typing import Literal
3+
4+
import pandas as pd
5+
from bioservices.mygeneinfo import MyGeneInfo
6+
7+
__all__ = [
8+
"convert",
9+
"determine_gene_type",
10+
]
11+
12+
T_MG_SCOPE = Literal["entrezgene", "ensembl.gene", "symbol"]
13+
T_MG_TRANSLATE = Literal["entrez_gene_id", "ensembl_gene_id", "gene_symbol"]
14+
T_MG_RETURN = list[dict[T_MG_TRANSLATE, str]]
15+
16+
17+
def _get_conversion(info: MyGeneInfo, values: list[str], scope: T_MG_SCOPE, fields: str, taxon: str) -> T_MG_RETURN:
18+
value_str = ",".join(map(str, values))
19+
results = info.get_queries(query=value_str, dotfield=True, scopes=scope, fields=fields, species=taxon)
20+
if not isinstance(results, list):
21+
raise TypeError(f"Expected results to be a list, but got {type(results)}")
22+
if not isinstance(results[0], dict):
23+
raise TypeError(f"Expected each result to be a dict, but got {type(results[0])}")
24+
25+
data: T_MG_RETURN = []
26+
for result in results:
27+
ensembl = result.get("query" if scope == "ensembl.gene" else "ensembl.gene")
28+
entrez = result.get("query" if scope == "entrezgene" else "entrezgene")
29+
symbol = result.get("query" if scope == "symbol" else "symbol")
30+
data.append({"ensembl_gene_id": ensembl, "entrez_gene_id": entrez, "gene_symbol": symbol})
31+
return data
32+
33+
34+
def convert(ids: int | str | Sequence[int] | Sequence[str] | Sequence[int | str], taxon: int | str, cache: bool = True):
35+
"""Convert between genomic identifiers.
36+
37+
This function will convert between the following components:
38+
- Entrez Gene ID
39+
- Ensembl Gene ID
40+
- Gene Symbol
41+
42+
:param ids: IDs to be converted
43+
:param taxon: Taxonomic identifier
44+
:param: scope: The type of identifier provided in `ids`
45+
:param cache: Should local caching be used for queries
46+
:return: DataFrame with columns "entrez_gene_id", "ensembl_gene_id", and "gene_symbol"
47+
"""
48+
my_geneinfo = MyGeneInfo(cache=cache)
49+
chunk_size = 1000
50+
id_list = list(map(str, [ids] if isinstance(ids, (int, str)) else ids))
51+
chunks = list(range(0, len(id_list), chunk_size))
52+
53+
data_type = determine_gene_type(id_list)
54+
if not all(v == data_type[id_list[0]] for v in data_type.values()):
55+
raise ValueError("All items in ids must be of the same type (Entrez, Ensembl, or symbols).")
56+
57+
scope = next(iter(data_type.values()))
58+
fields = ",".join({"ensembl.gene", "entrezgene", "symbol"} - {scope})
59+
taxon_str = str(taxon)
60+
return pd.DataFrame(
61+
[
62+
row
63+
for i in chunks
64+
for row in _get_conversion(
65+
info=my_geneinfo,
66+
values=id_list[i : i + chunk_size],
67+
scope=scope,
68+
fields=fields,
69+
taxon=taxon_str,
70+
)
71+
]
72+
)
73+
74+
75+
def determine_gene_type(items: str | list[str], /) -> dict[str, T_MG_SCOPE]:
76+
"""Determine the genomic data type.
77+
78+
:param items: A string or list of strings representing gene identifiers.
79+
The function will determine whether each identifier is an Entrez Gene ID,
80+
Ensembl Gene ID, or a gene symbol based on its format.
81+
82+
:return: A dictionary mapping each input item to its determined type, which can be one of:
83+
- "entrez_gene_id": If the item consists solely of digits.
84+
- "ensembl_gene_id": If the item starts with "ENS" and is
85+
followed by a specific format (length greater than 11 and the last 11 characters are digits).
86+
- "gene_symbol": If the item does not match the above criteria, it is assumed to be a gene symbol.
87+
"""
88+
items = [items] if isinstance(items, str) else items
89+
90+
determine: dict[str, Literal["entrezgene", "ensembl.gene", "symbol"]] = {}
91+
for i in items:
92+
i_str = str(i).split(".")[0] if isinstance(i, float) else str(i)
93+
if i_str.isdigit():
94+
determine[i_str] = "entrezgene"
95+
elif i_str.startswith("ENS") and (len(i_str) > 11 and all(i.isdigit() for i in i_str[-11:])):
96+
determine[i_str] = "ensembl.gene"
97+
else:
98+
determine[i_str] = "symbol"
99+
100+
return determine

main/como/proteomics/Crux.py

Lines changed: 19 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
# ruff: noqa
22

3+
from typing import Type
34
import asyncio
45
import multiprocessing
56
import os
@@ -8,10 +9,11 @@
89
from multiprocessing.sharedctypes import Synchronized
910
from pathlib import Path
1011

12+
from bioservices.biodbnet import BioDBNet
1113
import numpy as np
1214
import pandas as pd
1315
import tqdm
14-
from fast_bioservices import BioDBNet
16+
1517

1618
from como.proteomics.FileInformation import FileInformation, clear_print
1719

@@ -128,7 +130,9 @@ def mzml_to_sqt(self) -> None:
128130
)
129131

130132
# Replace all "comet.*" in output directory with the name of the file being processed
131-
comet_files = [str(file) for file in os.listdir(file_information.sqt_base_path) if str(file).startswith("comet.")]
133+
comet_files = [
134+
str(file) for file in os.listdir(file_information.sqt_base_path) if str(file).startswith("comet.")
135+
]
132136
for file_name in comet_files:
133137
# Determine the old file path
134138
old_file_path: Path = Path(file_information.sqt_base_path, file_name)
@@ -240,15 +244,19 @@ def collect_uniprot_ids_and_ion_intensity(self) -> None:
240244

241245
# Assign the file_information intensity dataframe to the gathered values
242246
self._file_information[i].intensity_df = pd.DataFrame(average_intensities_dict)
243-
self._file_information[i].intensity_df = self._file_information[i].intensity_df.groupby("uniprot", as_index=False).mean()
247+
self._file_information[i].intensity_df = (
248+
self._file_information[i].intensity_df.groupby("uniprot", as_index=False).mean()
249+
)
244250

245251
async def _convert_uniprot_wrapper(self) -> None:
246252
"""This function is a multiprocessing wrapper around the convert_ids function"""
247253
values = [self.async_convert_uniprot(self._file_information[i]) for i in range(len(self._file_information))]
248254

249255
# Create a progress bar of results
250256
# From: https://stackoverflow.com/a/61041328/
251-
progress_bar = tqdm.tqdm(desc="Starting UniProt to Gene Symbol conversion... ", total=len(self._file_information))
257+
progress_bar = tqdm.tqdm(
258+
desc="Starting UniProt to Gene Symbol conversion... ", total=len(self._file_information)
259+
)
252260
for i, result in enumerate(asyncio.as_completed(values)):
253261
await result # Get result from asyncio.as_completed
254262
progress_bar.set_description(f"Working on {i + 1} of {len(self._file_information)}")
@@ -271,13 +279,11 @@ async def async_convert_uniprot(self, file_information: FileInformation) -> None
271279
loop = asyncio.get_event_loop()
272280
# async with self._semaphore:
273281
# gene_symbols: pd.DataFrame = await loop.run_in_executor(None, self._biodbnet.db2db, "UniProt Accession", "Gene Symbol", input_values)
274-
gene_symbols: pd.DataFrame = await loop.run_in_executor(
275-
None,
276-
self._biodbnet.db2db,
277-
input_values,
278-
"UniProt Accession",
279-
"Gene Symbol",
282+
gene_symbols = self._biodbnet.db2db(
283+
input_db="UniProt Accession", output_db="Gene Symbol", input_values=input_values
280284
)
285+
if not isinstance(gene_symbols, pd.DataFrame):
286+
raise TypeError(f"Expected gene_symbols to be a DataFrame, but got {type(gene_symbols)}")
281287

282288
# The index is UniProt IDs. Create a new column of these values
283289
gene_symbols["uniprot"] = gene_symbols.index
@@ -372,7 +378,9 @@ def split_abundance_values(self) -> None:
372378
# Create a new dataframe to split the S# columns from
373379
split_frame: pd.DataFrame = dataframe.copy()
374380
# Get the current S{i} columns in
375-
abundance_columns: list[str] = [column for column in split_frame.columns if re.match(rf"{cell_type}_S{i}R\d+", column)]
381+
abundance_columns: list[str] = [
382+
column for column in split_frame.columns if re.match(rf"{cell_type}_S{i}R\d+", column)
383+
]
376384
take_columns: list[str] = ["symbol"] + abundance_columns
377385
average_intensity_name: str = f"{cell_type}_S{i}"
378386

main/como/proteomics_gen.py

Lines changed: 11 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -71,15 +71,18 @@ async def load_gene_symbol_map(gene_symbols: list[str], entrez_map: Path | None
7171
biodbnet.services.settings.TIMEOUT = 60
7272
for i in range(0, len(gene_symbols), step_size):
7373
# Operations: Goes from gene_symbols=["A", "B;C", "D"] -> gene.split=[["A"], ["B", "C"], ["D"]] -> chain.from_iterable["A", "B", "C", "D"]
74-
chunk: list[str] = list(itertools.chain.from_iterable(gene.split(";") for gene in gene_symbols[i : i + step_size]))
75-
dataframes.append(
76-
biodbnet.db2db(
77-
input_values=chunk,
78-
input_db=Input.GENE_SYMBOL.value,
79-
output_db=[Output.GENE_ID.value, Output.ENSEMBL_GENE_ID.value],
80-
taxon=9606,
81-
)
74+
chunk: list[str] = list(
75+
itertools.chain.from_iterable(gene.split(";") for gene in gene_symbols[i : i + step_size])
8276
)
77+
result = biodbnet.db2db(
78+
input_db="Gene Symbol",
79+
output_db=["Gene ID", "Ensembl Gene ID"],
80+
input_values=chunk,
81+
taxon=9606,
82+
)
83+
if not isinstance(result, pd.DataFrame):
84+
raise TypeError(f"Got {type(result)}, expected pd.DataFrame")
85+
dataframes.append(result)
8386
df = pd.concat(dataframes, axis="columns")
8487
print(df)
8588
df.loc[df["gene_id"].isna(), ["gene_id"]] = np.nan

main/como/rnaseq_gen.py

Lines changed: 3 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,6 @@
1010
from typing import NamedTuple, TextIO, cast
1111

1212
import anndata as ad
13-
import boolean
1413
import numpy as np
1514
import numpy.typing as npt
1615
import pandas as pd
@@ -19,13 +18,13 @@
1918
import sklearn.neighbors
2019
from anndata.compat import XDataArray
2120
from anndata.experimental.backed import Dataset2D
22-
from fast_bioservices.pipeline import ensembl_to_gene_id_and_symbol, gene_symbol_to_ensembl_and_gene_id
2321
from loguru import logger
2422
from scipy import sparse
2523
from zfpkm import zFPKM, zfpkm_plot
2624

2725
from como.data_types import FilteringTechnique, LogLevel, RNAType
2826
from como.migrations import gene_info_migrations
27+
from como.pipelines.identifier import convert
2928
from como.project import Config
3029
from como.utils import log_and_raise_error, read_file, set_up_logging
3130

@@ -173,17 +172,15 @@ async def _build_matrix_results(
173172
raise TypeError("AnnData.var is expected to be a pandas.DataFrame")
174173

175174
matrix.var = matrix.var.reset_index(drop=False, names=["gene_symbol"])
176-
conversion = await gene_symbol_to_ensembl_and_gene_id(symbols=matrix.var["gene_symbol"].tolist(), taxon=taxon)
175+
conversion = convert(ids=matrix.var["gene_symbol"].tolist(), taxon=taxon)
177176
else:
178177
if "ensembl_gene_id" not in matrix.columns:
179178
log_and_raise_error(
180179
message="'ensembl_gene_id' column not found in the provided DataFrame.",
181180
error=ValueError,
182181
level=LogLevel.CRITICAL,
183182
)
184-
conversion: pd.DataFrame = await ensembl_to_gene_id_and_symbol(
185-
ids=matrix["ensembl_gene_id"].tolist(), taxon=taxon
186-
)
183+
conversion: pd.DataFrame = convert(ids=matrix["ensembl_gene_id"].tolist(), taxon=taxon)
187184
# If the entrez gene id column is empty, it is indicative that the incorrect taxon id was provided
188185
if conversion["entrez_gene_id"].eq("-").all():
189186
logger.critical(

0 commit comments

Comments
 (0)