diff --git a/.gitignore b/.gitignore index 3db890198..d2c7af591 100644 --- a/.gitignore +++ b/.gitignore @@ -43,6 +43,7 @@ __pycache__ .pypirc .gitignore +.python-version chauthorinfo.sh pre-commit.sh @@ -55,3 +56,5 @@ Documentation/* *log *BAK *.sublime* +*.code-workspace +.vscode/ diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 000000000..a7484d773 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,19 @@ +FROM python:3.10.5-slim-buster + +WORKDIR /app + +COPY ./scripts/motif_search/ . + +COPY ./prody ./prody + +RUN apt update -y && \ + apt install build-essential -y + +RUN pip install --upgrade pip setuptools wheel + +RUN pip install -r requirements.txt + +RUN rm requirements.txt + +CMD ["flask", "run", "--host", "0.0.0.0", "--port", "5000"] + diff --git a/prody/database/__init__.py b/prody/database/__init__.py index 4091df195..5178ec852 100644 --- a/prody/database/__init__.py +++ b/prody/database/__init__.py @@ -69,6 +69,33 @@ .. _GOA: https://www.ebi.ac.uk/GOA/ +Swiss-Prot +================================ +The following classes and functions can be used to search and retrieve data from the Swiss-Prot database: + * :class:`.SwissProt` - class to handle Swiss-Prot data from Expasy + * :func:`.getCurrentRelease` - gets current Swiss-Prot release version + * :func:`.downloadRelease` - downloads current Swiss-Prot database files + * :func:`.saveRelease` - saves new Swiss-Prot release version + * :func:`.updateRelease` - updates Swiss-Prot local database + * :func:`.getLocalRelease` - checks local Swiss-Prot release version + * :func:`.checkForUpdates` - checks wheather there is newer Swiss-Prot version than current local one + + RefSeq + ================================ +The following classes and functions can be used to search and retrieve data from the RefSeq database: + * :class:` .RefSeq` - class to handle RefSeq data + * :func:` getCurrentRelease` - func desc + * :func:` getInstalledFiles` - func desc + * :func:` saveInstalledFiles` - func desc + * :func:` getLocalFiles` - func desc + * :func:` getFiles` - func desc + * :func:` pepareDownloadFileList` - func desc + * :func:` downloadRelease` - func desc + * :func:` saveRelease` - func desc + * :func:` updateRelease` - func desc + * :func:` getLocalRelease` - func desc + * :func:` checkForUpdates` - func desc + Interpro ==== @@ -115,6 +142,18 @@ from .quartataweb import * __all__.extend(quartataweb.__all__) +from . import swissprot +from .swissprot import * +__all__.extend(swissprot.__all__) + +from . import refseq +from .refseq import * +__all__.extend(refseq.__all__) + +from . import pdb +from .pdb import * +__all__.extend(pdb.__all__) + from . import interpro from .interpro import * __all__.extend(interpro.__all__) diff --git a/prody/database/pdb.py b/prody/database/pdb.py new file mode 100644 index 000000000..10ba9bfea --- /dev/null +++ b/prody/database/pdb.py @@ -0,0 +1,28 @@ +# -*- coding: utf-8 -*- + +"""Protein Data Bank database operations.""" + +import os + +from prody import LOGGER +from prody.utilities.motifhelpers import downloadFile + +__all__ = ["ProteinDataBank"] + + +class ProteinDataBank: + """Protein Data Bank""" + + @classmethod + def downloadRelease(cls, data_dir: str) -> None: + """Download latest pdb database. + + Args: + data_dir (str): Directory to download to. + """ + url = "https://ftp.wwpdb.org/pub/pdb/derived_data/" + file = "pdb_seqres.txt.gz" + LOGGER.timeit("downloadRelease") + directory = os.path.join(data_dir, cls.__name__) + downloadFile(url, directory, file) + LOGGER.report(msg="downloadRelease completed in %.2fs.", label="downloadRelease") diff --git a/prody/database/refseq.py b/prody/database/refseq.py new file mode 100644 index 000000000..52d7a0742 --- /dev/null +++ b/prody/database/refseq.py @@ -0,0 +1,180 @@ +# -*- coding: utf-8 -*- + +"""RefSeq database operations.""" + +import os +import re +from concurrent.futures import ThreadPoolExecutor + +import requests +import requests_cache +from datetime import timedelta +from prody import LOGGER +from prody.utilities.motifhelpers import downloadFile + +__all__ = ["RefSeq"] + +requests_cache.install_cache(expire_after=timedelta(weeks=1)) + + +class RefSeq: + """RefSeq database.""" + + RELEASE_FILE = "release.txt" + CHECK_SUMS = "files_installed.txt" + + @classmethod + def getCurrentRelease(cls) -> str: + """Get current RefSeq db release version. + + Raises: + Exception: No release version found. + + Returns: + str: RefSeq db release version. + """ + rs_release = "" + url = "https://ftp.ncbi.nlm.nih.gov/refseq/release/RELEASE_NUMBER" + try: + with requests_cache.disabled(): + response = requests.get(url) + except requests.exceptions.RequestException as exception: + LOGGER.error(str(exception)) + else: + rs_release = response.text.strip() + if not re.match(r"\d+", rs_release): + LOGGER.error("Could't determine release version.") + LOGGER.debug("RefSeq current release: {}".format(rs_release)) + return rs_release + + @classmethod + def getInstalledFiles(cls) -> str: + """Reads installed files with corresponding check sums.""" + LOGGER.debug("Downloading installed files with check sums.") + current_release = cls.getCurrentRelease() + url = "https://ftp.ncbi.nlm.nih.gov/refseq/release/release-catalog/release{}.files.installed".format(current_release) + try: + with requests_cache.disabled(): + response = requests.get(url) + except requests.exceptions.RequestException as exception: + LOGGER.error(str(exception)) + else: + result = "" + lines = response.text.split("\n") + for line in lines: + if re.search(r"complete\.(\d+\.){1,2}protein\.faa\.gz", line): + result += line + "\n" + return result + + @classmethod + def saveInstalledFiles(cls, data_dir: str) -> None: + """Saves installed files list with check sums locally.""" + installed_files = cls.getInstalledFiles() + path = os.path.join(data_dir, cls.__name__) + os.makedirs(path, exist_ok=True) + with open("{}/{}".format(path, cls.CHECK_SUMS), "w", encoding="utf-8") as file: + file.write(installed_files) + LOGGER.debug("{} file saved locally.".format(cls.CHECK_SUMS)) + + @classmethod + def getLocalFiles(cls, data_dir: str) -> dict: + """Lists local RefSeq FASTA protein files and corresponding check sums. + + Returns: + dict: file names and corresponding check sums. + """ + LOGGER.debug("Getting local RefSeq protein FASTA files") + path = os.path.join(data_dir, cls.__name__) + os.makedirs(path, exist_ok=True) + filename = os.path.join(path, cls.CHECK_SUMS) + if os.path.exists(filename): + with open(filename, "r", encoding="utf-8") as file: + text = file.read() + results = re.findall(r"^(\d+)\s+(\S+)$", text, re.M) + return {result[1]: result[0] for result in results} + else: + return {} + + @classmethod + def getFiles(cls) -> dict: + """Lists all FASTA protein files on RefSeq ftp server. + + Returns: + dict: FASTA protein files with check sums. + """ + LOGGER.debug("Getting protein FASTA file list from RefSeq.") + installed_files = cls.getInstalledFiles() + file_matcher = re.compile(r"^(\d+)\s+(\S+)$", re.M) + files = re.findall(file_matcher, installed_files) + return {file[1]: file[0] for file in files} if files else {} + + @classmethod + def pepareDownloadFileList(cls, data_dir: str) -> list: + """Prepare file list to be downloaded""" + LOGGER.debug("Preparing file list to be downloaded.") + remote_files = cls.getFiles() + local_files = cls.getLocalFiles(data_dir) + download_list = [] + for filename in remote_files.keys(): + if filename in local_files.keys(): + if remote_files[filename] != local_files[filename]: + download_list.append(filename) + else: + download_list.append(filename) + return download_list + + @classmethod + def downloadRelease(cls, data_dir: str) -> None: + """Download latest RefSeq database release.""" + + def report_done(future): + LOGGER.info("Downloaded file {}".format(future.result())) + + files = cls.pepareDownloadFileList(data_dir) + url = "https://ftp.ncbi.nlm.nih.gov/refseq/release/complete/" + LOGGER.timeit("downloadRelease") + with ThreadPoolExecutor(max_workers=100) as executor: + for file in files: + directory = os.path.join(data_dir, cls.__name__) + future = executor.submit(downloadFile, url, directory, file) + future.add_done_callback(report_done) + LOGGER.report(msg="downloadRelease completed in %.2fs.", label="downloadRelease") + + @classmethod + def saveRelease(cls, data_dir: str) -> None: + """Write current release version to disk.""" + current_release = cls.getCurrentRelease() + path = os.path.join(data_dir, cls.__name__) + os.makedirs(path, exist_ok=True) + with open("{}/{}".format(path, cls.RELEASE_FILE), "w", encoding="utf-8") as file: + file.write(current_release) + LOGGER.debug("RefSeq release {} saved.".format(current_release)) + + @classmethod + def updateRelease(cls, data_dir) -> None: + """Update release to the most recent one.""" + LOGGER.debug("Updating RefSeq release version.") + cls.downloadRelease(data_dir) + cls.saveRelease(data_dir) + + @classmethod + def getLocalRelease(cls, data_dir) -> str: + """Get release version from local disk.""" + path = os.path.join(data_dir, cls.__name__) + os.makedirs(path, exist_ok=True) + version = "" + release = os.path.join(path, cls.RELEASE_FILE) + if os.path.exists(release): + with open(release, "r", encoding="utf-8") as file: + version = file.readline() + LOGGER.debug("RefSeq local release: {}".format(version)) + else: + LOGGER.debug("RefSeq local release not found.") + return version + + @classmethod + def checkForUpdates(cls, data_dir) -> None: + """Check if local version is the recent one.""" + LOGGER.debug("RefSeq checking for updates.") + if cls.getCurrentRelease() != cls.getLocalRelease(data_dir): + cls.updateRelease(data_dir) diff --git a/prody/database/swissprot.py b/prody/database/swissprot.py new file mode 100644 index 000000000..a97c6ae4f --- /dev/null +++ b/prody/database/swissprot.py @@ -0,0 +1,110 @@ +# -*- coding: utf-8 -*- + +"""Swiss-Prot database operations.""" + +import os +import re +import requests +import requests_cache +from datetime import timedelta +from concurrent.futures import ThreadPoolExecutor + +from prody import LOGGER +from prody.utilities.motifhelpers import downloadFile + +__all__ = ["SwissProt"] + +requests_cache.install_cache(expire_after=timedelta(weeks=1)) + + +class SwissProt: + """Swiss-Prot database.""" + + RELEASE_FILE = "release.txt" + HEADERS = {"User-Agent": "Python pattern search agent", "Contact": "mkonstanty@gmail.com"} + + @classmethod + def getCurrentRelease(cls) -> str: + """Get current swiss-prot db release version. + + Raises: + Exception: No release version found. + + Returns: + str: Swiss-prot db release version. + """ + sp_release = "" + url = "https://ftp.expasy.org/databases/swiss-prot/release/reldate.txt" + try: + with requests_cache.disabled(): + response = requests.get(url, headers=cls.HEADERS) + except requests.exceptions.RequestException as exception: + LOGGER.error(str(exception)) + else: + sp_release = re.search(r"Swiss-Prot Release (\d{4}_\d{2})", response.text) + if not sp_release: + LOGGER.error("Could't determine release version.") + LOGGER.debug("Swiss-Prot current release: {}".format(sp_release[1])) + return sp_release[1] + + @classmethod + def downloadRelease(cls, data_dir: str, types=None) -> None: + """Download latest swiss-prot database release. + + Args: + types (list, optional): Database file types. Defaults to None. + """ + + def report_done(future): + LOGGER.info("Downloaded file {}".format(future.result())) + + types = types if types else ["xml", "dat", "fasta"] + files = ["uniprot_sprot.{}.gz".format(type) for type in types] + url = "https://ftp.expasy.org/databases/swiss-prot/release/" + LOGGER.timeit("downloadRelease") + with ThreadPoolExecutor(max_workers=3) as executor: + for file in files: + directory = os.path.join(data_dir, cls.__name__) + future = executor.submit(downloadFile, url, directory, file, headers=cls.HEADERS) + future.add_done_callback(report_done) + LOGGER.report(msg="downloadRelease completed in %.2fs.", label="downloadRelease") + + @classmethod + def saveRelease(cls, data_dir: str) -> None: + """Write current release version to disk.""" + current_release = cls.getCurrentRelease() + path = os.path.join(data_dir, cls.__name__) + os.makedirs(path, exist_ok=True) + with open("{}/{}".format(path, cls.RELEASE_FILE), "w", encoding="utf-8") as file: + file.write(current_release) + LOGGER.debug("Swiss-Prot release {} saved.".format(current_release)) + + @classmethod + def updateRelease(cls, data_dir) -> None: + """Update release to the most recent one.""" + LOGGER.debug("Updating Swiss-Prot release version.") + cls.downloadRelease(data_dir, types=["fasta"]) + cls.saveRelease(data_dir) + + @classmethod + def getLocalRelease(cls, data_dir: str) -> str: + """Get release version from local disk.""" + LOGGER.debug("Getting Swiss-Prot local release version.") + path = os.path.join(data_dir, cls.__name__) + os.makedirs(path, exist_ok=True) + version = "" + release = os.path.join(path, cls.RELEASE_FILE) + if os.path.exists(release): + with open(release, "r", encoding="utf-8") as file: + version = file.readline() + LOGGER.debug("Swiss-Prot local release: {}".format(version)) + else: + LOGGER.debug("Swiss-Prot local release not found.") + return version + + @classmethod + def checkForUpdates(cls, data_dir) -> None: + """Check if local version is the recent one.""" + LOGGER.debug("Swiss-Prot checking for updates.") + if cls.getCurrentRelease() != cls.getLocalRelease(data_dir): + cls.updateRelease(data_dir) diff --git a/prody/sequence/__init__.py b/prody/sequence/__init__.py index 9545bb480..4243de576 100644 --- a/prody/sequence/__init__.py +++ b/prody/sequence/__init__.py @@ -46,10 +46,19 @@ * :func:`.showShannonEntropy` - plot Shannon entropy * :func:`.showMSAOccupancy` - plot row (sequence) or column occupancy * :func:`.showMutinfoMatrix` - show mutual information matrix + + +Searching +======== + * :func:`.getPdbCodesFromMotif` - get PDB code from MOTIF """ __all__ = [] +from . import motif +from .motif import * +__all__.extend(motif.__all__) + from . import msa from .msa import * __all__.extend(msa.__all__) diff --git a/prody/sequence/motif.py b/prody/sequence/motif.py new file mode 100644 index 000000000..444c55d4f --- /dev/null +++ b/prody/sequence/motif.py @@ -0,0 +1,491 @@ +# -*- coding: utf-8 -*- + +"""This module prase protein sequences using PROSITE motifs.""" + +import os +import re +import requests +import requests_cache +from datetime import timedelta +import json +from typing import Optional + +from prody import LOGGER +from prody.utilities.motifhelpers import prositeToRegEx, openFile + + +__all__ = [ + "getUniprot", + "getPdbCodesFromMotif", + "saveMotifResults", + "expasySearchMotif", + "spOnlineMotifSearch", + "pdbOnlineMotifSearch", + "localMotifSearch", + "spLocalMotifSearch", + "pdbLocalMotifSearch", + "motifSearch", +] + +MOTIF_PATTERN = ( + r"^([ARDNCEQGHILKMFPSTWYVx]*" + r"(\[[ARDNCEQGHILKMFPSTWYV]+\])*" + r"({{[{{ARDNCEQGHILKMFPSTWYV}}]+}})*" + r"([ARDNCEQGHILKMFPSTWYVx]\(\d(,\d)?\))*)+$" +) +MOTIF_MATCHER = re.compile(MOTIF_PATTERN) + +requests_cache.install_cache(expire_after=timedelta(weeks=1)) + + +def getUniprot(id: str) -> str: + """Get details for accession key from UniProt + + Args: + id (str): Uniprot accession key + + Returns: + str: details + """ + LOGGER.debug("Get {} details from Uniprot.".format(id)) + headers = {"User-Agent": "Python pattern search agent", "Contact": "mkonstanty@gmail.com"} + url = "https://www.uniprot.org/uniprot/{}.txt".format(id) + try: + response = requests.get(url, headers=headers) + except requests.exceptions.RequestException as exception: + LOGGER.error(str(exception)) + result = response.text.replace("\n", "<\\br>") + return result + + +def getPdbCodesFromMotif(motif: str, results: Optional[int] = None, page: Optional[int] = None) -> list: + """Get PDB code from motif. + + Args: + motif (str): FASTA motif + results (Optional[int], optional): How many results to return. Defaults to None. + page (Optional[int], optional): Which result page to return. Defaults to None. + + Returns: + list: PDB ids + """ + LOGGER.debug("Get PDB code from motif {}".format(motif)) + url = "https://search.rcsb.org/rcsbsearch/v2/query" + headers = {"Content-Type": "application/json;charset=utf-8"} + if results and page: + paginate = {"paginate": {"start": results * (page - 1), "rows": results * page - 1}} + else: + paginate = {"return_all_hits": True} + data = { + "return_type": "entry", + "query": { + "type": "terminal", + "service": "seqmotif", + "parameters": {"value": motif, "pattern_type": "prosite", "target": "pdb_protein_sequence"}, + }, + "request_options": paginate, + } + try: + LOGGER.timeit("getPdbCodesFromMotif") + response = requests.post(url, headers=headers, data=json.dumps(data)) + LOGGER.debug("Result taken from cache: {}".format(response.from_cache)) + LOGGER.report(msg="getPdbCodesFromMotif completed in %.2fs.", label="getPdbCodesFromMotif") + except requests.exceptions.RequestException as exception: + LOGGER.error(str(exception)) + else: + response.encoding = "UTF-8" + result = json.loads(response.text) + result["result_set"] + return [x["identifier"] for x in result["result_set"]] + + +def saveMotifResults(database: str, motif: str, data: list[dict], directory: str) -> str: + """Save search data to a file. + + Args: + database (str): name of the database + motif (str): FASTA motif + data (list[dict]): List of proteins to be saved. + directory (str): directory to save results + + Returns: + str: saved file name + """ + os.makedirs(directory, exist_ok=True) + filename = "{}_{}.txt".format(motif.replace('-', ''), database) + with open("{}/{}".format(directory, filename), "w", encoding="utf-8") as file: + for protein in data: + for key, value in protein.items(): + file.write("{}: {}\n".format(key, value)) + file.write("\n") + LOGGER.info("Search written to a file {}.".format(filename)) + return filename + + +def _expasySPMotifToStructure(response: str) -> list: + """Parse HTML response from Expasy. + + Args: + response (str): HTML response from Expasy + + Returns: + list: List of result dicts + """ + LOGGER.timeit("_expasySPMotifToStructure") + RESULT_PATTERN = ( + r">sp\|" # database + r"(\w+)\|(\w+).*\n" # accession and id + r"(.*)\n" # description + r"([ARDNCEQGHILKMFPSTWYVx\n]+)\W+" # sequence + r"(\d{,3}) - (\d{,3}):\W+(\w+)" # start, end and match + ) + motifs = re.compile(RESULT_PATTERN, re.M) + results = re.findall(motifs, response) + structure = [] + keys = ("accession", "id", "description", "sequence", "start", "end", "match") + for result in results: + structure.append(dict(zip(keys, result))) + LOGGER.report(msg="_expasySPMotifToStructure completed in %.2fs.", label="_expasySPMotifToStructure") + return structure + + +def _expasyPDBMotifToStructure(response: str) -> list: + """Parse HTML response from Expasy. + + Args: + response (str): HTML response from Expasy + + Returns: + list: List of result dicts + """ + LOGGER.timeit("_expasyPDBMotifToStructure") + structure = [] + RESULT_PATTERN = ( + r">pdb\|(\w+).*\n" # pdb id + r".*?((\w+\|\w+[,.])+)\n" # accession and uniprot id + r"\W+(\d{,3}) - (\d{,3}):\W+(\w+)" # start, end and match + ) + motifs = re.compile(RESULT_PATTERN, re.M) + results = re.findall(motifs, response) + for result in results: + pdb_id, tmp, _, start, end, match = result + accessions_ids = dict(re.findall(r"(\w+)\|(\w+)", tmp)) + accessions = ", ".join(accessions_ids.keys()) + ids = ", ".join(accessions_ids.values()) + structure.append( + { + "PDB id": pdb_id, + "Accession ids": accessions, + "UniProt ids": ids, + "start": start, + "end": end, + "match": match, + } + ) + LOGGER.report(msg="_expasyPDBMotifToStructure completed in %.2fs.", label="_expasyPDBMotifToStructure") + return structure + + +def expasySearchMotif(motif: str, database: str) -> list: + """Search motif in remote Swiss-Prot database. + + https://prosite.expasy.org/scanprosite/scanprosite_doc.html#rest_intro + + Args: + motif (str): Motif in PROSITE format. + database (str): 'sp' (UniProtKB/Swiss-Prot) or 'tr' (UniProtKB/TrEMBL) or 'pdb' (PDB). + + Returns: + list: List of result dicts + """ + + def err(*args) -> list: + return [] + + LOGGER.debug("Searching motif {} in {} database via Expasy".format(motif, database)) + headers = {"User-Agent": "Python pattern search agent", "Contact": "mkonstanty@gmail.com"} + url = "https://prosite.expasy.org/cgi-bin/prosite/PSScan.cgi" + payload = {"sig": motif, "db": database} + try: + LOGGER.timeit("expasySearchMotif") + result = requests.get(url, params=payload, headers=headers) + LOGGER.debug("Result taken from cache: {}".format(result.from_cache)) + LOGGER.report(msg="expasySearchMotif completed in %.2fs.", label="expasySearchMotif") + except Exception as exception: + LOGGER.error("Remote search for motif {} in Swiss-Prot database failed: {}".format(motif, exception)) + else: + return {"sp": _expasySPMotifToStructure, "pdb": _expasyPDBMotifToStructure,}.get( + database, err + )(result.text) + + +def spOnlineMotifSearch(motif: str, *args) -> list: + """Search motif in remote Swiss-Prot database. + + Args: + motif (str): Motif in PROSITE format. + + Returns: + list: List of result dicts + """ + LOGGER.debug("Search for motif in Swiss-Prot online.") + return expasySearchMotif(motif, "sp") + + +def rsOnlineMotifSearch(motif: str, *args) -> list: + """ "Search motif in remote RefSeq database. + + Args: + motif (str): Motif in PROSITE format. + + Raises: + NotImplementedError: not implemented yet + + Returns: + list: List of result dicts + """ + LOGGER.debug("Search for motif in RefSeq online.") + # return remoteSearchMotif(motif, "refseq") + raise NotImplementedError + + +def pdbOnlineMotifSearch(motif: str, *args) -> list: + """Search motif in remote PDB database. + + + Args: + motif (str): Motif in PROSITE format. + + Returns: + list: list: List of result dicts + """ + LOGGER.debug("Search for motif in Protein Data Bank online.") + return expasySearchMotif(motif, "pdb") + + +def _defaultParse(pattern: str, protein: list) -> list: + """Parse protein sequence. + + Args: + pattern (str): motif regex + protein (list): FASTA protein + + Returns: + Optional[list]: Parsed protein. + """ + description = protein[0].replace(">", "") + sequence = "".join(protein[1:]).replace("\n", "") + motif_matcher = re.compile(pattern) + results = re.finditer(motif_matcher, sequence) + elements = [] + for result in results: + elements.append( + { + "description": description, + "sequence": sequence, + "start": result.start(), + "end": result.end(), + "match": result.group(0), + } + ) + return elements + + +def _spParse(pattern: str, protein: list) -> list: + """Parse Swiss-Prot protein sequence. + + Args: + pattern (str): motif regex + protein (list): FASTA protein + + Returns: + Optional[list]: Parsed protein. + """ + line = protein[0].replace(">", "").replace("\n", "") + sequence = "".join(protein[1:]).replace("\n", "") + tmp = line.split(" ") + _, accession, id = tmp[0].split("|") + description = " ".join(tmp[1:]) if id else line + motif_matcher = re.compile(pattern) + results = re.finditer(motif_matcher, sequence) + elements = [] + for result in results: + elements.append( + { + "accession": accession, + "id": id, + "description": description, + "sequence": sequence, + "start": result.start(), + "end": result.end(), + "match": result.group(0), + } + ) + return elements + + +def _rsParse(pattern: str, protein: list) -> list: + """Parse RefSeq protein sequence. + + Args: + pattern (str): motif regex + protein (list): FASTA protein + + Returns: + Optional[list]: Parsed protein. + """ + raise NotImplementedError + + +def _pdbParse(pattern: str, protein: list) -> list: + """Parse Protein Data Bank protein sequence. + + Args: + pattern (str): motif regex + protein (list): FASTA protein + + Returns: + Optional[list]: Parsed protein. + """ + elements = [] + line = protein[0].replace(">", "") + sequence = "".join(protein[1:]) + pdb_id, mol, length, *description = re.split(r"\s+", line) + if "protein" not in mol: + return elements + length = length.split(":")[1] + motif_matcher = re.compile(pattern) + results = re.finditer(motif_matcher, sequence) + for result in results: + elements.append( + { + "PDB id": pdb_id, + "Description": " ".join(description), + "Protein length": length, + "sequence": sequence, + "start": result.start(), + "end": result.end(), + "match": result.group(0), + } + ) + return elements + + +def localMotifSearch(filename: str, motif: str, parse: str) -> list: + """Search for motif in local database file. + + Args: + filename (str): path to filename + motif (str): PROSITE motif + + Returns: + list: Parsed proteins + """ + data = [] + LOGGER.timeit("localMotifSearch") + pattern = prositeToRegEx(motif) + protein = [] + try: + file = openFile(filename) + line = file.readline() + if line.startswith(">"): + protein.append(line) + else: + LOGGER.error("Only FASTA files supported.") + for line in file: + if line.startswith(">"): + data.extend( + {"sp": _spParse, "rs": _rsParse, "pdb": _pdbParse,}.get( + parse, _defaultParse + )(pattern, protein) + ) + protein.clear() + protein.append(line) + else: + protein.append(line) + except Exception as e: + LOGGER.error(str(e)) + finally: + file.close() + LOGGER.report(msg="localMotifSearch completed in %.2fs.", label="localMotifSearch") + return data + + +def spLocalMotifSearch(motif: str, data_dir: str) -> list: + """Search for motif in local Swiss-Prot database file. + + Args: + motif (str): PROSITE motif + data_dir (str): local Swiss-Prot database directory + Returns: + list: Parsed proteins + """ + LOGGER.debug("Running local Swiss-Prot motif search...") + filename = os.path.join(data_dir, "SwissProt/uniprot_sprot.fasta.gz") + return localMotifSearch(filename, motif, "sp") + + +def rsLocalMotifSearch(motif: str, data_dir: str) -> list: + """Search for motif in local RefSeq database file. + + Args: + motif (str): PROSITE motif + data_dir (str): local RefSeq database directory + upload_dir (str): directory for local custom databases + + Raises: + NotImplementedError: not implemented yet + + Returns: + list: list: Parsed proteins + """ + LOGGER.debug("Running local RefSeq motif search...") + raise NotImplementedError + + +def pdbLocalMotifSearch(motif: str, data_dir: str) -> list: + """Search for motif in local PDB database file. + + Args: + motif (str): PROSITE motif + data_dir (str): local RefSeq database directory + upload_dir (str): directory for local custom databases + + Raises: + NotImplementedError: not implemented yet + + Returns: + list: Parsed proteins + """ + LOGGER.debug("Running local ProteinDataBank motif search...") + filename = os.path.join(data_dir, "ProteinDataBank/pdb_seqres.txt.gz") + return localMotifSearch(filename, motif, "pdb") + + +def motifSearch(database: str, motif: str, data_dir: str) -> Optional[list]: + """Search for motif in protein database. + + Args: + database (str): selected database + motif (str): motif to be searched + data_dir (str): local builin database directory + upload_dir (str): local custom database directory + + Returns: + Optional[list]: Parsed proteins + """ + + def custom(motif: str, data_dir: str) -> Optional[list]: + LOGGER.debug("Running local motif search...") + filename = os.path.join(data_dir, "local", database) + return localMotifSearch(filename, motif, "default") + + LOGGER.debug("Motif search:\t{}\t{}".format(database, motif)) + return { + "sp-online": spOnlineMotifSearch, + "rs-online": rsOnlineMotifSearch, + "pdb-online": pdbOnlineMotifSearch, + "sp-local": spLocalMotifSearch, + "rs-local": rsLocalMotifSearch, + "pdb-local": pdbLocalMotifSearch, + }.get(database, custom)(motif, data_dir) diff --git a/prody/tests/sequence/test_motif.py b/prody/tests/sequence/test_motif.py new file mode 100644 index 000000000..0002205e4 --- /dev/null +++ b/prody/tests/sequence/test_motif.py @@ -0,0 +1,102 @@ +from prody import LOGGER, _expasyMotifToStructure, expasySearchMotif +from prody.tests import TestCase + +LOGGER.verbosity = None + + +EXPASY_SCAN_PROSITE_TOOL_RESPONSE = """
include splice variants (Swiss-Prot) +Output format: Text + +Hits for USERPAT1'P-x(2)-G-E-S-G(2)-[AS]' motif on UniProtKB/Swiss-Prot sequences: +UniProtKB/Swiss-Prot (Release 2022_01 of 23-Feb-22) contains 566,996 entries. + +found: 11 hits in 11 sequences + +Graphical view (graphical view with feature detection) +shaded alignment of hits + + : + >USERPAT1 (user pattern) : + Pattern: P-x(2)-G-E-S-G(2)-[AS] + Approximate number of expected random matches in ~ 100'000 sequences (50'000'000 residues): 0.44 + + +>sp|Q3SYW2|CO2_BOVIN (750 aa) +Complement C2 (EC 3.4.21.43) (C3/C5 convertase) [Cleaved into: Complement C2b fragment; Complement C2a fragment]. [Bos taurus (Bovine)] +MDPLMAVLCLLPLYPGLATAALSCPKNVNISGGSFTLSNGWNPGSILTYSCPLGHYPYPVVTRLCKSNGQWQIPRSTRST +KAICKPVRCPAPVSFENGVYIPRLGSHPVGGNLSFECEDGFTLRGSAVRQCRPNGMWDGETAVCDNGASHCPNPGISVGA +VRTGSRFGLGDKVRYRCSSNLVLTGSAERECQDDGVWSGTEAICRQPYSYDFPEDVAPALGTSFSHLLATTNPIQQKKKQ +NLGRKIQIQRSGHLNLYLLLDASQSVSKDDFEIFKDSASRMVDRIFSFEIKVSVAIITFASKPKIIMSVLEDRSRDVTEV +ENSLRNINYKDHENGTGTNIYEALHAVYIMMNNQMNRPHMNPGAWQEIRHAIILLTDGKSNMGGSPKVAVDNIKEVLNIN +QKRKDYLDIYAIGVGSLHVDWKELNNLGSKKDGERHAFILKDVQALSQVFEHMLDVSQLTDPICGVGNMSANASAQERTP +WHVTIKPKSQETCRGALISDQWVLTAAHCFRNAEDRTLWRVSVGDPNFQGSKEFQIEEAVISPGFNVFSKKSQGIPEFYG +DDIALLKLTQKVKMSTHARPICLPCTVGANLALRKLPGSTCRDHEKELLNQVSIPAHFVALNGDKLNINLKTGSEWTNCV +KVVLKDKTTFPNLTDVREVVTDQFLCSGTQGDDSPCKGESGGAVFLERRLRFFQVGLVSWGLYNPCGGSSKNSRKPAPHG +KVPRDFHINLFRLQPWLRQHLEGILNFVPL + 675 - 683: PckGESGGA + + +>sp|Q863A0|CO2_GORGO (752 aa) +Complement C2 (EC 3.4.21.43) (C3/C5 convertase) [Cleaved into: Complement C2b fragment; Complement C2a fragment]. [Gorilla gorilla gorilla (Western lowland gorilla)] +MGPLMVLFCLLFVYTGLADSAPSCPQNVNISGGTFTLSHGWAPGSLLTYSCPQGLYPSPASRLCKSSGQWQTPGATRSLS +KAVCKPVRCPAPVSFENGIYTPRLGSYPVGGNVSFECEDGFILRGSPVRQCRPNGMWDGETAVCDNGAGHCPNPGISLGA +VRTGFRFGHGDKVRYRCSSNLVLTGSSERECQGNGVWSGTEPICRQPYSYDFPEDVAPALGTSFSHMLGATNPTQKTKES +LGRKIQIQRSGHLNLYLLLDCSQSVSENDFLIFKESASLMVDRIFSFEINVSVAIITFASKPKVLMSVLNDNSRDMTEVI +SSLENANYKDHENGTGTNTYAALNSVYLMMNNQMRILGMETMAWQEIRHAIILLTDGKSNMGGSPKTAVDRIREILNINQ +KRNDYLDIYAIGVGKLDVDWRELNELGSKKDGERHAFILQDTKALHQVFEHMLDVSKLTDTICGVGNMSANASDQERTPW +HVTIKPKSQETCRGALISDQWVLTAAHCFRDGNDHSLWRVNVGDPKSQWGKEFLIEKAVISPGFDVFAKKNQGILEFYGD +DIALLKLAQKVKMSTHARPICLPCTMEANLALRRPQGSTCRDHENELLNKQSVPAHFVALNGSKLNINLKMGVEWTSCAE +VVSQEKTMFPNLTDVREVVTDQFLCSGTQEDESPCKGESGGAVFLERRFRFFQVGLVSWGLYNPCLGSADKNSRKRAPRS +KVPPPRDFHINLFRMQPWLRQHLGDVLNFLPL + 674 - 682: PckGESGGA""" + +EXPASY_SCAN_PROSITE_TOOL_RESULT = [ + { + "database": "sp", + "accession": "Q3SYW2", + "id": "CO2_BOVIN", + "description": "Complement C2 (EC 3.4.21.43) (C3/C5 convertase) [Cleaved into: Complement C2b fragment; Complement C2a fragment]. [Bos taurus (Bovine)]", + "sequence": "MDPLMAVLCLLPLYPGLATAALSCPKNVNISGGSFTLSNGWNPGSILTYSCPLGHYPYPVVTRLCKSNGQWQIPRSTRST" + "KAICKPVRCPAPVSFENGVYIPRLGSHPVGGNLSFECEDGFTLRGSAVRQCRPNGMWDGETAVCDNGASHCPNPGISVGA" + "VRTGSRFGLGDKVRYRCSSNLVLTGSAERECQDDGVWSGTEAICRQPYSYDFPEDVAPALGTSFSHLLATTNPIQQKKKQ" + "NLGRKIQIQRSGHLNLYLLLDASQSVSKDDFEIFKDSASRMVDRIFSFEIKVSVAIITFASKPKIIMSVLEDRSRDVTEV" + "ENSLRNINYKDHENGTGTNIYEALHAVYIMMNNQMNRPHMNPGAWQEIRHAIILLTDGKSNMGGSPKVAVDNIKEVLNIN" + "QKRKDYLDIYAIGVGSLHVDWKELNNLGSKKDGERHAFILKDVQALSQVFEHMLDVSQLTDPICGVGNMSANASAQERTP" + "WHVTIKPKSQETCRGALISDQWVLTAAHCFRNAEDRTLWRVSVGDPNFQGSKEFQIEEAVISPGFNVFSKKSQGIPEFYG" + "DDIALLKLTQKVKMSTHARPICLPCTVGANLALRKLPGSTCRDHEKELLNQVSIPAHFVALNGDKLNINLKTGSEWTNCV" + "KVVLKDKTTFPNLTDVREVVTDQFLCSGTQGDDSPCKGESGGAVFLERRLRFFQVGLVSWGLYNPCGGSSKNSRKPAPHG" + "KVPRDFHINLFRLQPWLRQHLEGILNFVPL", + "start": "675", + "end": "683", + "match": "PckGESGGA", + }, + { + "database": "sp", + "accession": "Q863A0", + "id": "CO2_GORGO", + "description": "Complement C2 (EC 3.4.21.43) (C3/C5 convertase) [Cleaved into: Complement C2b fragment; Complement C2a fragment]. [Gorilla gorilla gorilla (Western lowland gorilla)]", + "sequence": "MGPLMVLFCLLFVYTGLADSAPSCPQNVNISGGTFTLSHGWAPGSLLTYSCPQGLYPSPASRLCKSSGQWQTPGATRSLS" + "KAVCKPVRCPAPVSFENGIYTPRLGSYPVGGNVSFECEDGFILRGSPVRQCRPNGMWDGETAVCDNGAGHCPNPGISLGA" + "VRTGFRFGHGDKVRYRCSSNLVLTGSSERECQGNGVWSGTEPICRQPYSYDFPEDVAPALGTSFSHMLGATNPTQKTKES" + "LGRKIQIQRSGHLNLYLLLDCSQSVSENDFLIFKESASLMVDRIFSFEINVSVAIITFASKPKVLMSVLNDNSRDMTEVI" + "SSLENANYKDHENGTGTNTYAALNSVYLMMNNQMRILGMETMAWQEIRHAIILLTDGKSNMGGSPKTAVDRIREILNINQ" + "KRNDYLDIYAIGVGKLDVDWRELNELGSKKDGERHAFILQDTKALHQVFEHMLDVSKLTDTICGVGNMSANASDQERTPW" + "HVTIKPKSQETCRGALISDQWVLTAAHCFRDGNDHSLWRVNVGDPKSQWGKEFLIEKAVISPGFDVFAKKNQGILEFYGD" + "DIALLKLAQKVKMSTHARPICLPCTMEANLALRRPQGSTCRDHENELLNKQSVPAHFVALNGSKLNINLKMGVEWTSCAE" + "VVSQEKTMFPNLTDVREVVTDQFLCSGTQEDESPCKGESGGAVFLERRFRFFQVGLVSWGLYNPCLGSADKNSRKRAPRS" + "KVPPPRDFHINLFRMQPWLRQHLGDVLNFLPL", + "start": "674", + "end": "682", + "match": "PckGESGGA", + }, +] + +EXPASY_ERROR = '
ERROR: Invalid characters in pattern "aaaaaa".
' + + +class TestMotif(TestCase): + def testExpasyTextToStructure(self): + self.assertEqual(_expasyMotifToStructure(EXPASY_SCAN_PROSITE_TOOL_RESPONSE), EXPASY_SCAN_PROSITE_TOOL_RESULT) + + def testEmptyExpasyStructure(self): + self.assertEqual(_expasyMotifToStructure(EXPASY_ERROR), []) diff --git a/prody/utilities/__init__.py b/prody/utilities/__init__.py index 9bb73a255..92f8f5cae 100644 --- a/prody/utilities/__init__.py +++ b/prody/utilities/__init__.py @@ -73,6 +73,7 @@ from .seqtools import * from .TreeConstruction import * from .eigtools import * +from .motifhelpers import * from . import catchall from .catchall import * diff --git a/prody/utilities/motifhelpers.py b/prody/utilities/motifhelpers.py new file mode 100644 index 000000000..7fafff48d --- /dev/null +++ b/prody/utilities/motifhelpers.py @@ -0,0 +1,166 @@ +# -*- coding: utf-8 -*- + +"""This module defines functions and constants related to MOTIF.""" + +from typing import Any, Optional +import re +import os +import requests +import requests_cache +from datetime import timedelta +from .logger import LOGGER + +__all__ = [ + "prositeToRegEx", + "validateMotif", + "downloadFile", + "openFile", + "getLocalDBs", + "getGenericDBs", +] + +CHUNK_SIZE = 16 * 1024 + +requests_cache.install_cache(expire_after=timedelta(weeks=1)) + + +def prositeToRegEx(motif: str) -> str: + """Change PROSITE Motif to python regular expression. + + Args: + pattern (str): PROSITE Motif + + Returns: + str: regular expression + """ + pattern = ( + motif.replace("{", "[^") + .replace("}", "]") + .replace("(", "{") + .replace(")", "}") + .replace("-", "") + .replace("x", ".") + .replace(">", "$") + .replace("<", "*") + ) + return pattern + + +def validateMotif(prosite: str) -> bool: + """_summary_ + + Args: + prosite (str): _description_ + + Returns: + bool: _description_ + """ + motif = prosite.replace("-", "") + motif_pattern = r"^([ARDNCEQGHILKMFPSTWYVx]*(\[[ARDNCEQGHILKMFPSTWYV]+\])*({{[{{ARDNCEQGHILKMFPSTWYV}}]+}})*([ARDNCEQGHILKMFPSTWYVx]\(\d(,\d)?\))*)+$" # noqa + motif_matcher = re.compile(motif_pattern) + matcher = re.match(motif_matcher, motif) + return True if matcher else False + + +def downloadFile(url: str, dir: str, file: str, **kwargs: Any) -> str: + """Downloads a file via http protocol. + + Args: + url (str): file url for download + dir (str): directory to download to + file (str): filename to download + kwargs (Any): keyword arguments for GET request + """ + kwargs.update({"stream": True}) + LOGGER.info("Downloading file {}.".format(file)) + os.makedirs(dir, exist_ok=True) + try: + with requests_cache.disabled(): + with requests.get(url + file, **kwargs) as request: + with open(os.path.join(dir, file), "wb") as output: + for chunk in request.iter_content(chunk_size=CHUNK_SIZE): + output.write(chunk) + except requests.exceptions.RequestException as exception: + LOGGER.error(str(exception)) + return file + + +def _getCompressionType(filename: str) -> Optional[str]: + """Attempts to guess the compression (if any) on a file using the first few bytes. + + Args: + filename (str): file to analize compression type + + Returns: + Optional[str]: compression type + """ + magic_dict = { + "gz": (b"\x1f", b"\x8b", b"\x08"), + "bz2": (b"\x42", b"\x5a", b"\x68"), + "zip": (b"\x50", b"\x4b", b"\x03", b"\x04"), + } + max_len = max(len(x) for x in magic_dict) + with open(filename, "rb") as unknown_file: + file_start = unknown_file.read(max_len) + compression_type = None + for file_type, magic_bytes in magic_dict.items(): + if file_start.startswith(magic_bytes): + compression_type = file_type + return compression_type + + +def openFile(filename: str) -> Any: + """Open file and return handler. + + Args: + filename (str): file to open + + Raises: + NotImplementedError: compression not supported + + Returns: + Any: file handler + """ + compression = _getCompressionType(filename) + if compression == "gz": + import gzip + + return gzip.open(filename, "rt") + elif not compression: + return open(filename, "r") + else: + raise NotImplementedError + + +def getLocalDBs(upload_dir: str) -> list: + """Get all local custom databases. + + Args: + upload_dir (str): Directory to search + + Returns: + list: databases + """ + return os.listdir(upload_dir) + + +def getGenericDBs(data_dir: str) -> dict: + """Get all buildin databases. + + Args: + data_dir (str): Directory to search + + Returns: + dict: databases + """ + result = {} + sp = "{}/SwissProt/release.txt".format(data_dir) + if os.path.isfile(sp): + result.update({"sp-local": "Swiss-Prot (local, buildin)"}) + rs = "{}/RefSeq/release.txt".format(data_dir) + if os.path.isfile(rs): + result.update({"rs-local": "RefSeq (local, buildin)"}) + pdb = "{}/ProteinDataBank/pdb_seqres.txt.gz".format(data_dir) + if os.path.isfile(pdb): + result.update({"pdb-local": "Protein Data Bank (local, buildin)"}) + return result diff --git a/pyproject.toml b/pyproject.toml index b0d88ed49..ffc2b96d2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,6 +48,8 @@ dependencies = [ "biopython", "pyparsing<=3.1.1", "scipy", + "requests", + "requests_cache", ] [project.urls] diff --git a/scripts/motif_search/app.py b/scripts/motif_search/app.py new file mode 100755 index 000000000..0d24a93b5 --- /dev/null +++ b/scripts/motif_search/app.py @@ -0,0 +1,116 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import os +from flask import Flask, render_template, request, send_file, abort, redirect, url_for, flash +from markupsafe import escape +from werkzeug.utils import secure_filename +from prody.sequence.motif import motifSearch, saveMotifResults, getUniprot +from prody.database.swissprot import SwissProt +from prody.database.refseq import RefSeq +from prody.database.pdb import ProteinDataBank +from prody.utilities.motifhelpers import validateMotif, getLocalDBs, getGenericDBs + + +app = Flask(__name__) +app.config.from_pyfile("./config.py") +app.static_folder = app.config.get("STATIC_DIR") + + +@app.route("/", methods=["GET", "POST"]) +def index(): + data_dir = app.config.get("DATA_DIR", "./data") + upload_dir = os.path.join(data_dir, "local") + os.makedirs(upload_dir, exist_ok=True) + if request.method == "POST": + if "file" not in request.files: + flash("No file was uploaded.") + return redirect(request.url) + file = request.files["file"] + if file.filename == "": + flash("Please select a file before upload.") + return redirect(request.url) + filename = secure_filename(file.filename) + file.save(os.path.join(upload_dir, filename)) + flash("File {} uploaded.".format(file)) + return redirect(request.url) + generic_dbs = getGenericDBs(data_dir) + local_dbs = getLocalDBs(upload_dir) + return render_template("index.html", generic_dbs=generic_dbs, local_dbs=local_dbs) + + +@app.route("/searchMotif") +def searchMotif(): + static_dir = app.config.get("STATIC_DIR", "./static") + data_dir = app.config.get("DATA_DIR", "./data") + upload_dir = os.path.join(data_dir, "local") + os.makedirs(upload_dir, exist_ok=True) + motif = str(escape(request.args.get("motif", default=None, type=str))) + if not validateMotif(motif): + flash('"{}" is not PROSITE motif.'.format(motif)) + return redirect(url_for("index")) + database = str(escape(request.args.get("database", default=None, type=str))) + try: + data = motifSearch(database, motif, data_dir) + except: + abort(404) + if data: + filename = saveMotifResults(database, motif, data, static_dir) + flash("Results saved to {}".format(filename)) + else: + flash("No {} pattern found in {}".format(motif, database)) + return redirect(url_for("index")) + local_dbs = getLocalDBs(upload_dir) + generic_dbs = getGenericDBs(data_dir) + return render_template( + "search.html", + generic_dbs=generic_dbs, + local_dbs=local_dbs, + data=data, + filename=filename, + motif=motif, + ) + + +@app.route("/result/| t |
Return +
to the Home Page +Found {{ data|length }} matches for {{ motif }} motif.
+ + + +{% endblock %} +{% block content %} +