initial commit

Author: pvanheus

.vscode/settings.json

@@ -0,0 +1,3 @@
+{
+    "python.pythonPath": "/home/pvh/anaconda3/envs/parse_accessions/bin/python"
+}

README.md

@@ -0,0 +1,12 @@
+This project contains two sets of code:
+
+Firstly, parser for the accession format described here: https://www.ncbi.nlm.nih.gov/Sequin/acc.html called 
+`scrape_accession_rules.py`. This can be run as 
+
+```
+scrape_accession_rules.py ncbi_accession_rules.json
+```
+
+to scrape the rules at that website and save them in JSON format. The file generated by this parser is
+used in the second script, `parse_ncbi_accession.py` which can be used as a Python module but also
+as a command line tool e.g. `parse_ncbi_accession.py ncbi_accession_rules.json AE014297`.

ncbi_accession_rules.json

@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+
+import json
+import re
+
+class RuleMatcher(object):
+    def matches(this, accession):
+        pass
+
+class RangeMatcher(RuleMatcher):
+    
+    def __init__(self, prefix, suffix_length):
+        self.re = re.compile(r'{prefix}[A-Z]{suffix_length}')
+    
+    def matches(this, accession):
+        return self.re.match(accession) is not None
+
+def make_range_matcher(spec):
+    if '-' in spec:
+        (start, end) = spec.split('-')
+    elif '_' in spec:
+        (start, end) = spec.split('_')
+    else:
+        raise ValueError('require specification with - or _ in it, got: ' + spec)
+    assert end.endswith('Z'), f'Range specifier has unknown format: {spec}'
+    for i in range(len(start)):
+        if start[i] != end[i]:
+            break
+    prefix_length = i
+    assert prefix_length < len(start), f'Prefix not found in range specifier {spec}'
+    suffix_length = len(start) - prefix_length
+    return (prefix_length + suffix_length, RangeMatcher(start[:prefix_length], suffix_length))
+
+
+def build_accession_parser(rules_file):
+    rules_data = json.load(rules_file)
+    rules_by_prefix_len = {}
+    for prefix_list, database, type_description in rules_data:
+        for prefix in prefix_list:
+            prefix_length = len(prefix)
+            if '-' in prefix or '_' in prefix:
+                (prefix_length, matcher) = make_range_matcher(prefix)
+                if prefix_length not in rules_by_prefix_len:
+                    rules_by_prefix_len[prefix_length] = []
+                rules_by_prefix_len[prefix_length].append((matcher, database, type_description))
+            else:
+                if prefix_length not in rules_by_prefix_len:
+                    rules_by_prefix_len[prefix_length] = []
+                rules_by_prefix_len[prefix_length].append((prefix, database, type_description))
+    return rules_by_prefix_len
+
+letter_re = re.compile('[A-Z]+')
+def match_accession(accession, rules_by_prefix_len):
+    if accession[2] == '_':
+        return ('NCBI', 'RefSeq')
+    letter_match = letter_re.match(accession)
+    if letter_match is None:
+        raise ValueError('an accession number must start with at least one capital letter, this does not: ' + accession)
+    letter_prefix = letter_match.group(0)
+    letter_match_length = len(letter_prefix)
+    accession_type = ''
+    if letter_match_length == 0:
+        # this should never happen
+        raise ValueError('an accession number must start with at least one capital letter, this does not: ' + accession)
+    if letter_match_length < 3:
+        accession_type = 'nucleotide'
+    elif letter_match_length < 4:
+        accession_type = 'protein'
+    elif letter_match_length == 4 or letter_match_length == 6:
+        accession_type = 'WGS'
+    elif letter_match_length == 5:
+        accession_type = 'MGA'
+    else:
+        raise ValueError('an accession number must start with less than 7 capital letters, this does not: ' + accession_type)
+    
+    rules = rules_by_prefix_len[letter_match_length]
+    for rule in rules:
+        (matcher, database, type_description) = rule
+        if (isinstance(matcher, RuleMatcher) and matcher.matches(accession)) or letter_prefix == matcher:
+            return (database, type_description)
+
+if __name__ == "__main__":
+    import argparse
+    parser = argparse.ArgumentParser()
+    parser.add_argument('rules_data_file', type=argparse.FileType())
+    parser.add_argument('accession')
+    args = parser.parse_args()
+    rules_by_prefix_len = build_accession_parser(args.rules_data_file)
+    print(match_accession(args.accession, rules_by_prefix_len))

parse_ncbi_accession.py

@@ -0,0 +1,55 @@
+#!/usr/bin/env python
+
+import argparse
+import json
+from typing import List, TextIO
+
+from bs4 import BeautifulSoup
+import requests
+
+
+# def iterate_id_range(spec: str) -> List[str]:
+#     (start, end) = spec.split('-')
+#     assert end.endswith('Z'), f'Range specifier has unknown format: {spec}'
+#     for i in range(len(start)):
+#         if start[i] != end[i]:
+#             break
+#     prefix_length = i
+#     assert prefix_length < len(start), f'Prefix not found in range specifier {spec}'
+#     suffix_length = len(start) - prefix_length
+
+def parse_rules(text: str) -> List[List[str]]:
+    """parse_rules(text: str):
+       parse the rules from the NCBI webpage"""
+    soup = BeautifulSoup(text, 'html.parser')
+    data = []
+    for table in soup.find_all('table', cellpadding='3'):
+        for row in table.find_all('tr'):
+            (prefix, database, type_description, _) = row.text.split('\n')
+            if prefix.strip() == 'Prefix':
+                continue
+            if ',' in prefix:
+                prefix = [p.strip() for p in prefix.split(',')]
+            else:
+                prefix = [prefix]
+            data.append([prefix, database, type_description])
+    return data
+
+
+def fetch(url='https://www.ncbi.nlm.nih.gov/Sequin/acc.html') -> str:
+    response = requests.get(url)
+    if response.status_code == requests.codes['ok']:
+        return response.text
+    else:
+        raise requests.exceptions.RequestException(f'Failed to download {url}', response=response)
+
+
+def save_data(data: List[List[str]], output: TextIO) -> None:
+    json.dump(data, output, indent=2)
+    output.close()
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description='Save rules from NCBI website')
+    parser.add_argument('output_file', type=argparse.FileType('w'))
+    args = parser.parse_args()
+    save_data(parse_rules(fetch()), args.output_file)