Repository: vanheeringen-lab/genomepy

Issue:
exactly matching NCBI assembly reports
We currently find the closest matching assembly report, but that can deviate in versions a bit (as seen with GENCODE genomes, here #183).

Discussion:
codeclimate[bot]: Code Climate has analyzed commit de57933d and detected **0 issues** on this pull request.

The test coverage on the diff in this pull request is **86.8%** (**50%** is the threshold).

This pull request will bring the total coverage in the repository to **92.4%**.

View more on [Code Climate](https://codeclimate.com/github/vanheeringen-lab/genomepy/pull/194).


Patch:
diff --git a/genomepy/annotation/mygene.py b/genomepy/annotation/mygene.py
index ac395a7..4e1c9c8 100644
--- a/genomepy/annotation/mygene.py
+++ b/genomepy/annotation/mygene.py
@@ -161,55 +161,6 @@ def _parse_mygene_input(field, product=None):
     return field, product
 
 
-# def ensembl_genome_info(self) -> Optional[Tuple[str, str, str]]:
-#     """
-#     Return Ensembl genome information for this genome.
-#
-#     Requires accession numbers to match (excluding patch numbers)
-#
-#     Returns
-#     -------
-#     (str, str, str)
-#         Ensembl name, accession, taxonomy_id
-#     """
-#     _check_property(self.readme_file, "README.txt")
-#
-#     metadata, _ = read_readme(self.readme_file)
-#     if metadata.get("provider") == "Ensembl":
-#         return metadata["name"], metadata["assembly_accession"], metadata["tax_id"]
-#
-#     if metadata.get("assembly_accession", "na") == "na":
-#         logger.warning("Cannot find a matching genome without an assembly accession.")
-#         return
-#
-#     asm_acc = metadata["assembly_accession"]
-#     search_result = nearest_assembly(asm_acc, "Ensembl")
-#     if search_result is None:
-#         return
-#
-#     return safe(search_result[0]), search_result[2], search_result[3]
-#
-#
-# def _query_mygene(
-#     self,
-#     query: Iterable[str],
-#     field: str = "genomic_pos",
-# ) -> pd.DataFrame:
-#     # mygene.info only queries the most recent version of the Ensembl database
-#     # We can only safely continue if the local genome matched the Ensembl genome.
-#     # Even if the local genome was installed via Ensembl, we still need to check
-#     # if it is the same version
-#     ensembl_info = ensembl_genome_info(self)
-#     if ensembl_info is None:
-#         return pd.DataFrame()
-#
-#     tax_id = ensembl_info[2]
-#     if not str(tax_id).isdigit():
-#         raise ValueError("No taxomoy ID found")
-#
-#     return query_mygene(query, tax_id, field)
-
-
 def _filter_query(query: pd.DataFrame) -> pd.DataFrame:
     """per queried gene, keep the best matching, non-NaN, mapping"""
     if "notfound" in query:
diff --git a/genomepy/providers/__init__.py b/genomepy/providers/__init__.py
index d6e786e..281eee8 100644
--- a/genomepy/providers/__init__.py
+++ b/genomepy/providers/__init__.py
@@ -1,16 +1,15 @@
 """Provider class, modules & related functions"""
 import os
-from typing import List, Optional
+from typing import Optional
 
 import pandas as pd
 from loguru import logger
 
 from genomepy.files import read_readme
-from genomepy.providers.base import ASM_FORMAT
 from genomepy.providers.ensembl import EnsemblProvider
 from genomepy.providers.gencode import GencodeProvider
 from genomepy.providers.local import LocalProvider
-from genomepy.providers.ncbi import NcbiProvider
+from genomepy.providers.ncbi import NcbiProvider, download_assembly_report
 from genomepy.providers.ucsc import UcscProvider
 from genomepy.providers.url import UrlProvider
 from genomepy.utils import get_genomes_dir, safe
@@ -22,7 +21,6 @@ __all__ = [
     "search",
     "map_locations",
     "download_assembly_report",
-    "nearest_assembly",
     "online_providers",
 ]
 
@@ -134,106 +132,6 @@ def search(term, provider: str = None):
             yield ret
 
 
-def _closest_patch_lvl(reference, targets):
-    ref_patch = int(reference.split(".")[1]) if "." in reference else 0
-    nearest = [999, []]
-    for target in targets:
-        tgt_patch = int(target.split(".")[1]) if "." in target else 0
-        distance = abs(tgt_patch + 0.1 - ref_patch)  # tiebreaker: newer > older patches
-        if distance == nearest[0]:
-            nearest[1].append(target)
-        if distance < nearest[0]:
-            nearest = [distance, [target]]
-    return nearest[1]
-
-
-def _best_accession(reference: str, targets: list):
-    """Return the nearest accession ID from a list of IDs"""
-    if len(targets) == 1:
-        return targets[0]
-
-    # GCA/GCF
-    matching_prefix = [t for t in targets if t.startswith(reference[0:3])]
-    if len(matching_prefix) > 0:
-        targets = matching_prefix
-
-    # patch levels
-    # e.g. GCA_000002035.4 & GCA_000002035.3
-    targets = _closest_patch_lvl(reference, targets)
-
-    if len(targets) > 1:
-        logger.info(
-            f"Multiple matching accession numbers found. Selecting the closest ({targets[0]})."
-        )
-    return targets[0]
-
-
-def _best_search_result(asm_acc: str, results: List[list]) -> list:
-    """Return the best search result based on accession IDs."""
-    results = [res for res in results if res[2] is not None]
-
-    if len(results) > 1:
-        accessions = [res[2] for res in results]
-        bes_acc = _best_accession(asm_acc, accessions)
-        results = [res for res in results if res[2] == bes_acc]
-
-    if len(results) > 0:
-        return results[0]
-
-
-def nearest_assembly(asm_acc: str, provider: str) -> list:
-    """
-    Return the search result of the assembly nearest to
-    the given accession ID, in the specified provider.
-    """
-    if not asm_acc.startswith(("GCA_", "GCF_")):
-        raise ValueError("asm_acc must be an assembly accession ID.")
-
-    search_results = list(search(asm_acc, provider=provider))
-    best_result = _best_search_result(asm_acc, search_results)
-    if best_result is None:
-        logger.warning(f"No assembly similar to {asm_acc} on {provider}.")
-    return best_result
-
-
-def download_assembly_report(asm_acc: str, fname: str = None):
-    """
-    Retrieve the NCBI assembly report.
-
-    Returns the assembly_report as a pandas DataFrame if fname is not specified.
-
-    Parameters
-    ----------
-    asm_acc : str
-        Assembly accession (GCA or GCF)
-    fname : str, optional
-        Save assembly_report to this filename.
-
-    Returns
-    -------
-    pandas.DataFrame
-        NCBI assembly report.
-    """
-    search_result = nearest_assembly(asm_acc, "NCBI")
-    if search_result is None:
-        return
-    ncbi_acc = search_result[2]
-    ncbi_name = search_result[0]
-
-    # NCBI FTP location of assembly report
-    assembly_report = (
-        f"https://ftp.ncbi.nlm.nih.gov/genomes/all/{ncbi_acc[0:3]}/"
-        + f"{ncbi_acc[4:7]}/{ncbi_acc[7:10]}/{ncbi_acc[10:13]}/"
-        + f"{ncbi_acc}_{ncbi_name}/{ncbi_acc}_{ncbi_name}_assembly_report.txt"
-    )
-    asm_report = pd.read_csv(assembly_report, sep="\t", comment="#", names=ASM_FORMAT)
-
-    if fname:
-        asm_report.to_csv(fname, sep="\t", index=False)
-    else:
-        return asm_report
-
-
 def map_locations(
     frm: str, to: str, genomes_dir: Optional[str] = None
 ) -> Optional[pd.DataFrame]:
diff --git a/genomepy/providers/base.py b/genomepy/providers/base.py
index a9d30c0..ad85b27 100644
--- a/genomepy/providers/base.py
+++ b/genomepy/providers/base.py
@@ -18,19 +18,6 @@ from genomepy.files import extract_archive, get_file_info, update_readme
 from genomepy.online import download_file
 from genomepy.utils import get_genomes_dir, get_localname, lower, mkdir_p, rm_rf, safe
 
-ASM_FORMAT = [
-    "Sequence-Name",
-    "Sequence-Role",
-    "Assigned-Molecule",
-    "Assigned-Molecule-Location/Type",
-    "GenBank-Accn",
-    "Relationship",
-    "RefSeq-Accn",
-    "Assembly-Unit",
-    "Sequence-Length",
-    "UCSC-style-name",
-]
-
 
 class BaseProvider:
     """
@@ -111,7 +98,7 @@ class BaseProvider:
             if tid.isdigit():
                 return int(tid)
 
-    def assembly_accession(self, name: str) -> str:
+    def assembly_accession(self, name: str) -> str or None:
         """
         Return the assembly accession number (GCA* or GCF*) for a genome.
 
diff --git a/genomepy/providers/gencode.py b/genomepy/providers/gencode.py
index b8668fd..abd7b42 100644
--- a/genomepy/providers/gencode.py
+++ b/genomepy/providers/gencode.py
@@ -1,16 +1,16 @@
 import os
 from time import sleep
 
-import pandas as pd
 from loguru import logger
 
 from genomepy.caching import cache
 from genomepy.exceptions import GenomeDownloadError
 from genomepy.files import update_readme
 from genomepy.online import check_url, connect_ftp_link
-from genomepy.providers.base import ASM_FORMAT, BaseProvider, download_annotation
+from genomepy.providers.base import BaseProvider, download_annotation
+from genomepy.providers.ncbi import download_assembly_report
 from genomepy.providers.ucsc import UcscProvider
-from genomepy.utils import get_genomes_dir, get_localname, mkdir_p
+from genomepy.utils import get_genomes_dir, get_localname
 
 
 class GencodeProvider(BaseProvider):
@@ -147,7 +147,7 @@ class GencodeProvider(BaseProvider):
             # download exact assembly report to rename the scaffolds
             acc = self.assembly_accession(name)
             fname = os.path.join(genomes_dir, localname, "assembly_report.txt")
-            exact_assembly_report(name, acc, fname)
+            download_assembly_report(acc, fname)
 
             download_annotation(genomes_dir, link, localname)
             logger.info("Annotation download successful")
@@ -286,14 +286,3 @@ def get_genomes(ftp_link):
 
     ftp.quit()
     return genomes
-
-
-def exact_assembly_report(name, acc, fname):
-    assembly_report = (
-        f"https://ftp.ncbi.nlm.nih.gov/genomes/all/{acc[0:3]}/"
-        + f"{acc[4:7]}/{acc[7:10]}/{acc[10:13]}/"
-        + f"{acc}_{name}/{acc}_{name}_assembly_report.txt"
-    )
-    asm_report = pd.read_csv(assembly_report, sep="\t", comment="#", names=ASM_FORMAT)
-    mkdir_p(os.path.dirname(fname))
-    asm_report.to_csv(fname, sep="\t", index=False)
diff --git a/genomepy/providers/ncbi.py b/genomepy/providers/ncbi.py
index a0e2fc4..2d7d13e 100644
--- a/genomepy/providers/ncbi.py
+++ b/genomepy/providers/ncbi.py
@@ -1,5 +1,6 @@
 import os
 import re
+import urllib.error
 from typing import List
 from urllib.request import urlopen
 
@@ -9,9 +10,22 @@ from tqdm.auto import tqdm
 
 from genomepy.caching import cache
 from genomepy.exceptions import GenomeDownloadError
-from genomepy.online import check_url
-from genomepy.providers.base import BaseProvider, ASM_FORMAT
-from genomepy.utils import safe
+from genomepy.online import check_url, read_url
+from genomepy.providers.base import BaseProvider
+from genomepy.utils import mkdir_p, safe
+
+ASM_FORMAT = [
+    "Sequence-Name",
+    "Sequence-Role",
+    "Assigned-Molecule",
+    "Assigned-Molecule-Location/Type",
+    "GenBank-Accn",
+    "Relationship",
+    "RefSeq-Accn",
+    "Assembly-Unit",
+    "Sequence-Length",
+    "UCSC-style-name",
+]
 
 
 class NcbiProvider(BaseProvider):
@@ -104,13 +118,10 @@ class NcbiProvider(BaseProvider):
         """
         # (down)load the assembly report
         asm_fname = os.path.join(out_dir, "assembly_report.txt")
-        if os.path.exists(asm_fname):
-            asm_report = pd.read_csv(asm_fname, sep="\t", comment="#", header=None)
-        else:
-            url = self._ftp_or_html_link(name, "_assembly_report.txt", True)
-            asm_report = pd.read_csv(url, sep="\t", comment="#", header=ASM_FORMAT)
+        if not os.path.exists(asm_fname):
             # save assembly report
-            asm_report.to_csv(asm_fname, sep="\t", index=False)
+            download_assembly_report(self.assembly_accession(name), asm_fname)
+        asm_report = pd.read_csv(asm_fname, sep="\t", comment="#", header=None)
 
         # create mapping of chromosome accessions to chromosome names
         tr = asm_report.set_index(6)[0].to_dict()
@@ -216,3 +227,80 @@ def get_genomes(assembly_url):
                 name = safe(line[15])  # overwrites older asm_names
                 genomes[name] = dict(zip(header, line))
     return genomes
+
+
+def download_assembly_report(acc: str, fname: str = None):
+    """
+    Retrieve the NCBI assembly report.
+
+    Returns the assembly_report as a pandas DataFrame if fname is not specified.
+
+    Parameters
+    ----------
+    acc : str
+        Assembly accession (GCA or GCF)
+    fname : str, optional
+        Save assembly_report to this filename.
+
+    Returns
+    -------
+    pandas.DataFrame
+        NCBI assembly report.
+    """
+    msg = "Could not download the assembly report from NCBI. "
+    if not isinstance(acc, str) or not acc.startswith(("GCA", "GCF")):
+        logger.warning(msg)
+        return None
+    assembly_report = _assembly_report_url(acc)
+    if assembly_report is None:
+        logger.warning(msg + f"Assembly accession '{acc}' not found.")
+        return None
+    asm_report = pd.read_csv(assembly_report, sep="\t", comment="#", names=ASM_FORMAT)
+
+    if fname:
+        mkdir_p(os.path.dirname(fname))
+        asm_report.to_csv(fname, sep="\t", index=False)
+    else:
+        return asm_report
+
+
+def _assembly_report_url(acc: str) -> str or None:
+    ftp_basedir = "https://ftp.ncbi.nlm.nih.gov/genomes/all"
+    asm_dir = f"{ftp_basedir}/{acc[0:3]}/{acc[4:7]}/{acc[7:10]}/{acc[10:13]}/"
+    try:
+        text = read_url(asm_dir)
+    except urllib.error.HTTPError:
+        return None
+
+    # try to find exact accession patch level
+    hits = re.findall(acc + '_.*/"', text)  # len 0-1
+    if len(hits) == 0:
+        # patch level not found, pick the closest instead
+        hits = re.findall(acc + '.*/"', text)
+        available_accessions = ["_".join(h.split("_")[0:2]) for h in hits]
+        closest = _closest_patch_lvl(acc, available_accessions)
+        hits = [h for h in hits if closest + "_" in h]  # len 0-1
+    if len(hits) == 0:
+        return None
+
+    name = hits[0][:-2]
+    assembly_report = asm_dir + name + "/" + name + "_assembly_report.txt"
+    return assembly_report
+
+
+def _closest_patch_lvl(reference: str, targets: list) -> str:
+    ref_patch = _patch_lvl(reference)
+    closest = ""
+    closest_dist = 999
+    for target in targets:
+        tgt_patch = _patch_lvl(target)
+        dist = abs(tgt_patch - 0.1 - ref_patch)  # tiebreaker: newer > older patches
+        if dist < closest_dist:
+            closest = target
+            closest_dist = dist
+    return closest
+
+
+def _patch_lvl(acc: str) -> int:
+    """str(GCA_000000000.6) -> int(6)"""
+    return int(acc.split(".")[1]) if "." in acc else 0
diff --git a/genomepy/providers/ucsc.py b/genomepy/providers/ucsc.py
index c4eadcf..3305cbd 100644
--- a/genomepy/providers/ucsc.py
+++ b/genomepy/providers/ucsc.py
@@ -123,7 +123,7 @@ class UcscProvider(BaseProvider):
                     yield name
                     break  # max one hit per genome
 
-    def assembly_accession(self, name: str) -> str:
+    def assembly_accession(self, name: str) -> str or None:
         """
         Return the assembly accession (`GCA_`/`GCF_`) for a genome.
 
@@ -149,8 +149,6 @@ class UcscProvider(BaseProvider):
         if acc:
             return acc
 
-        return "na"
-
     def annotation_links(self, name, **kwargs) -> List[str]:
         """
         Return a sorted list of available gene annotation types for a genome
@@ -619,6 +617,7 @@ def scrape_accession(htmlpath: str) -> str:
         # contains additional info, such as '(latest)' or '(suppressed)'. Unused for now.
         valid_accessions = re.findall(r"assembly accession:.*?GC[AF]_.*?<", text)
         text = " ".join(valid_accessions)
+        # this selects GCA > GCF if both are found
         match = accession_regex.search(text)
         if match:
             return match.group(0)


Test patch:
diff --git a/tests/test_09_provider_ncbi.py b/tests/test_09_provider_ncbi.py
index 6f6291e..2142655 100644
--- a/tests/test_09_provider_ncbi.py
+++ b/tests/test_09_provider_ncbi.py
@@ -2,6 +2,10 @@ import os
 from shutil import copyfile
 from tempfile import TemporaryDirectory
 
+import pandas as pd
+
+import genomepy
+
 
 def test_ncbiprovider(ncbi):
     assert ncbi.name == "NCBI"
@@ -64,3 +68,12 @@ def test__get_genomes(ncbi):
         assert field in genome
     assert genome["species_taxid"] == "2097"
     assert genome["taxid"] == "243273"
+
+
+def test_download_assembly_report():
+    assembly_report = "tests/data/sacCer3/assembly_report.txt"
+    genomepy.providers.download_assembly_report("GCA_000146045", assembly_report)
+    report = pd.read_csv(assembly_report, sep="\t", comment="#")
+
+    assert isinstance(report, pd.DataFrame)
+    assert list(report.columns) == genomepy.providers.ncbi.ASM_FORMAT
diff --git a/tests/test_11_provider.py b/tests/test_11_provider.py
index 2d937e7..43d75ac 100644
--- a/tests/test_11_provider.py
+++ b/tests/test_11_provider.py
@@ -60,169 +60,6 @@ def test_search():
     assert "8364" in str(metadata[3])
 
 
-def test__best_search_result():
-    # ncbi_human = list(genomepy.search("GCA_000001405", provider="NCBI"))
-    ncbi_human = [
-        [
-            "GRCh38",
-            "NCBI",
-            "GCF_000001405.26",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p1",
-            "NCBI",
-            "GCF_000001405.27",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p2",
-            "NCBI",
-            "GCF_000001405.28",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p3",
-            "NCBI",
-            "GCF_000001405.29",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p4",
-            "NCBI",
-            "GCF_000001405.30",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p5",
-            "NCBI",
-            "GCF_000001405.31",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p6",
-            "NCBI",
-            "GCF_000001405.32",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p7",
-            "NCBI",
-            "GCF_000001405.33",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p8",
-            "NCBI",
-            "GCF_000001405.34",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p9",
-            "NCBI",
-            "GCF_000001405.35",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p10",
-            "NCBI",
-            "GCF_000001405.36",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p11",
-            "NCBI",
-            "GCF_000001405.37",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p12",
-            "NCBI",
-            "GCF_000001405.38",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh37.p13",
-            "NCBI",
-            "GCF_000001405.25",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-        [
-            "GRCh38.p13",
-            "NCBI",
-            "GCF_000001405.39",
-            9606,
-            True,
-            "Homo sapiens",
-            "Genome Reference Consortium",
-        ],
-    ]
-
-    # list of search results
-    best_result = genomepy.providers._best_search_result("GCA_000001405.39", ncbi_human)
-    assert best_result[2] == "GCF_000001405.39"
-
-    # zero search results
-    assert genomepy.providers._best_search_result("GCA_000001405.39", []) is None
-
-    # one search results
-    best_result = genomepy.providers._best_search_result(
-        "GCA_000001405.26", [ncbi_human[0]]
-    )
-    assert best_result[2] == ncbi_human[0][2]
-
-
-def test_download_assembly_report():
-    assembly_report = "tests/data/sacCer3/assembly_report.txt"
-    genomepy.providers.download_assembly_report("GCA_000146045", assembly_report)
-    report = pd.read_csv(assembly_report, sep="\t", comment="#")
-
-    assert isinstance(report, pd.DataFrame)
-    assert list(report.columns) == genomepy.providers.ASM_FORMAT
-
-
 def test_map_location():
     mapping = genomepy.providers.map_locations(
         frm="sacCer3", to="ensembl", genomes_dir="tests/data"
diff --git a/tests/test_13_annotation.py b/tests/test_13_annotation.py
index ae1b6f7..774fafb 100644
--- a/tests/test_13_annotation.py
+++ b/tests/test_13_annotation.py
@@ -47,7 +47,7 @@ def test__get_name_and_dir():
         assert genome_dir == genomepy.utils.cleanpath(td["expected_dir"]), td
 
     # does not exist
-    with pytest.raises(ValueError):
+    with pytest.raises(FileNotFoundError):
         genomepy.annotation._get_name_and_dir("tests/data/GRCz11/what.who")
     # no bed/gtf/fa
     with pytest.raises(NotImplementedError):
@@ -83,7 +83,7 @@ def test_annotation_init(caplog, annot):
     assert "Could not find 'empty.annotation.gtf" in caplog.text
 
     # Genome doesn't exist
-    with pytest.raises(ValueError):
+    with pytest.raises(FileNotFoundError):
         genomepy.Annotation("never_existed", genomes_dir="tests/data")
 
 
diff --git a/tests/test_14_functions.py b/tests/test_14_functions.py
index 0583319..f7778c8 100644
--- a/tests/test_14_functions.py
+++ b/tests/test_14_functions.py
@@ -20,7 +20,7 @@ def test_list_available_genomes():
     g = genomepy.functions.list_available_genomes("Ensembl")
     metadata = next(g)
     assert isinstance(metadata, list)
-    assert metadata[0:2] == ["athCun1", "Ensembl"]
+    assert metadata[0:2] == ["mCalJac1.pat.X", "Ensembl"]
 
 
 def test_list_installed_genomes():