Difference between revisions of "Team:EPF Lausanne/Software"

"""
Copyright (C) 2015 iGEM Team EPF_Lausanne

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>
"""

import numpy as np
from scipy.integrate import *

class Solver:
    """
    Class that allows the solution of a system of non-linear ODEs. The system is specified by the function fun

        dy/dt = fun(t,y)

    where t is a number and y and dy/dt are numpy arrays or lists.

    The solution is performed with the dopri5 method, an explicit Runge-Kutta method of order (4)5.
    The method is due to Dormand & Prince, and is implemented by E. Hairer and G. Wanner.

    See
        http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.ode.html
    for more details.

    NOTE:
    Our Solver can take a function of the form

        f(t,y,pars)

    where PARS are parameters. PARS can be eventually passed to the constructor of the Solver.
    """

    def __init__(self,dt,fun,t0,T,y0,pars=[]):
        self.dt = dt # Time step
        self.fun = fun # Function representing the ODE
        self.t0 = t0 # Initial time
        self.T = T # Final time
        self.y0 = y0 # Initial condition

        self.pars = pars # Parameters of the system of ODEs

    def solve(self):
        """
        Solve the system of ODEs

            dy/dt = fun(t,y)

        on the interval [t0,T], with the initial condition y(0)=y0.

        Returns two lists, time and solution, containing time points and the solution at these time points.
        """

        # Choose integrator type
        r = ode(self.fun).set_integrator('dopri5')

        # Initialize the integrator
        r.set_initial_value(self.y0, self.t0)

        # Set parameters for the ODE function
        r.set_f_params(*self.pars)

        # Initialize solution list and time points list
        solution = np.asarray(self.y0)
        time = np.asarray(self.t0)

        while r.successful() and r.t < self.T:
            r.integrate(r.t + self.dt) # Perform one integration step, i.e. obtain the solution y at time t+dt

            time = np.append(time,r.t) # Append the new time
            solution = np.vstack((solution,r.y)) # Append the new solution

        return time, solution # Return time and solution vectors

    def solve_for_t(self,t):
        """
        Solve the system of ODEs

            dy/dt = fun(t,y)

        on the interval [t0,T], with the initial condition y(0)=y0.

        Returns two lists, solution and time, containing time points and the solution at these time points.

        The solution is computed at the points specified in t, i.e. the time step dt is ignored.
        """

        # Choose integrator type: dopri5 in this case
        r = ode(self.fun).set_integrator('dopri5')

        # Initialize the integrator
        r.set_initial_value(self.y0, self.t0)

        # Set parameters for the ODE function
        r.set_f_params(*self.pars)

        # Initialize solution list and time points list
        solution = []
        time = []

        for tt in t:
            r.integrate(tt) # Perform one integration step

            time.append(tt) # Append the new time
            solution.append(r.y) # Append the new solution

        return np.asarray(time), np.asarray(solution) # Return time and solution vectors


if __name__ == "__main__":
    """
    Our test functions:
       rapid_equilibrium (standard function)
       rapid_equilibrium_from_string() (returns a function compiled from a string)
    """

    import matplotlib.pylab as plt
    from test import * # Import test functions for the ODE integrator

    dt = 0.1

    t0 = 0
    T = 100
    y0 = [1,0,0]

    # Store the funtion compiled from a string
    rapid_equilibrium_s = rapid_equilibrium_from_string()

    mysolver = Solver(dt,rapid_equilibrium,t0,T,y0)
    mysolver_string = Solver(dt,rapid_equilibrium_s,t0,T,y0)

    t,y = mysolver.solve()
    tt,yt = mysolver.solve_for_t(np.linspace(t0,T,10))
    ts,ys = mysolver_string.solve()

    plt.plot(t,y)
    plt.plot(tt,yt,'x')
    plt.plot(ts,ys)
    plt.show()

"""
Copyright (C) 2015 iGEM Team EPF_Lausanne

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>
"""

import numpy as np
import warnings as w
from integrator import Solver
from scipy.optimize import minimize

class ODEFit:
    """
    Class that allows the fit of the solution (experimental data) of a system of ODEs

        dy/dt = fun(t,y,pars)

    where PARS are the parameters to fit.

    This class relay on the ODE Solver class in order to solve the system of ODEs.
    """

    def __init__(self,fun,t,y,y0,pars0):
        """
        The constructor need the function FUN defining the system of ODEs.

        t and y are experimental data (y can contain a list with the data of different species).

        y0 is the initial guess for the initial conditions.
        pars0 is the initial guess for the parameters.


        NOTE:
        The initial condition y0 is free, i.e. it is a parameter as well.
        """

        # Function defining the system of ODEs
        self.fun = fun

        # Initial guess for parameters and initial conditions
        self.pars0 =  pars0 + y0 # Merge all parameters to optimize in a single list

        # Data to fit (convert to Numpy arrays)
        self.t = np.asarray(t)
        self.y = np.asarray(y)

        # Number of ODEs composing the system
        self.n = self.y.shape[-1]

        # Fitted parameters (initialized to initial guess)
        self.pars = self.pars0

        # Index of separation between real parameters and (free) initial conditions
        self.index = len(self.pars)-self.n

        # Keep trak of FIT routine call (if True, self.pars contains the optimized parameters)
        self.fitted = False

    def _solve_for_pars(self,pars):
        """
        Solve the system of ODEs for a given set of parameters.

        Solution is computed at experimental data points.
        """

        # pars[:self.index] : fun parameters
        # pars[self.index:] : inital conditions (last entries)
        mysolver = Solver(0.1,self.fun,self.t[0],self.t[-1],pars[self.index:],pars[:self.index])

        # Solve the system of ODEs at experimental data points
        t,y = mysolver.solve_for_t(self.t)

        return np.asarray(t),np.asarray(y)

    def _resid(self,y_model,y_data):
        """
        Compute the square of the difference between the ODE solution and the data to fit.
        """
        return ((y_model-y_data)**2).sum()

    def _score(self,pars):
        """
        The score is defined as the sum of residuals for each ODE composing the system.

        The score function is minimized in order to optimize the parameters (least squares of residuals).
        """

        # Solve the system of ODEs for a given set of parameters
        t,y_model = self._solve_for_pars(pars)

        score = 0

        # Compute the score of the current set of parameters
        for i in range(self.n):
            score += self._resid(self.y[:,i],y_model[:,i]) # The score is the sum of residuals for each ODE

        return score

    def fit(self):
        """
        Fit data in order to obtain the best parameter estimation.

        The best parameter set is obtained by minimizing the score function, i.e. the sum of residuals.
        """

        opt = minimize(self._score,self.pars0)

        # Extract the optimized parameters
        self.pars = opt.x

        # The fit function has been called
        self.fitted = True

        if opt.success == True: # Check the success of the minimization process
            return opt.x[:self.index], opt.x[self.index:]
        else:
            raise RuntimeError(opt.message)

    def solve_fitted(self,dt=0.1):
        """
        Solve the system of ODEs with the optimal parameters with a dense output.

        dt is the discretization step, which can be changed manually.
        """

        if self.fitted == False: # Warning: the FIT function has never been called
            w.warn("Warning: model not fitted, parameters are the initial ones.")

        # Solve for the optimal set of parameters
        mysolver = Solver(dt,self.fun,self.t[0],self.t[-1],self.pars[self.index:],self.pars[:self.index])

        # Solve for a dense grid, not for data points
        return mysolver.solve()

if __name__ == "__main__":
    from integrator import *
    from test import *
    import numpy as np
    from numpy.random import randn
    import matplotlib.pylab as plt

    def randomize(lst,scale=100.):
        """
        Randomize the data contained in lst (displacement along y axis)
        """
        for i in range(len(lst)):
            lst[i] += randn() / scale

        return lst

    def random_data(dt,fun,t0,T,y0,sel=10):
        """
        Create random data for the system of ODEs
        """

        mysolver = Solver(dt,fun,t0,T,y0)

        # Solve the system of ODEs
        t,y = mysolver.solve()

        t_stripped = []
        y_stripped = []

        # Select one out of SEL points from the solution
        for i in range(len(y)):
            if i % sel == 0:
                t_stripped.append(t[i])
                y_stripped.append(y[i])

        y_random = []

        # Randomize the stripped set of data points (mimiks experimental noise)
        for i in y_stripped:
            y_random.append(randomize(i))

        return t_stripped, y_random

    # Test function representing the system of ODEs
    test_function = rapid_equilibrium

    dt = 0.1

    t0 = 0
    T = 100
    y0 = [1,0,0]

    # Creates noisy data
    t,y = random_data(dt,test_function,t0,T,y0)

    y = np.asarray(y)

    # Initial guess for the parameters
    pars0 = [1,0.01]

    # Crate a ODEFit object: data points are the noisy solution
    myfit = ODEFit(test_function,t,y,y0,pars0)

    # Fit noisy data points
    pars = myfit.fit()[0] # Extracts only the real parameters and not initial conditions

    # Print parameters
    print(pars)

    # Solve the system of ODEs using optimal parameters
    tt,yy = myfit.solve_fitted()

    plt.plot(t,y,'x') # Plot noisy data
    plt.plot(tt,yy) # Plot fiting function
    plt.show()

"""
Copyright (C) 2015 iGEM Team EPF_Lausanne

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>



HumanBlaster.py -- blasts gRNAs against human genome

Requirements :
    - Install selenium on your computer:
        http://selenium-python.readthedocs.org/en/latest/installation.html
    - Save the gRNAs in a file named gRNAs like this : name:sequence
"""

from selenium import webdriver
from selenium.webdriver.common.keys import Keys
from selenium.webdriver.support.ui import Select
from selenium.webdriver.common.by import By
from selenium.webdriver.support.ui import WebDriverWait
from selenium.webdriver.support import expected_conditions as EC
import os

def skip(inp, nb):
    """
    Skip NB lines on INP file.
    """
    for _ in range(nb):
        inp.readline()

def seqReader(inp, initSkip):
    """
    Analyse the blast sequence.

    This function will display the percentage of similarity with a sequence in
    the human genome.

    You can change the sensibility of the script by changing the conditions
        cond = (int(k[6]) == 1 and int(k[7]) >= 23)
        cond = (int(k[7]) == 23 and int(k[6]) <= 1)
    """

    skip(inp, initSkip) # Skip initial lines

    line = inp.readline()

    name = str(line.split(":")[1][1:-1])

    skip(inp, 4) # Skip four lines

    l = inp.readline()

    war = 0
    while l != "" and l[0] != "#":
        k = l.split("\t")
        if len(k) < 7:
            return ""
        if (name[:2] == "CC"):
            cond = (int(k[6]) == 1 and int(k[7]) >= 23)
        else:
            cond = (int(k[7]) == 23 and int(k[6]) <= 1)

        if cond:
            print (name + ": " + k[1].split("|")[3] + " with " + k[2] + "% similarity, " + k[4] + " mismatches.")

        l = inp.readline()

    return l

with open("gRNAs", 'r') as input1, open("FASTA_gRNAs", 'w') as output:
    text = ""
    for i in input1:
        name = i.split(":")[0]
        seq = i.split(":")[1][:-1]
        text += ">" + name + "_TGG" + "\n" + seq + "TGG\n"
        text += ">" + name + "_GGG" + "\n" + seq + "GGG\n"
        text += ">" + name + "_AGG" + "\n" + seq + "AGG\n"
        text += ">" + name + "_CGG" + "\n" + seq + "CGG\n"
        output.write(">" + name + "_TGG" + "\n" + seq + "TGG\n")
        output.write(">" + name + "_GGG" + "\n" + seq + "GGG\n")
        output.write(">" + name + "_AGG" + "\n" + seq + "AGG\n")
        output.write(">" + name + "_CGG" + "\n" + seq + "CGG\n")

        text += ">" + "CCG_" + name + "\n" + "CCG" + seq + "\n"
        text += ">" + "CCC_" + name + "\n" + "CCC" + seq + "\n"
        text += ">" + "CCA_" + name + "\n" + "CCA" + seq + "\n"
        text += ">" + "CCT_" + name + "\n" + "CCT" + seq + "\n"
        output.write(">" + "CCG_" + name + "\n" + "CCG" + seq + "\n")
        output.write(">" + "CCC_" + name + "\n" + "CCC" + seq + "\n")
        output.write(">" + "CCA_" + name + "\n" + "CCA" + seq + "\n")
        output.write(">" + "CCT_" + name + "\n" + "CCT" + seq + "\n")



    driver = webdriver.Firefox() # Will open Firefox
    driver.get("http://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome") #Address
    assert "Nucleotide BLAST: Search nucleotide databases using a nucleotide query" in driver.title
    elem = driver.find_element_by_id("seq") #Selects the place where the sequences will be written
    elem.send_keys(text) #Writes them
    elem = driver.find_element_by_id("qtitle")
    elem.send_keys("Our Blast") #Change here if you want to have another name for the BLAST
    driver.find_element_by_id("Rhc").click() #Chooses human genome
    driver.find_element_by_id("b1").click() #Starts the BLAST
    #Will wait 500000 seconds or until the download page is ready
    try:
        element = WebDriverWait(driver, 500000).until(EC.presence_of_element_located((By.XPATH, "/html/body/div[1]/div[2]/div[3]/div[2]/form/table/tbody/tr[2]/td[1]/div[2]/a[5]")))
    finally:
        driver.close() #Will leave if it takes too long (13 hours max, should be enough)
    download = driver.find_element_by_id("showDownload").click() #Will start the download
    download = driver.find_element_by_id("hitText").click()

print("You need to save the file in this folder using the name resultat.txt")
continuation = input('When you are done press 1\n')

if continuation:
    with open("resultat.txt", 'r') as inp:
        l = seqReader(inp, 2)
        while (l != "" and l[0] == "#"):
            l = seqReader(inp, 1)