Add files via upload

Author: eukaryo

directpath-comparison-experiment.py

@@ -1,6 +1,7 @@
 
 import re
 from czno import getMinBarrierPathDijkstra, basePairList, energyOfStr
+from previous_direct_methods import MorganHiggs1998GreedyDirect, Voss2004GreedyDirect
 
 def HammingDistance(structure1, structure2):
     bp1 = set(basePairList(structure1))
@@ -33,14 +34,24 @@ def PathwayToBarrier(sequence, pathway):
     return max([energyOfStr(sequence, s) for s in pathway])
 
 
-def SingleExperiment(sequence, structure1, structure2):
+def DirectPathSingleExperiment(sequence, structure1, structure2):
 
     result = []
 
-    pathway = getMinBarrierPathDijkstra(sequence, structure1, structure2)
-    result.append(max([energyOfStr(sequence, s) for s in pathway]))
+    pathway1 = getMinBarrierPathDijkstra(sequence, structure1, structure2)
+    result.append(PathwayToBarrier(sequence, pathway1))
+    
+    pathway2 = MorganHiggs1998GreedyDirect(sequence, structure1, structure2)
+    result.append(PathwayToBarrier(sequence, pathway2))
+
+    pathway3 = Voss2004GreedyDirect(sequence, structure1, structure2)
+    result.append(PathwayToBarrier(sequence, pathway3))
 
-    #TODO: ここで、先行手法をsubprocessで叩いて結果を得てエネルギーバリアを求める
+    best_value4 = float("inf")
+    for i in range(1000):
+        pathway4 = Voss2004GreedyDirect(sequence, structure1, structure2, 10, i)
+        best_value4 = min(best_value4, PathwayToBarrier(sequence, pathway4))
+    result.append(best_value4)
 
     return result
 
@@ -52,6 +63,6 @@ def SingleExperiment(sequence, structure1, structure2):
             for j in range(i+1, len(data[1])):
                 hamdist = HammingDistance(data[1][i], data[1][j])
                 if 5 <= hamdist and hamdist <= 20:
-                    result = SingleExperiment(data[0], data[1][i], data[1][j])
+                    result = DirectPathSingleExperiment(data[0], data[1][i], data[1][j])
                     print(str(len(data[0]))+" "+str(hamdist)+" "+" ".join([str(x) for x in result]))
 

previous_direct_methods.py

@@ -0,0 +1,50 @@
+
+import random
+from czno import energyOfStr, basePairList, isExclusive
+
+def basePairSetToDotNotation(sequence, bp):
+    answer = ["." for _ in sequence]
+    for x in bp:
+        assert x[0] < x[1]
+        answer[x[0]]="("
+        answer[x[1]]=")"
+    return "".join(answer)
+
+
+
+def MorganHiggs1998GreedyDirect(sequence, structure1, structure2):
+    A = set(basePairList(structure1))
+    B = set(basePairList(structure2))
+    answer = [structure1]
+    while A != B:
+        NeedToClose = sorted(list(B - A))
+        NeedToOpen = sorted(list(A - B))
+        if len(NeedToClose) >= 1:
+            NextClose = min([(sum([isExclusive(x,y) for y in NeedToOpen]), x) for x in NeedToClose])[1]
+            NextOpens = [y for y in NeedToOpen if isExclusive(NextClose, y)]
+            for x in NextOpens:
+                A.remove(x)
+                answer.append(basePairSetToDotNotation(sequence, A))
+            A.add(NextClose)
+            answer.append(basePairSetToDotNotation(sequence, A))
+        else:
+            for x in NeedToOpen:
+                A.remove(x)
+                answer.append(basePairSetToDotNotation(sequence, A))
+    return answer
+
+def Voss2004GreedyDirect(sequence, structure1, structure2, k = 1, seed = 12345):
+    random.seed(seed)
+    A = set(basePairList(structure1))
+    B = set(basePairList(structure2))
+    answer = [structure1]
+    while A != B:
+        NeedToClose = sorted(list(B - A))
+        NeedToOpen = sorted(list(A - B))
+        PossibleClose = [x for x in NeedToClose if sum([isExclusive(x,y) for y in NeedToOpen]) == 0]
+        PossibleChange = NeedToOpen + PossibleClose
+        PossibleNextState = [basePairSetToDotNotation(sequence, A^set([x])) for x in PossibleChange]
+        PossibleNextStateAndEnergy = [(energyOfStr(sequence, x), x) for x in PossibleNextState]
+        NextState = random.choice(sorted(PossibleNextStateAndEnergy)[0:k])[1]
+        answer.append(NextState)
+    return answer

random_dataset_generator.py

@@ -41,11 +41,7 @@ def basePairList(secStr):
 
 def randomSeqGenerator(length, seed):
     random.seed(seed)
-    base = ["A","C","G","U"]
-    seq = ""
-    for i in range(length):
-        seq += base[random.randint(0,3)]
-    return seq
+    return "".join([random.choice(["A","C","G","U"]) for _ in range(length)])
 
 def isoLenGenerator(length, dist_lower, dist_upper, n):
     result = [set([]) for i in range(dist_lower, dist_upper+1)]

research.sh

@@ -5,3 +5,5 @@ python3 iso_len_time_experiment.py dataset-iso-len.txt > result-isolen-time.txt
 python3 random_dataset_generator.py isodist > dataset-iso-dist.txt
 python3 iso_dist_time_experiment.py dataset-iso-dist.txt > result-isodist-time.txt
 python3 s151Rfam-localminima-generator.py > s151-localminima-dataset.txt
+tar cvJf s151-localminima-dataset.txt.tar.xz s151-localminima-dataset.txt
+tar xvJf s151-localminima-dataset.txt.tar.xz