Demo entry 6345749

sdsd

   

Submitted by anonymous on Feb 06, 2017 at 12:47
Language: Python. Code size: 14.6 kB.

#! env/bin python
# Este programa define una clase Fasta con algunos metodos utiles para el
# analisis de este formato. Autora: Beatriz Beamud Aranguren

# ! usr/bin python

import numpy as np
import matplotlib.pyplot as plt
import re

class Fasta:
    def __init__(self, fich):
        self.fich = fich
        self.head, self.seq = self.leer_fasta()
        self.rev = self.reverse()
        self.complem = self.complementary()
        self.pauta_1, self.pauta_2, self.pauta_3, self.pauta_n1, \
        self.pauta_n2, self.pauta_n3 = self.pautas_lectura()

    # Para obtener la cabecera y la secuencia de un archivo fasta
    def leer_fasta(self):
        fichero = open(self.fich, 'r')
        head = fichero.readline()
        seq = ''
        for linea in fichero:
            # Quitamos los retornos de carro
            if linea[-1] == "\n":
                linea = linea[:len(linea) - 1]
            if linea[-1] ==  "\r":
                linea = linea[:len(linea) - 1]
            seq = seq + linea
        # Pasamos todos los nt a mayuscula
        seq.upper()
        fichero.close()
        return head, seq

    # Obtenemos la secuencia reversa
    def reverse(self):
        longitud = len(self.seq)
        rev = ''
        # Cuenta regresiva
        for i in range(longitud - 1, 0, -1):
            base = self.seq[i]
            rev += base
        return rev

    # Obtenemos la secuencia complementaria
    def complementary(self):
        complem = []
        for i in range(0, len(self.rev)):
            if self.rev[i] == 'A':
                complem.append('T')
            if self.rev[i] == 'T':
                complem.append('A')
            if self.rev[i] == 'C':
                complem.append('G')
            if self.rev[i] == 'G':
                complem.append('C')
        complem = ''.join(complem)
        return complem

    # Obtenemos las 6 pautas de lectura, positivas y negativas, de una seq
    def pautas_lectura(self):
        pauta_1 = []
        pauta_2 = []
        pauta_3 = []
        pauta_n1 = []
        pauta_n2 = []
        pauta_n3 = []
        for i in range(0, len(self.seq), 3):
            start1 = i
            start2 = i + 1
            start3 = i + 2
            pauta_1.append(self.seq[start1:start1 + 3])
            pauta_2.append(self.seq[start2:start2 + 3])
            pauta_3.append(self.seq[start3:start3 + 3])
        for i in range(0, len(self.complem), 3):
            start1 = i
            start2 = i + 1
            start3 = i + 2
            pauta_n1.append(self.complem[start1:start1 + 3])
            pauta_n2.append(self.complem[start2:start2 + 3])
            pauta_n3.append(self.complem[start3:start3 + 3])
        return pauta_1, pauta_2, pauta_3, pauta_n1, pauta_n2, pauta_n3

    # Calcula la freq de codones para las 6 pautas de lectura
    def frecTri(self):
        freq_pautas = []
        bases = ["T", "C", "G", "A"]
        # Con estas 3 iteraciones se crean todas las posibilidades de codones
        codones = [a + b + c for a in bases for b in bases for c in bases]
        for i in [self.pauta_1, self.pauta_2, self.pauta_3, self.pauta_n1,
                  self.pauta_n2, self.pauta_n3]:
            freq_pauta = {}
            for codon in codones:
                count_i = i.count(codon)
                # Solo se anyaden los que esten presentes en la cadena de DNA
                if count_i > 0:
                    freq_pauta[codon] = count_i
            freq_pautas.append(freq_pauta)
        return freq_pautas

    # Funcion que obtiene la base mas repetida consecutivamente
    def maxRep(self):
        bases = ['A', 'C', 'G', 'T']
        maximo = 0
        for i in bases:
            # Lista vacia donde se anyadira las diferentes veces que aparece
            # la letra conseq para luego tomar el valor maximo
            conseq = [0]
            # Suponemos que la letra esta al menos una vez
            contador = 1
            for c in range(1, len(self.seq)):
                # si esta la base que nos interesa se itera sobre esta
                # para ver si esta consecutivamente
                if self.seq[c] == i:
                    if self.seq[c - 1] == self.seq[c]:
                        contador = contador + 1
                        conseq.append(contador)
                    else:
                        # si encuentra otra base, el contador volvera a 1.
                        contador = 1
                        # se anyade este valor a la lista por si la base
                        # esta solo 1 vez y nunca entrara en el if anterior
                        conseq.append(contador)
                else:
                    contador = 0
                    conseq.append(contador)
            # Si el valor maximo de las repeticiones de cada base es mayor
            # que la anterior, se reeemplaza el valor maximo
            if max(conseq) > maximo:
                maximo = max(conseq)
                base = i
                posicion = conseq.index(max(conseq))
        print "La base que mas veces aparece es la %s y se encuentra %s " \
              "veces consecutivas. Se encuentra en la posicion %i - %i." \
              %(base, maximo, (posicion-maximo), posicion)

    # Esta funcion obtiene el triplete mas repetido consecutivamente
    def maxRepTri(self):
        pauta = 0
        for i in [self.pauta_1, self.pauta_2, self.pauta_3, self.pauta_n1,
                  self.pauta_n2, self.pauta_n3]:
            contador = 1
            maximo = 1
            posicion = 0
            for c in range(1, len(i)):
                posicion += 3
                if i[c] == i[c -1 ]:
                    contador += 1
                    if contador > maximo:
                        posicion_final = (c * 3 + 3)
                        maximo = contador
                        triplete = i[c]
                else:
                    contador = 1
            pauta += 1
            print "El triplete que mas veces aparece en la pauta %s es " \
                  "%s y se encuentra %s veces consecutivas. Se encuentra en " \
                  "la posicion %i - %i." % (pauta, triplete, maximo,
                                            (posicion_final - maximo * 3),
                                            posicion_final)
        print "NOTA: Las pautas 4, 5 y 6 corresponden a las pautas negativas"

    # Esta funcion ordena las frecuencias de un diccionario
    def ordenar_freq(self, dicc):
        print ("Introduce si quiere ordenar de mayor a menor "
                           "(M) o al contrario (m)")
        orden = raw_input("Si no quiere ordenar cualquier otra tecla:\n")
        # Siendo v:value, k:key
        lista = [(v, k) for k, v in dicc.items()]
        if orden in ['m','M']:
            lista.sort()
        if orden == 'M':
            # Para ordenar de mayor a menor
            lista.reverse()
        # Tenemos una lista de tuplas. Lo ponemos en el orden correcto.
        codones = [k for v, k in lista]
        frecuencias = [v for v, k in lista]
        return codones, frecuencias

    # Esta funcion genera un histograma de codones vs frecuencia
    def hist_tripl(self, numero_p, codones, frecuencias):
        # Obtengo el numero de valores del eje x
        x = np.arange(len(codones))
        plt.bar(x, frecuencias, align='center')
        # Asignamos el nombre del codon a cada barra
        plt.xticks(x, codones)
        plt.title(("Histograma pauta", numero_p))
        plt.xlabel("Codones")
        plt.ylabel("Frecuencia aparicion")
        plt.show()

    # Esta funcion construye una matriz en la que se comparan dos secuencias
    def comparar_secuencias(self, seq2):
        # Nos interesa obtener una matriz cuadrada con la dimension de la
        #  secuencia de mayor longitud.
        if len(self.seq) >= len(seq2.seq):
            len_max = len(self.seq)
            dif_len = len_max - len(seq2.seq)
            # A la seq mas pequenya se le ponen tantos '-' como bases falten
            seq2.seq = seq2.seq + ('-' * dif_len)
            # Con np.empty creamos una matriz nula de [dimesion, dimension]
            resultado = np.empty([len(self.seq), len(self.seq)])
        else:
            resultado = np.empty([len(seq2.seq), len(seq2.seq)])
            len_max = len(seq2.seq)
            dif_len = len_max - len(self.seq)
            self.seq = self.seq + ('-' * dif_len)
        # por cada col y fila de la matriz si la posicion de la
        # columna(seq1)= posicion de la fila(seq2), se anyade un 1
        for c in range(len(resultado)):
            for r in range(len(resultado)):
                if self.seq[c] == seq2.seq[r]:
                    resultado[c][r] = 1
        return resultado

    # Esta funcion nos permite realizar un dot plot a partir de una matriz
    def dot_plot(self, matriz):
        # creo la figura
        plt.figure()
        #  represento la matriz
        plt.matshow(matriz)
        plt.title("Dotplot")
        plt.show()

    # Esta funcion elimina los caracteres no deseados de la seq mediante re
    def limpiar_seq_re(self):
        # Eliminamos caracteres en blanco
        seq_limpia = re.sub('\s*', "", self.seq)
        # Eliminamos digitos
        seq_limpia = re.sub('\d*', "", self.seq)
        return seq_limpia

    def extraer_tri_re(self):
        pauta = 1
        seq_p1 = self.seq
        seq_p2 = self.seq[1:]
        seq_p3 = self.seq[2:]
        seq_n1 = self.complem
        seq_n2 = self.complem[1:]
        seq_n3 = self.complem[2:]
        tripletes_totales = set()
        print "\n Los tripletes de cada pauta son los siguientes:\n"
        for i in [seq_p1, seq_p2, seq_p3, seq_n1, seq_n2, seq_n3]:
            tripletes = re.findall("[ACGT][ACGT][ACGT]", i)
            print 'Pauta: ', pauta
            print tripletes
            for t in tripletes:
                tripletes_totales.add(t)
            pauta += 1
        print "\nNOTA: Las pautas 4, 5 y 6 corresponden a las pautas " \
              "negativas\n"
        print "\nTripletes totales unicos:\n"
        print tripletes_totales

    def pos_tri_re(self):
        pauta = 1
        seq_p1 = self.seq
        seq_p2 = self.seq[1:]
        seq_p3 = self.seq[2:]
        seq_n1 = self.complem
        seq_n2 = self.complem[1:]
        seq_n3 = self.complem[2:]
        tripletes_totales = set()
        print "\n Los tripletes de cada pauta son los siguientes:\n"
        for i in [seq_p1, seq_p2, seq_p3, seq_n1, seq_n2, seq_n3]:
            print 'Pauta: ', pauta
            diccionario = {}
            tripletes = re.findall("[ACGT][ACGT][ACGT]", i)
            pos_tripletes = re.finditer('[ACGT][ACGT][ACGT]', i)
            indices = []
            for i in pos_tripletes:
                it = i.start()
                indices.append(it)
            for j in range(0, len(tripletes)):
                clave = tripletes[j]
                valor = indices[j]
                if clave in diccionario:
                    diccionario[clave] = str(diccionario[clave]) + "," + str(valor)
                else:
                    diccionario[clave] = valor
            print diccionario


def menu_usuario():
    print "Que desea realizar?"
    print '----------------------------------------------------------'
    print " a.) Obtener la frecuencia absoluta de tripletes"
    print " b.) Obtener la base mas repetida y su posicion"
    print " c.) Obtener el triplete mas repetido y su posicion"
    print " d.) Obtener un histograma con la frecuencia de tripletes"
    print " e.) Obter un dotplot con una secuencia o con varias"
    print " f.) Depurar secuencias mediante re"
    print " g.) Extraer todos los tripletes de todas las pautas mediante re"
    print " h.) Extraer los tripletes y sus posiciones mediante re"
    print " z.) Salir"
    print '----------------------------------------------------------'
    respuesta = raw_input()
    return respuesta


def main():
    fich = raw_input("\nIntroduce el fichero FASTA que quiere analizar:\n")
    fasta = Fasta(fich)
    respuesta = menu_usuario()

    if respuesta == 'a':
        freq_pautas = fasta.frecTri()
        for i in range(len(freq_pautas)):
            print 'Pauta ', i+1
            print freq_pautas[i]

    if respuesta == 'b':
        fasta.maxRep()

    if respuesta == 'c':
        fasta.maxRepTri()

    if respuesta == 'd':
        freq_pautas = fasta.frecTri()
        i=1
        for pauta in freq_pautas:
            codones, frecuencias = fasta.ordenar_freq(pauta)
            hist_pauta = fasta.hist_tripl(i, codones, frecuencias)
            seguir = raw_input("Quieres ver la siguiente pauta de lectura"
                               " (s/n)?\n")
            if seguir in ['n', 'N']:
                break
            i += 1

    if respuesta == 'e':
        opcion = raw_input("Quieres comparar con la misma secuencia (s)"
                           " o con otra(n):\n")
        if opcion in ['s', 'S']:
            fasta2=Fasta(fich)
        else:
            seq2 = raw_input("Introduce la secuencia con la que la quieres "
                         "comparar:\n")
            fasta2 = Fasta(seq2)
        matriz = fasta.comparar_secuencias(fasta2)
        fasta.dot_plot(matriz)

    if respuesta == 'f':
        seq_limpia =fasta.limpiar_seq_re()
        print "La secuencia quedaria asi:\n", seq_limpia
        guardar = raw_input("Desea guardarla (s/n)?\n")
        if guardar == 's':
            output = raw_input("Introduce el nombre del archivo de salida:\n")
            output = open(output + '.fa', 'w')
            # Escribimos cabecera, sin depurar.
            output.write(fasta.head)
            for line in seq_limpia:
                output.write(line)
            output.close()

    if respuesta == 'g':
        fasta.extraer_tri_re()

    if respuesta =='h':
        print "Se imprimira un diccionario por cada pauta siendo la clave" \
              " el triplete y el valor todas las posiones donde esta"
        fasta.pos_tri_re()

    if respuesta == 'z':
        exit()

    repetir = raw_input("\nDesea realizar otro analisis (s/n)?\n")
    if repetir in ['s', 'S']:
        main()

if __name__ == '__main__':
    main()

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).