<commit_before>from django.core.management.base import BaseCommand, CommandError
from django.conf import settings
from django.db import connection
from django.utils.text import slugify
from django.db import IntegrityError

from build.management.commands.base_build import Command as BaseBuild
from protein.models import (Protein, ProteinConformation, ProteinState, ProteinAnomaly, ProteinAnomalyType,
    ProteinSegment)
from residue.models import ResidueGenericNumber, ResidueNumberingScheme, Residue, ResidueGenericNumberEquivalent
from common.models import WebLink, WebResource, Publication
from structure.models import (Structure, StructureType, StructureSegment, StructureStabilizingAgent,PdbData,
    Rotamer, StructureSegmentModeling, StructureCoordinates, StructureCoordinatesDescription, StructureEngineering,
    StructureEngineeringDescription, Fragment)
from construct.functions import *
#from structure.functions import BlastSearch
from Bio.PDB import PDBParser,PPBuilder
from Bio import pairwise2

from structure.assign_generic_numbers_gpcr import GenericNumbering
from ligand.models import Ligand, LigandType, LigandRole, LigandProperities
from interaction.models import *
from interaction.views import runcalculation,parsecalculation

import logging
import os
import re
import yaml
import time
from collections import OrderedDict
import json
from urllib.request import urlopen
from Bio.PDB import parse_pdb_header


## FOR VIGNIR ORDERED DICT YAML IMPORT/DUMP
_mapping_tag = yaml.resolver.BaseResolver.DEFAULT_MAPPING_TAG
def dict_constructor(loader, node):
    return OrderedDict(loader.construct_pairs(node))

def represent_ordereddict(dumper, data):
    value = []

    for item_key, item_value in data.items():
        node_key = dumper.represent_data(item_key)
        node_value = dumper.represent_data(item_value)

        value.append((node_key, node_value))

    return yaml.nodes.MappingNode(u'tag:yaml.org,2002:map', value)

yaml.add_representer(OrderedDict, represent_ordereddict)
yaml.add_constructor(_mapping_tag, dict_constructor)

class Command(BaseBuild):
    help = 'Reads source data and creates pdb structure records'

    def add_arguments(self, parser):
        parser.add_argument('-p', '--proc',
            type=int,
            action='store',
            dest='proc',
            default=1,
            help='Number of processes to run')
        parser.add_argument('-f', '--filename',
            action='append',
            dest='filename',
            help='Filename to import. Can be used multiple times')
        parser.add_argument('-u', '--purge',
            action='store_true',
            dest='purge',
            default=False,
            help='Purge existing records')

    # source file directory
    structure_data_dir = os.sep.join([settings.DATA_DIR, 'structure_data', 'structures'])
    pdb_data_dir = os.sep.join([settings.DATA_DIR, 'structure_data', 'pdbs'])

    # read source files
    filenames = os.listdir(structure_data_dir)

    ### USE below to fix seg ends
    xtal_seg_end_file = os.sep.join([settings.DATA_DIR, 'structure_data', 'annotation', 'mod_xtal_segends.yaml'])
    with open(xtal_seg_end_file, 'r') as f:
        xtal_seg_ends = yaml.load(f)

    xtal_anomalies_file = os.sep.join([settings.DATA_DIR, 'structure_data', 'annotation', 'all_anomalities.yaml'])
    with open(xtal_anomalies_file, 'r') as f2:
        xtal_anomalies = yaml.load(f2)

    s = ProteinSegment.objects.all()
    segments = {}
    for segment in s:
        segments[segment.slug] = segment

    #Must delete all these first for bulk to work
    # PdbData.objects.all().delete()


    def handle(self, *args, **options):
        # delete any existing structure data
        if options['purge']:
            try:
                self.purge_structures()
            except Exception as msg:
                print(msg)
                self.logger.error(msg)

        # where filenames specified?
        if options['filename']:
            self.filenames = options['filename']

        try:
            self.logger.info('CREATING STRUCTURES')
            # run the function twice (once for representative structures, once for non-representative)
            iterations = 2
            for i in range(1,iterations+1):
                self.prepare_input(options['proc'], self.filenames, i)

            self.logger.info('COMPLETED CREATING STRUCTURES')
        except Exception as msg:
            print(msg)
            self.logger.error(msg)

    def purge_structures(self):
        Structure.objects.all().delete()
        ResidueFragmentInteraction.objects.all().delete()
        ResidueFragmentInteractionType.objects.all().delete()
        StructureLigandInteraction.objects.all().delete()
        #Remove previous Rotamers/Residues to prepare repopulate
        Fragment.objects.all().delete()
        Rotamer.objects.all().delete()
        PdbData.objects.all().delete()

    def create_rotamers(self, structure, pdb_path,d):
        wt_lookup = {} #used to match WT seq_number to WT residue record
        pdbseq = {} #used to keep track of pdbseq residue positions vs index in seq
        ref_positions = {} #WT postions in alignment
        mapped_seq = {} # index in contruct, tuple of AA and WT [position,AA]
        debug = False
        preferred_chain = structure.preferred_chain

        if len(preferred_chain.split(','))>1: #if A,B
            preferred_chain = preferred_chain.split(',')[0]


        AA = {'ALA':'A', 'ARG':'R', 'ASN':'N', 'ASP':'D',
     'CYS':'C', 'GLN':'Q', 'GLU':'E', 'GLY':'G',
     'HIS':'H', 'ILE':'I', 'LEU':'L', 'LYS':'K',
     'MET':'M', 'PHE':'F', 'PRO':'P', 'SER':'S',
     'THR':'T', 'TRP':'W', 'TYR':'Y', 'VAL':'V'}

        atom_num_dict = {'E':9, 'S':6, 'Y':12, 'G':4, 'A':5, 'V':7, 'M':8, 'L':8, 'I':8, 'T':7, 'F':11, 'H':10, 'K':9, 
                         'D':8, 'C':6, 'R':11, 'P':7, 'Q':9, 'N':8, 'W':14}


        entry_name = d['construct_crystal']['uniprot']

        # print(d['xml_segments'])
        # print(d['deletions'])
        deletions = []
        ## Remove WT ranges that arent in xtal (no need to try to map them)
        #print(d)
        if 'deletions' in d:
            for del_range in d['deletions']:
                if del_range['start']==146 and structure.pdb_code.index=='4K5Y':
                    #Manual fix for faulty 4K5Y annotation
                    continue
                for i in range(del_range['start'],del_range['end']+1):
                    deletions.append(i)
            #print("Annotation missing WT residues",d['deletions'])


        removed = []
        ## Remove segments that arent receptor (tags, fusion etc)
        if 'xml_segments' in d:
            for seg in d['xml_segments']:
                #print(seg)
                if seg[1]:
                    if seg[1][0]!=entry_name:
                        if seg[0].split("_")[1]==preferred_chain:
                            #print(seg[2],seg[3]+1)
                            #for i in range(seg[2],seg[3]+1):
                            # print(seg)
                            for i in seg[6]:
                                removed.append(i)
        
        if len(deletions)>len(d['wt_seq'])*0.9:
            #if too many deltions
            removed = []
            deletions = []

        s = PDBParser(PERMISSIVE=True, QUIET=True).get_structure('ref', pdb_path)[0]
        chain = s[preferred_chain] #select only one chain (avoid n-mer receptors)
        ppb=PPBuilder()
        seq = ''
        i = 1

        check_1000 = 0
        prev_id = 0
        bigjump = False
        for pp in ppb.build_peptides(chain): #remove >1000 pos (fusion protein / gprotein)
            for i,res in enumerate(pp,1 ):
                id = res.id
                residue_id = res.get_full_id()
                # print(i,id[1],AA[res.resname])
                if id[1] in removed:
                    chain.detach_child(id)
                    # print("removed")
                    continue
                # if id[1]<600:
                #     check_1000 += 1
                #     #need check_1000 to catch structures where they lie in 1000s (4LDE, 4LDL, 4LDO, 4N4W, 4QKX)
                # if structure.pdb_code.index in ["4RWD","3SN6","4L6R"] and id[1]>1000:
                #     last_valid = 0
                #     bigjump = True
                #     removed.append(id[1])<commit_msg>Changes to build_structures script to ignore SIFTS annotation to remove "known" non-wt residues, as it's often wrong...<commit_after>from django.core.management.base import BaseCommand, CommandError
from django.conf import settings
from django.db import connection
from django.utils.text import slugify
from django.db import IntegrityError

from build.management.commands.base_build import Command as BaseBuild
from protein.models import (Protein, ProteinConformation, ProteinState, ProteinAnomaly, ProteinAnomalyType,
    ProteinSegment)
from residue.models import ResidueGenericNumber, ResidueNumberingScheme, Residue, ResidueGenericNumberEquivalent
from common.models import WebLink, WebResource, Publication
from structure.models import (Structure, StructureType, StructureSegment, StructureStabilizingAgent,PdbData,
    Rotamer, StructureSegmentModeling, StructureCoordinates, StructureCoordinatesDescription, StructureEngineering,
    StructureEngineeringDescription, Fragment)
from construct.functions import *
#from structure.functions import BlastSearch
from Bio.PDB import PDBParser,PPBuilder
from Bio import pairwise2

from structure.assign_generic_numbers_gpcr import GenericNumbering
from ligand.models import Ligand, LigandType, LigandRole, LigandProperities
from interaction.models import *
from interaction.views import runcalculation,parsecalculation

import logging
import os
import re
import yaml
import time
from collections import OrderedDict
import json
from urllib.request import urlopen
from Bio.PDB import parse_pdb_header


## FOR VIGNIR ORDERED DICT YAML IMPORT/DUMP
_mapping_tag = yaml.resolver.BaseResolver.DEFAULT_MAPPING_TAG
def dict_constructor(loader, node):
    return OrderedDict(loader.construct_pairs(node))

def represent_ordereddict(dumper, data):
    value = []

    for item_key, item_value in data.items():
        node_key = dumper.represent_data(item_key)
        node_value = dumper.represent_data(item_value)

        value.append((node_key, node_value))

    return yaml.nodes.MappingNode(u'tag:yaml.org,2002:map', value)

yaml.add_representer(OrderedDict, represent_ordereddict)
yaml.add_constructor(_mapping_tag, dict_constructor)

class Command(BaseBuild):
    help = 'Reads source data and creates pdb structure records'

    def add_arguments(self, parser):
        parser.add_argument('-p', '--proc',
            type=int,
            action='store',
            dest='proc',
            default=1,
            help='Number of processes to run')
        parser.add_argument('-f', '--filename',
            action='append',
            dest='filename',
            help='Filename to import. Can be used multiple times')
        parser.add_argument('-u', '--purge',
            action='store_true',
            dest='purge',
            default=False,
            help='Purge existing records')

    # source file directory
    structure_data_dir = os.sep.join([settings.DATA_DIR, 'structure_data', 'structures'])
    pdb_data_dir = os.sep.join([settings.DATA_DIR, 'structure_data', 'pdbs'])

    # read source files
    filenames = os.listdir(structure_data_dir)

    ### USE below to fix seg ends
    xtal_seg_end_file = os.sep.join([settings.DATA_DIR, 'structure_data', 'annotation', 'mod_xtal_segends.yaml'])
    with open(xtal_seg_end_file, 'r') as f:
        xtal_seg_ends = yaml.load(f)

    xtal_anomalies_file = os.sep.join([settings.DATA_DIR, 'structure_data', 'annotation', 'all_anomalities.yaml'])
    with open(xtal_anomalies_file, 'r') as f2:
        xtal_anomalies = yaml.load(f2)

    s = ProteinSegment.objects.all()
    segments = {}
    for segment in s:
        segments[segment.slug] = segment

    #Must delete all these first for bulk to work
    # PdbData.objects.all().delete()


    def handle(self, *args, **options):
        # delete any existing structure data
        if options['purge']:
            try:
                self.purge_structures()
            except Exception as msg:
                print(msg)
                self.logger.error(msg)

        # where filenames specified?
        if options['filename']:
            self.filenames = options['filename']

        try:
            self.logger.info('CREATING STRUCTURES')
            # run the function twice (once for representative structures, once for non-representative)
            iterations = 2
            for i in range(1,iterations+1):
                self.prepare_input(options['proc'], self.filenames, i)

            self.logger.info('COMPLETED CREATING STRUCTURES')
        except Exception as msg:
            print(msg)
            self.logger.error(msg)

    def purge_structures(self):
        Structure.objects.all().delete()
        ResidueFragmentInteraction.objects.all().delete()
        ResidueFragmentInteractionType.objects.all().delete()
        StructureLigandInteraction.objects.all().delete()
        #Remove previous Rotamers/Residues to prepare repopulate
        Fragment.objects.all().delete()
        Rotamer.objects.all().delete()
        PdbData.objects.all().delete()

    def create_rotamers(self, structure, pdb_path,d):
        wt_lookup = {} #used to match WT seq_number to WT residue record
        pdbseq = {} #used to keep track of pdbseq residue positions vs index in seq
        ref_positions = {} #WT postions in alignment
        mapped_seq = {} # index in contruct, tuple of AA and WT [position,AA]
        debug = False
        preferred_chain = structure.preferred_chain

        if len(preferred_chain.split(','))>1: #if A,B
            preferred_chain = preferred_chain.split(',')[0]


        AA = {'ALA':'A', 'ARG':'R', 'ASN':'N', 'ASP':'D',
     'CYS':'C', 'GLN':'Q', 'GLU':'E', 'GLY':'G',
     'HIS':'H', 'ILE':'I', 'LEU':'L', 'LYS':'K',
     'MET':'M', 'PHE':'F', 'PRO':'P', 'SER':'S',
     'THR':'T', 'TRP':'W', 'TYR':'Y', 'VAL':'V'}

        atom_num_dict = {'E':9, 'S':6, 'Y':12, 'G':4, 'A':5, 'V':7, 'M':8, 'L':8, 'I':8, 'T':7, 'F':11, 'H':10, 'K':9, 
                         'D':8, 'C':6, 'R':11, 'P':7, 'Q':9, 'N':8, 'W':14}


        entry_name = d['construct_crystal']['uniprot']

        # print(d['xml_segments'])
        # print(d['deletions'])
        deletions = []
        ## Remove WT ranges that arent in xtal (no need to try to map them)
        #print(d)
        if 'deletions' in d:
            for del_range in d['deletions']:
                if del_range['start']==146 and structure.pdb_code.index=='4K5Y':
                    #Manual fix for faulty 4K5Y annotation
                    continue
                for i in range(del_range['start'],del_range['end']+1):
                    deletions.append(i)
            #print("Annotation missing WT residues",d['deletions'])


        removed = []
        ## Remove segments that arent receptor (tags, fusion etc)
        if 'xml_segments' in d:
            for seg in d['xml_segments']:
                #print(seg)
                if seg[1]:
                    if seg[1][0]!=entry_name:
                        if seg[0].split("_")[1]==preferred_chain:
                            #print(seg[2],seg[3]+1)
                            #for i in range(seg[2],seg[3]+1):
                            # print(seg)
                            for i in seg[6]:
                                removed.append(i)
        # Reset removed, since it causes more problems than not
        removed = []
        if len(deletions)>len(d['wt_seq'])*0.9:
            #if too many deltions
            removed = []
            deletions = []

        s = PDBParser(PERMISSIVE=True, QUIET=True).get_structure('ref', pdb_path)[0]
        chain = s[preferred_chain] #select only one chain (avoid n-mer receptors)
        ppb=PPBuilder()
        seq = ''
        i = 1

        check_1000 = 0
        prev_id = 0
        bigjump = False
        for pp in ppb.build_peptides(chain): #remove >1000 pos (fusion protein / gprotein)
            for i,res in enumerate(pp,1 ):
                id = res.id
                residue_id = res.get_full_id()
                # print(i,id[1],AA[res.resname])
                if id[1] in removed:
                    chain.detach_child(id)
                    # print("removed")
                    continue
                # if id[1]<600:
                #     check_1000 += 1
                #     #need check_1000 to catch structures where they lie in 1000s (4LDE, 4LDL, 4LDO, 4N4W, 4QKX)
                # if structure.pdb_code.index in ["4RWD","3SN6","4L6R"] and id[1]>1000:
                #     last_valid = 0
                #     bigjump = True
                #     removed.append(id[1])
{'fraction_non_alphanumeric': 0.11781676872350448, 'fraction_numerical': 0.012953367875647668, 'mean_word_length': 2.4826737748615955, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 6}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3097598', 'n_tokens_mistral': 5388, 'n_tokens_neox': 4992, 'n_words': 1229}