From 35acf20650a96f8065b77aaacb8e5023eca261b4 Mon Sep 17 00:00:00 2001
From: HarshitGupta11 <50410275+HarshitGupta11@users.noreply.github.com>
Date: Thu, 20 Jun 2019 19:14:12 +0530
Subject: [PATCH] Add files via upload

---
 create_synteny_matrix_v2.py | 119 ++++++++++++++++++++++++++++++++++++
 test_model.py               |   2 +-
 2 files changed, 120 insertions(+), 1 deletion(-)
 create mode 100644 create_synteny_matrix_v2.py

diff --git a/create_synteny_matrix_v2.py b/create_synteny_matrix_v2.py
new file mode 100644
index 0000000..2efe511
--- /dev/null
+++ b/create_synteny_matrix_v2.py
@@ -0,0 +1,119 @@
+import numpy as np
+import requests
+import edlib as ed
+import pandas as pd
+import time
+import sys
+import progressbar
+from skbio.alignment import local_pairwise_align_ssw
+from skbio import DNA,TabularMSA,RNA
+
+def update(gene_seq,gene):
+    while(1):
+        try:
+            server = "https://rest.ensembl.org"
+            ext = "/sequence/id/"+str(gene)+"?type=cds;multiple_sequences=1"
+
+            r = requests.get(server+ext, headers={ "Content-Type" : "application/json"})
+
+            if not r.ok:
+                r.raise_for_status()
+                sys.exit()
+            r=r.json()
+            if len(r)==1:
+                r=dict(r[0])
+                gene_seq[gene]=str(r["seq"])
+                return
+            else:
+                maxi=0
+                maxlen=0
+                for i in range(len(r)):
+                    m=r[i]
+                    m=dict(m)
+                    if len(m["seq"])>maxlen:
+                        maxi=i
+                r=dict(r[maxi])
+                gene_seq[gene]=str(r["seq"])
+                return
+        except Exception as e:
+            print("\nError:",e)
+            continue
+
+def create_synteny_matrix_mul(gene_seq,g1,g2,n):
+    for gene in g1:
+        if gene=="NULL_GENE":
+            continue
+        try:
+            temp=gene_seq[gene]
+        except:
+            print("Updating gene sequences for gene:",gene)
+            update(gene_seq,gene)
+    for gene in g2:
+        if gene=="NULL_GENE":
+            continue
+        try:
+            temp=gene_seq[gene]
+        except:
+            print("Updating gene sequences for gene:",gene)
+            update(gene_seq,gene)
+    #print(n)
+    sm=np.zeros((n,n,2))
+    sml=np.zeros((n,n,2))
+    for i in range(n):
+        if g1[i]=="NULL_GENE":
+            continue
+        for j in range(n):
+            if g2[j]=="NULL_GENE":
+                continue
+            norm_len=max(len(gene_seq[g1[i]]),len(gene_seq[g2[j]]))
+            try:
+                result = ed.align(gene_seq[g1[i]],gene_seq[g2[j]], mode="NW", task="distance")
+                sm[i][j][0]=result["editDistance"]/(norm_len)
+                result = ed.align(gene_seq[g1[i]],gene_seq[g2[j]][::-1], mode="NW", task="distance")
+                sm[i][j][1]=result["editDistance"]/(norm_len)
+            except:
+                return np.zeros((n,n,2)),np.zeros((n,n,2))
+    return sm,sml
+
+def synteny_matrix(gene_seq,hdf,lsy,n,enable_break):
+    sg=[]
+    sl=[]
+    t=0
+    ind=[]
+    for index,row in progressbar.progressbar(hdf.iterrows()):
+        g1=str(row["gene_stable_id"])
+        g2=str(row["homology_gene_stable_id"])
+        x=[]
+        y=[]
+        t+=1
+        try:
+            temp=lsy[g1]
+            temp=lsy[g2]
+        except:
+            continue
+        for i in range(len(lsy[g1]['b'])-1,-1,-1):
+            x.append(lsy[g1]['b'][i])
+        x.append(g1)
+        for k in lsy[g1]['f']:
+            x.append(k)
+
+        for i in range(len(lsy[g2]['b'])-1,-1,-1):
+            y.append(lsy[g2]['b'][i])
+        y.append(g2)
+        for k in lsy[g2]['f']:
+            y.append(k)
+        
+        assert(len(x)==len(y))
+        assert(len(x)==(2*n+1))
+        smgtemp,smltemp=create_synteny_matrix_mul(gene_seq,x,y,2*n+1)
+        if np.all(smgtemp==0):
+            continue
+        sg.append(smgtemp)
+        sl.append(smltemp)
+        ind.append(index)
+        if t==5 and enable_break==1:
+            break
+    #print("Time Taken:",end-start)
+    #print("Average Time:",(end-start)/len(sg))
+    print(t)
+    return np.array(sg),np.array(sl),np.array(ind)
diff --git a/test_model.py b/test_model.py
index ddcd560..62ab0e6 100644
--- a/test_model.py
+++ b/test_model.py
@@ -6,7 +6,7 @@ import sys
 import pickle
 from get_data import get_data_genome
 from process_data import get_nearest_neighbors
-from create_synteny_matrix import create_synteny_matrix_mul
+from create_synteny_matrix_v2 import create_synteny_matrix_mul
 from ete3 import Tree
 import traceback