457 rows × 3 columns
\n", + "454 rows × 3 columns
\n", "" ], "text/plain": [ @@ -1095,16 +1029,16 @@ "3 6 170805979 0\n", "4 3 198295559 0\n", ".. ... ... ...\n", - "452 LRG_741 231167 0\n", - "453 LRG_763 176286 0\n", - "454 LRG_792 42144 0\n", - "455 LRG_793 38439 0\n", - "456 LRG_93 22459 0\n", + "449 LRG_311 115492 0\n", + "450 LRG_721 33396 0\n", + "451 LRG_741 231167 0\n", + "452 LRG_763 176286 0\n", + "453 LRG_93 22459 0\n", "\n", - "[457 rows x 3 columns]" + "[454 rows x 3 columns]" ] }, - "execution_count": 10, + "execution_count": 28, "metadata": {}, "output_type": "execute_result" } @@ -1117,123 +1051,52 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Adding annotations to an AnnData object:" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [], - "source": [ - "import pandas as pd\n", - "import scanpy as sc\n", - "pbmc = sc.datasets.pbmc3k()" + "## Adding annotations to an AnnData object\n", + "\n", + "We provide a utility function `annotate` to merge a table of feature annotations from an scverse dataset and a gene annotation table. In this example we use subsampled scRNA-seq data from [Litviňuková et al. 2020](https://www.nature.com/articles/s41586-020-2797-4), where reads were mapped to the hg38 reference genome with [CellRanger v3.0.1](https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/3.0), which uses the Ensembl release 98 for gene annotations." ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 14, "metadata": {}, "outputs": [ { - "data": { - "text/html": [ - "| \n", - " | gene_ids | \n", - "
|---|---|
| index | \n", - "\n", - " |
| MIR1302-10 | \n", - "ENSG00000243485 | \n", - "
| FAM138A | \n", - "ENSG00000237613 | \n", - "
| OR4F5 | \n", - "ENSG00000186092 | \n", - "
| RP11-34P13.7 | \n", - "ENSG00000238009 | \n", - "
| RP11-34P13.8 | \n", - "ENSG00000239945 | \n", - "
26662 rows × 11 columns
\n", "" ], "text/plain": [ - " orig_gene_name gene_name gene_biotype gene_seq_start \n", - "gene_id \n", - "ENSG00000243485 MIR1302-10 MIR1302-2HG lncRNA 29554.0 \\\n", - "ENSG00000237613 FAM138A FAM138A lncRNA 34554.0 \n", - "ENSG00000186092 OR4F5 OR4F5 protein_coding 65419.0 \n", - "ENSG00000238009 RP11-34P13.7 lncRNA 89295.0 \n", - "ENSG00000239945 RP11-34P13.8 lncRNA 89551.0 \n", + " gene_ids-Harvard-Nuclei gene_id gene_name \n", + "AL627309.1 ENSG00000238009 ENSG00000238009 AL627309.1 \\\n", + "AC114498.1 ENSG00000235146 ENSG00000235146 AC114498.1 \n", + "AL669831.2 ENSG00000229905 ENSG00000229905 AL669831.2 \n", + "AL669831.5 ENSG00000237491 ENSG00000237491 LINC01409 \n", + "FAM87B ENSG00000177757 ENSG00000177757 FAM87B \n", + "... ... ... ... \n", + "AC007325.2 ENSG00000277196 ENSG00000277196 AC007325.2 \n", + "BX072566.1 ENSG00000277630 ENSG00000277630 BX072566.1 \n", + "AL354822.1 ENSG00000278384 ENSG00000278384 AL354822.1 \n", + "AC004556.1 ENSG00000276345 ENSG00000276345 AC004556.3 \n", + "AC240274.1 ENSG00000271254 ENSG00000271254 AC240274.1 \n", "\n", - " gene_seq_end seq_name seq_strand seq_coord_system \n", - "gene_id \n", - "ENSG00000243485 31109.0 1 1.0 chromosome \\\n", - "ENSG00000237613 36081.0 1 -1.0 chromosome \n", - "ENSG00000186092 71585.0 1 1.0 chromosome \n", - "ENSG00000238009 133723.0 1 -1.0 chromosome \n", - "ENSG00000239945 91105.0 1 -1.0 chromosome \n", + " gene_biotype gene_seq_start gene_seq_end seq_name \n", + "AL627309.1 lncRNA 89295.0 133723.0 1 \\\n", + "AC114498.1 lncRNA 587629.0 594768.0 1 \n", + "AL669831.2 lncRNA 760911.0 761989.0 1 \n", + "AL669831.5 lncRNA 778747.0 810065.0 1 \n", + "FAM87B lncRNA 817371.0 819837.0 1 \n", + "... ... ... ... ... \n", + "AC007325.2 protein_coding 138082.0 161852.0 KI270734.1 \n", + "BX072566.1 protein_coding 108007.0 139659.0 GL000213.1 \n", + "AL354822.1 protein_coding 51867.0 54893.0 GL000218.1 \n", + "AC004556.1 protein_coding 2585.0 11802.0 KI270721.1 \n", + "AC240274.1 protein_coding 4612.0 29626.0 KI270711.1 \n", "\n", - " description \n", - "gene_id \n", - "ENSG00000243485 MIR1302-2 host gene [Source:HGNC Symbol;Acc:HG... \\\n", - "ENSG00000237613 family with sequence similarity 138 member A [... \n", - "ENSG00000186092 olfactory receptor family 4 subfamily F member... \n", - "ENSG00000238009 novel transcript \n", - "ENSG00000239945 novel transcript \n", + " seq_strand seq_coord_system \n", + "AL627309.1 -1.0 chromosome \\\n", + "AC114498.1 1.0 chromosome \n", + "AL669831.2 1.0 chromosome \n", + "AL669831.5 1.0 chromosome \n", + "FAM87B 1.0 chromosome \n", + "... ... ... \n", + "AC007325.2 -1.0 scaffold \n", + "BX072566.1 -1.0 scaffold \n", + "AL354822.1 -1.0 scaffold \n", + "AC004556.1 1.0 scaffold \n", + "AC240274.1 -1.0 scaffold \n", "\n", - " gene_id_version canonical_transcript \n", - "gene_id \n", - "ENSG00000243485 ENSG00000243485.5 ENST00000473358 \n", - "ENSG00000237613 ENSG00000237613.2 ENST00000417324 \n", - "ENSG00000186092 ENSG00000186092.7 ENST00000641515 \n", - "ENSG00000238009 ENSG00000238009.6 ENST00000477740 \n", - "ENSG00000239945 ENSG00000239945.1 ENST00000495576 " + " description \n", + "AL627309.1 novel transcript \\\n", + "AC114498.1 novel transcript \n", + "AL669831.2 novel transcript \n", + "AL669831.5 long intergenic non-protein coding RNA 1409 [S... \n", + "FAM87B family with sequence similarity 87 member B [S... \n", + "... ... \n", + "AC007325.2 proline dehydrogenase 1, mitochondrial [Source... \n", + "BX072566.1 POTE ankyrin domain family member D-like [Sour... \n", + "AL354822.1 NULL \n", + "AC004556.1 39S ribosomal protein L23, mitochondrial [Sour... \n", + "AC240274.1 neuroblastoma breakpoint family member 1 [Sour... \n", + "\n", + " gene_id_version \n", + "AL627309.1 ENSG00000238009.6 \n", + "AC114498.1 ENSG00000235146.2 \n", + "AL669831.2 ENSG00000229905.1 \n", + "AL669831.5 ENSG00000237491.10 \n", + "FAM87B ENSG00000177757.2 \n", + "... ... \n", + "AC007325.2 ENSG00000277196.4 \n", + "BX072566.1 ENSG00000277630.4 \n", + "AL354822.1 ENSG00000278384.1 \n", + "AC004556.1 ENSG00000276345.1 \n", + "AC240274.1 ENSG00000271254.6 \n", + "\n", + "[26662 rows x 11 columns]" ] }, - "execution_count": 14, + "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "pbmc.var.head()" + "annotated_var = gf.annotate(adata.var, genes, var_on='gene_ids-Harvard-Nuclei')\n", + "annotated_var" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 28, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "adata.var = annotated_var.copy()" + ] } ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:genomic-annotations-dev]", + "display_name": "Python 3", "language": "python", - "name": "conda-env-genomic-annotations-dev-py" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -1424,12 +1396,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.16" - }, - "vscode": { - "interpreter": { - "hash": "ae6466e8d4f517858789b5c9e8f0ed238fb8964458a36305fca7bddc149e9c64" - } + "version": "3.8.8" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index 8ce0ee5..466dd88 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,6 +53,7 @@ doc = [ test = [ "pytest", "pytest-cov", + "anndata" ] [tool.coverage.run] diff --git a/src/genomic_features/__init__.py b/src/genomic_features/__init__.py index 679dc43..e6440ee 100644 --- a/src/genomic_features/__init__.py +++ b/src/genomic_features/__init__.py @@ -1,6 +1,7 @@ from importlib.metadata import version from . import ensembl, filters +from .annotate import annotate __all__ = ["ensembl"] diff --git a/src/genomic_features/annotate.py b/src/genomic_features/annotate.py new file mode 100644 index 0000000..330f568 --- /dev/null +++ b/src/genomic_features/annotate.py @@ -0,0 +1,114 @@ +import warnings + +import pandas as pd + + +def annotate( + adata_var: pd.DataFrame, + annotation_df: pd.DataFrame, + var_on: str = None, + annotation_on: str = None, +) -> pd.DataFrame: + """ + Annotate variables from AnnData object with a DataFrame of annotations. + + Parameters + ---------- + adata_var + The AnnData.var DataFrame. + annotation_df + The DataFrame of annotations (e.g. output of `EnsemblDB.genes()`). + var_on + The column name in `adata_var` to join var_on (useful if geneIDs are not stored in `adata.var_names`). + If `None`, tables are joined based var_on index of `adata_var`. + annotation_on + The column name in `annotation_df` to join var_on. If `None`, the column name is inferred from the + `annotation_df` DataFrame (i.e. the column with unique values ending in `_id`). + + Returns + ------- + annotated_var + The annotated AnnData.var DataFrame. + """ + # Pick column with IDs in annotation table + if annotation_on is None: + annotation_on = [ + i + for i in annotation_df.columns + if annotation_df[i].is_unique and i.endswith("_id") + ] + if len(annotation_on) == 0: + raise ValueError( + "No unique ID column found in annotation_df - specify ID column with `annotation_on` parameter" + ) + if len(annotation_on) > 1: + raise ValueError( + "Multiple unique ID columns found in annotation_df - specify ID column with `annotation_on` parameter" + ) + else: + annotation_on = annotation_on[0] + else: + if not annotation_df[annotation_on].is_unique: + raise ValueError(f"Column {annotation_on} does not contain unique IDs") + + adata_var_merge = adata_var.copy() + + # Check for common columns between adata_var and annotation_df + # (these will be overwritten by the merge) + common_cols = annotation_df.columns[ + annotation_df.columns.isin(adata_var_merge.columns) + ].tolist() + if len(common_cols) > 0: + warnings.warn( + f"Columns {common_cols} are present in both adata_var and annotation_df - these will be overwritten", + stacklevel=2, + ) + adata_var_merge.drop(common_cols, axis=1, inplace=True) + + # Pick IDs in adata_var table + if var_on is None: + assert ( + adata_var_merge.index.is_unique + ), "var_names do not contain unique IDs - specify ID column with `var_on` parameter" + index_name = adata_var_merge.index.name + adata_var_merge = adata_var_merge.reset_index(names="var_names") + else: + if var_on not in adata_var_merge.columns: + raise ValueError( + "Column specified by `var_on` does not exist in adata_var_merge" + ) + if not adata_var_merge[var_on].is_unique: + raise ValueError("Column specified by `var_on` does not contain unique IDs") + adata_var_merge["var_names"] = adata_var_merge[var_on].copy() + + # Merge + annotated_var = ( + adata_var_merge.reset_index() + .merge( + annotation_df, how="left", right_on=[annotation_on], left_on=["var_names"] + ) + .set_index("index") + ) + + # Retain order and var_names + if var_on is None: + annotated_var = annotated_var.set_index("var_names") + annotated_var = annotated_var.loc[adata_var_merge["var_names"]] + annotated_var.index.name = index_name + else: + annotated_var = annotated_var.loc[adata_var_merge.index] + annotated_var = annotated_var.drop("var_names", axis=1) + + # Check for missing genes + missing_vars = annotated_var.index[ + annotated_var[annotation_df.drop([annotation_on], axis=1).columns] + .isna() + .all(axis=1) + ] + if len(missing_vars) > 0: + warnings.warn( + f"Missing annotations for {len(missing_vars)}/{annotated_var.shape[0]} vars", + stacklevel=2, + ) + + return annotated_var diff --git a/tests/test_annotate_anndata.py b/tests/test_annotate_anndata.py new file mode 100644 index 0000000..481e73f --- /dev/null +++ b/tests/test_annotate_anndata.py @@ -0,0 +1,92 @@ +import numpy as np +import pandas as pd +import pytest +from anndata import AnnData + +import genomic_features as gf + + +@pytest.fixture(scope="module") +def genes(): + return gf.ensembl.annotation("Hsapiens", 108).genes() + + +@pytest.fixture(scope="module") +def adata(): + return AnnData( + X=np.array([[1, 2, 3], [4, 5, 6]]), + obs=pd.DataFrame(index=["cell1", "cell2"]), + var=pd.DataFrame( + index=["ENSG00000278862", "ENSG00000212694", "ENSG00000113721"] + ), + ) + + +@pytest.fixture(scope="module") +def annotated_adata_var(adata, genes): + return gf.annotate(adata.var, genes) + + +def test_output(annotated_adata_var): + assert isinstance(annotated_adata_var, pd.DataFrame) + + +def test_size(annotated_adata_var, adata): + assert annotated_adata_var.shape[0] == adata.n_vars + + +def test_var_names(annotated_adata_var, adata): + assert "var_names" not in annotated_adata_var + assert annotated_adata_var.index.equals(adata.var.index) + assert annotated_adata_var.index.name == adata.var.index.name + + +def test_anndata_changes(adata, genes): + gf.annotate(adata.var, genes) + assert "gene_name" not in adata.var + adata.var["gene_ids"] = adata.var_names.copy() + gf.annotate(adata.var, genes, var_on="gene_ids") + assert "gene_name" not in adata.var + + +def test_missing_genes(genes): + # Test missing genes get NAs + adata_custom_genes = AnnData( + X=np.array([[1, 2, 3, 5], [4, 5, 6, 8]]), + obs=pd.DataFrame(index=["cell1", "cell2"]), + var=pd.DataFrame( + index=["ENSG00000278862", "ENSG00000212694", "ENSG00000113721", "GENE1"] + ), + ) + with pytest.warns(UserWarning): + annotated_adata_var = gf.annotate(adata_custom_genes.var, genes) + assert annotated_adata_var.loc["GENE1"][genes.columns].isna().all() + + +def test_unique_ids(genes, adata): + adata_duplicate_genes = AnnData( + X=np.array([[1, 2, 3, 5], [4, 5, 6, 8]]), + obs=pd.DataFrame(index=["cell1", "cell2"]), + var=pd.DataFrame( + index=[ + "ENSG00000278862", + "ENSG00000212694", + "ENSG00000113721", + "ENSG00000212694", + ] + ), + ) + with pytest.raises(AssertionError): + gf.annotate(adata_duplicate_genes.var, genes) + with pytest.raises(ValueError): + gf.annotate(adata.var, genes, annotation_on="gene_name") + + +def test_clashing_columns(adata, genes): + # Test that common columns are dropped + adata.var["gene_biotype"] = "foo" + adata.var["seq_name"] = "chr1" + with pytest.warns(UserWarning, match=r"gene_biotype.+seq_name"): + annotated_var = gf.annotate(adata.var, genes) + assert "lncRNA" in annotated_var["gene_biotype"].values + assert "12" in annotated_var["seq_name"].values