Louis BECQUEY / biorseo

Go to a project
Toggle navigation Toggle navigation pinning
Louis BECQUEY
1fe1a11e 1 parent b07a46fc

minor commit

Inline Side-by-side
Showing 2 changed files with 104 additions and 100 deletions
......@@ -93,6 +93,11 @@ MOIP::MOIP(const RNA& rna, const vector<Motif>& insertionSites, float theta, boo
define_problem_constraints();
if (verbose_) cout << "A total of " << getNumConstraints(model_) << " constraints are used." << endl;
if (getNumConstraints(model_) > 2000) {
cerr << "Stopping 'cause too big for me..." << endl;
exit(-1);
}
// Define the motif objective function:
obj1 = IloExpr(env_);
for (uint i = 0; i < insertion_sites_.size(); i++) {
......@@ -170,8 +175,7 @@ SecondaryStructure MOIP::solve_objective(int o, double min, double max)
}
if (verbose_)
cout << "\t>Solution status: objective values (" << cplex_.getValue(obj1)
<< ", " << cplex_.getValue(obj2) << ')';
cout << "\t>Solution status: objective values (" << cplex_.getValue(obj1) << ", " << cplex_.getValue(obj2) << ')';
// Build a secondary Structure
SecondaryStructure best_ss = SecondaryStructure(rna_);
......@@ -301,10 +305,10 @@ void MOIP::define_problem_constraints(void)
count++;
}
}
if (count > 1) {
if (count > 0) {
model_.add(c3 <= (kxi - IloNum(2)));
if (verbose_) cout << "\t\t";
// if (verbose_) cout << x.atlas_id << '-' << j << ": ";
if (verbose_) cout << x.get_identifier() << '-' << j << ": ";
if (verbose_) cout << (c3 <= (kxi - IloNum(2))) << endl;
}
}
......
......@@ -227,102 +227,102 @@ filename = argv[1]
verbose = int(argv[3])
basename = filename[0:filename.index('.')]
# Retrieving possible 2D structrures from RNAsubopt
if (verbose):
print("Retrieving possible 2D structures from RNAsubopt...")
dbn = open(basename+"_temp.dbn", "w")
subprocess.call(["RNAsubopt", "-i", filename], stdout=dbn)
dbn.close()
dbn = open(basename+"_temp.dbn", "r")
dbn.readline()
s = dbn.readline().split(' ')[0]
structures = []
l = dbn.readline()
while l:
structures.append(l.split(' ')[0])
l = dbn.readline()
dbn.close()
subprocess.call(["rm", basename+"_temp.dbn"])
for ss in structures:
if (verbose):
print(ss)
if (verbose):
print()
# Extracting probable loops from these structures
if (verbose):
print("Extracting probable loops from these structures...")
HLs = []
ILs = []
for ss in structures:
loop_candidates = enumerate_loops(ss)
for loop_candidate in loop_candidates:
if len(loop_candidate) == 1 and loop_candidate not in HLs:
HLs.append(loop_candidate)
if len(loop_candidate) == 2 and loop_candidate not in ILs:
ILs.append(loop_candidate)
if (verbose):
print("TOTAL:")
print(len(HLs), "probable hairpin loops found")
print(len(ILs), "probable internal loops")
print()
# Retrieve subsequences corresponding to the possible loops
if (verbose):
print("Retrieving subsequences corresponding to the possible loops...")
loops = []
for i, l in enumerate(HLs):
loops.append(Loop(">HL%d" % (i+1), s[l[0][0]-1:l[0][1]], "h", l))
if (verbose):
print()
print(loops[-1].get_header(), "\t\t", l)
print(loops[-1].subsequence())
for i, l in enumerate(ILs):
loops.append(
Loop(">IL%d" % (i+1), s[l[0][0]-1:l[0][1]]+'*'+s[l[1][0]-1:l[1][1]], "i", l))
if (verbose):
print()
print(loops[-1].get_header(), "\t\t", l)
print(loops[-1].subsequence())
if (verbose):
print()
# Scanning loop subsequences against motif database
if (verbose):
print("Scanning loop subsequences against motif database (using %d threads)..." %
(cpu_count()-1))
pool = Pool(processes=cpu_count()-1)
insertion_sites = [x for y in pool.map(launchJar3d, loops) for x in y]
insertion_sites.sort(reverse=True)
if (verbose):
print(len(insertion_sites), "insertions found:")
# Writing results to file
c = 0
resultsfile = open(basename+".sites.csv", "w")
resultsfile.write("Motif,Rotation,Score,Start1,End1,Start2,End2\n")
for site in insertion_sites:
if site.score > 10:
c += 1
string = "FOUND with score %d:\t\t possible insertion of motif " % site.score + site.atlas_id
if site.rotation:
string += " (reversed)"
string += (" on " + site.loop.get_header() + " at positions")
if (verbose):
print(string, site.loop.subsequence(),
'*'.join([str(x) for x in site.position]))
resultsfile.write(site.atlas_id+',' +
str(bool(site.rotation))+",%d" % site.score+',')
positions = [','.join([str(y) for y in x]) for x in site.position]
if len(positions) == 1:
positions.append("-,-")
resultsfile.write(','.join(positions)+'\n')
if c < len(insertion_sites):
if (verbose):
print("... and %d more with score(s) below 10." %
(len(insertion_sites)-c))
resultsfile.close()
# # Retrieving possible 2D structrures from RNAsubopt
# if (verbose):
# print("Retrieving possible 2D structures from RNAsubopt...")
# dbn = open(basename+"_temp.dbn", "w")
# subprocess.call(["RNAsubopt", "-i", filename], stdout=dbn)
# dbn.close()
# dbn = open(basename+"_temp.dbn", "r")
# dbn.readline()
# s = dbn.readline().split(' ')[0]
# structures = []
# l = dbn.readline()
# while l:
# structures.append(l.split(' ')[0])
# l = dbn.readline()
# dbn.close()
# subprocess.call(["rm", basename+"_temp.dbn"])
# for ss in structures:
# if (verbose):
# print(ss)
# if (verbose):
# print()
# # Extracting probable loops from these structures
# if (verbose):
# print("Extracting probable loops from these structures...")
# HLs = []
# ILs = []
# for ss in structures:
# loop_candidates = enumerate_loops(ss)
# for loop_candidate in loop_candidates:
# if len(loop_candidate) == 1 and loop_candidate not in HLs:
# HLs.append(loop_candidate)
# if len(loop_candidate) == 2 and loop_candidate not in ILs:
# ILs.append(loop_candidate)
# if (verbose):
# print("TOTAL:")
# print(len(HLs), "probable hairpin loops found")
# print(len(ILs), "probable internal loops")
# print()
# # Retrieve subsequences corresponding to the possible loops
# if (verbose):
# print("Retrieving subsequences corresponding to the possible loops...")
# loops = []
# for i, l in enumerate(HLs):
# loops.append(Loop(">HL%d" % (i+1), s[l[0][0]-1:l[0][1]], "h", l))
# if (verbose):
# print()
# print(loops[-1].get_header(), "\t\t", l)
# print(loops[-1].subsequence())
# for i, l in enumerate(ILs):
# loops.append(
# Loop(">IL%d" % (i+1), s[l[0][0]-1:l[0][1]]+'*'+s[l[1][0]-1:l[1][1]], "i", l))
# if (verbose):
# print()
# print(loops[-1].get_header(), "\t\t", l)
# print(loops[-1].subsequence())
# if (verbose):
# print()
# # Scanning loop subsequences against motif database
# if (verbose):
# print("Scanning loop subsequences against motif database (using %d threads)..." %
# (cpu_count()-1))
# pool = Pool(processes=cpu_count()-1)
# insertion_sites = [x for y in pool.map(launchJar3d, loops) for x in y]
# insertion_sites.sort(reverse=True)
# if (verbose):
# print(len(insertion_sites), "insertions found:")
# # Writing results to file
# c = 0
# resultsfile = open(basename+".sites.csv", "w")
# resultsfile.write("Motif,Rotation,Score,Start1,End1,Start2,End2\n")
# for site in insertion_sites:
# if site.score > 10:
# c += 1
# string = "FOUND with score %d:\t\t possible insertion of motif " % site.score + site.atlas_id
# if site.rotation:
# string += " (reversed)"
# string += (" on " + site.loop.get_header() + " at positions")
# if (verbose):
# print(string, site.loop.subsequence(),
# '*'.join([str(x) for x in site.position]))
# resultsfile.write(site.atlas_id+',' +
# str(bool(site.rotation))+",%d" % site.score+',')
# positions = [','.join([str(y) for y in x]) for x in site.position]
# if len(positions) == 1:
# positions.append("-,-")
# resultsfile.write(','.join(positions)+'\n')
# if c < len(insertion_sites):
# if (verbose):
# print("... and %d more with score(s) below 10." %
# (len(insertion_sites)-c))
# resultsfile.close()
# Lauching biominserter to get 2D predictions
subprocess.call([bminDir+"/bin/biominserter",
......