Python_Biotec - 2_Bio-checkpoint


Cadena Inversa


In [1]:
#Función complementaria
def inversa(dna):
    for nuc in dna:
        if nuc not in ['A', 'T', 'C', 'G']:
            print 'No es DNA. Revise su cadena'
            return "" 
    
    inversa = dna[::-1]
    
    return inversa      

dna = raw_input('Escriba una cadena de DNA: ')
dna = dna.upper()

dna_inversa = inversa(dna)

if inversa:
    print 'La cadena inversa es: ' + dna_inversa


Escriba una cadena de DNA: TAGAGTAC
La cadena inversa es: CATGAGAT

Cadena Complementaria


In [2]:
#Función complementaria
def complementaria(dna):
    
    comp_dic = {'A':'T','T':'A','C':'G','G':'C'}    
    
    comp = ""
    
    for nuc in dna:
        if nuc not in ['A', 'T', 'C', 'G']:
            print 'No es DNA. Revise su cadena'
            return ""
        else: comp += comp_dic[nuc]
           
    return comp



dna = raw_input('Escriba una cadena de DNA: ')
dna = dna.upper()

dna_complementaria = complementaria(dna)

if dna_complementaria:
    print 'La cadena complementaria es: ' + dna_complementaria


Escriba una cadena de DNA: CCTATCGA
La cadena complementaria es: GGATAGCT

Inversa Complementaria


In [3]:
dna = raw_input('Escriba una cadena de DNA: ')
dna = dna.upper()

dna_inv_comp = inversa(complementaria(dna))

if dna_inv_comp:
    print 'La cadena inversa complementaria es: ' + dna_inv_comp


Escriba una cadena de DNA: gatatttac
La cadena inversa complementaria es: GTAAATATC

DNA to RNA


In [4]:
def dna2rna(dna):
    dna=dna.upper()
    rna=dna.replace('T', 'U')
    return rna 
    
    
    
dna=raw_input("¿Cómo es tu DNA?  ")
rna = dna2rna(dna)
print 'Tu RNA es: ' + rna


¿Cómo es tu DNA?  GATTACA
Tu RNA es: GAUUACA

Cuenta AA


In [9]:
def fasta2sec(filename):
    file = open (filename)
    if not file:
        print "El archivo no existe"
        return ""
    secuencia = ''
    #Iniciamos el bucle de lectura
    for linea in file:
        #Asociamos a la variable "linea" cada línea del archivo, pero sin el salto de línea
        linea = linea.strip('\n')
        
        #Comprobamos que la línea no está vacía, ni empieza por '>' (1ª linea de info)
        if linea != '' and linea[0] != '>' :
            #Añadimos a la variable "secuencia" cada iteración de linea
          secuencia = secuencia + linea
        
        
    return secuencia

filename = raw_input('Escribe la ruta de tu fichero: ')

secuencia = fasta2sec(filename)

print "Tu proteína tiene " + str(len(secuencia)) + " aminoácidos"


Escribe la ruta de tu fichero: dna.txt
Tu proteína tiene 290 aminoácidos

Tipos AA


In [6]:
def contar_tipos(secuencia):
    
    hydrophilic = ['C', 'F', 'I', 'L', 'M', 'V', 'W']
    hydrophobic = ['D', 'E', 'G', 'K', 'N', 'P', 'Q', 'R', 'S', 'T']
    neutral = ['A', 'H', 'Y']
    
    contadores = {'hphi': 0, 'hpho': 0, 'neu'  : 0}

    #Iteramos la secuencia por elemento
    for aa in secuencia:

        #Comprobamos dónde está cada aa, y sumamos 1 al contador correspondiente a su tipo
        if aa in hydrophilic:
            contadores['hphi'] +=1
        elif aa in hydrophobic:
            contadores['hpho'] +=1
        elif aa in neutral:
            contadores['neu'] +=1
        #Añadimos "else" por si acaso
        else: 
            print 'Algún dato no es AA. Revisa el fichero'
            #Preguntamos si el usuario quiere continuar
            pregunta = raw_input('¿Quieres continuar de todos modos? [s/n]')
            #Si la respuesta es sí...
            if pregunta == 'y' or pregunta == 'Y':
                return contadores 
            #Si la respuesta es no...
            elif pregunta =='n' or pregunta == 'N':
                return contadores
            #Si la respuesta es otra...
            return contadores
        
    return contadores
        
contadores = contar_tipos(secuencia)        
                
#Imprimimos la respuesta con los datos
print 'La cantidad de aa de cada tipo es: '
print '%.2f %% de hidrofílicos' %(float(contadores['hphi'])/len(secuencia)*100)
print '%.2f %% de hidrofóbicos' %(float(contadores['hpho'])/len(secuencia)*100)
print '%.2f %% de neutros' %(float(contadores['neu'])/len(secuencia)*100)


La cantidad de aa de cada tipo es: 
26.90 % de hidrofílicos
45.52 % de hidrofóbicos
27.59 % de neutros

Visualizar Tipos de AA


In [10]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

index = range(len(contadores))

barlist = plt.bar(index, contadores.values(), align='center')

# Un poco de color
barlist[0].set_color('#999999')
barlist[1].set_color('#ff3333')
barlist[2].set_color('#3333ff')

plt.xticks(index, contadores.keys())
plt.show()


Comparar DNA


In [102]:
def identidad(sec1, sec2):
    count = 0
    dif = []
    for i in range(len(sec1)):
        if sec1[i] == sec2[i]:
            count += 1
            dif.append(1)
        else:
            dif.append(0)
            
    # Devolvemos % de similitud y lista de diferencias)
    return (count*100./len(sec1), dif)
    
pan = fasta2sec('pan.fasta')
sapiens = fasta2sec('sapiens.fasta')

pan = pan[216-1:659]
sapiens = sapiens[196-1:639]

print pan[0:3]
print sapiens[0:3]
    
diffs =  identidad(pan,sapiens)

print 'Porcentaje de similitud %.2f %%' % diffs[0]


ATG
ATG
Porcentaje de similitud 99.32 %

Visualizar diferencias de DNA


In [108]:
from numpy import array

difs_array = array( [diffs[1]] )

plt.imshow(difs_array, aspect = 50)

ax = plt.gca()
ax.yaxis.set_visible(False)


##Poner Titulo

Definiciones


In [111]:
# Codigo genetico
gencode = {
  'ATA':'I', 'ATC':'I', 'ATT':'I', 'ATG':'M',
  'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACT':'T',
  'AAC':'N', 'AAT':'N', 'AAA':'K', 'AAG':'K',
  'AGC':'S', 'AGT':'S', 'AGA':'R', 'AGG':'R',
  'CTA':'L', 'CTC':'L', 'CTG':'L', 'CTT':'L',
  'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCT':'P',
  'CAC':'H', 'CAT':'H', 'CAA':'Q', 'CAG':'Q',
  'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGT':'R',
  'GTA':'V', 'GTC':'V', 'GTG':'V', 'GTT':'V',
  'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCT':'A',
  'GAC':'D', 'GAT':'D', 'GAA':'E', 'GAG':'E',
  'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGT':'G',
  'TCA':'S', 'TCC':'S', 'TCG':'S', 'TCT':'S',
  'TTC':'F', 'TTT':'F', 'TTA':'L', 'TTG':'L',
  'TAC':'Y', 'TAT':'Y', 'TAA':'_', 'TAG':'_',
  'TGC':'C', 'TGT':'C', 'TGA':'_', 'TGG':'W'}


# Lista de aminoacidos
aminoacidos = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 
               'I', 'K', 'L', 'M', 'N', 'P', 'Q',  
               'R', 'S', 'T', 'V', 'W', 'Y', '_']
               
# Diccionario que para cada aminoacido
# da como valor la lista de codones que lo codifican
# separados por el caracter ":"
aa2codon={'A':'GCA:GCC:GCG:GCT', 'C':'TGT:TGC', 'D':'GAT:GAC', 
          'E':'GAG:GAA', 'F':'TTT:TTC', 'G':'GGT:GGG:GGA:GGC' , 
          'H':'CAT:CAC', 'I':'ATC:ATA:ATT', 'K':'AAG:AAA', 
          'L':'CTC:CTG:CTA:CTT:TTA:TTG', 'M':'ATG', 'N':'AAC:AAT', 
          'P':'CCT:CCG:CCA:CCC', 'Q':'CAA:CAG', 
          'R':'AGG:AGA:CGA:CGC:CGG:CGT', 
          'S':'AGC:AGT:TCG:TCA:TCC:TCT', 
          'T':'ACC:ACA:ACG:ACT', 'V':'GTA:GTC:GTG:GTT',
          'W':'TGG', 'Y':'TAT:TAC', '_':'TAG:TAA:TGA'}

In [1]:
#Abrimos el archivo fasta
file = open('SP.fasta')

secuencia = fasta2sec('SP.fasta')


#Archivo de exones
file = open('CDScoords.txt')

#Alternativa para customizar el programa. Esto lee el nombre del primer gen en
#vez de que se lo demos en el código
'''for linea in file:
    pos = linea.split()
    a = pos[0]
    break'''

a = 'SPAC212.06c'

nueva_secuencia = ''
for linea in file:
    pos = linea.split()
    if pos[0] == a:
        nueva_secuencia += secuencia[int(pos[1])-1:int(pos[2])]
    else:
        print 'Para el gen %s: ' %a #Y ahora la función que hace todo
        codon_counter={'A':[0,0,0,0], 'C':[0,0], 'D':[0,0], 
          'E':[0,0], 'F':[0,0], 'G':[0,0,0,0] , 
          'H':[0,0], 'I':[0,0,0], 'K':[0,0], 
          'L':[0,0,0,0,0,0], 'M':[0], 'N':[0,0], 
          'P':[0,0,0,0], 'Q':[0,0], 
          'R':[0,0,0,0,0,0], 
          'S':[0,0,0,0,0,0], 
          'T':[0,0,0,0], 'V':[0,0,0,0],
          'W':[0], 'Y':[0,0], '_':[0,0,0]}
        #print nueva_secuencia
        for i in range(0,len(nueva_secuencia),3) :
            codon = nueva_secuencia[i:i+3]
            aa = gencode[codon]
            posibles_codones = aa2codon[aa].split(':')
            #print posibles_codones
            for x in range(0,len(posibles_codones)):
                #print aa, codon, x
                if codon == posibles_codones[x]:
                    codon_counter[aa][x] += 1
                    break
        #print codon_counter
        for y in aminoacidos:
            print 'Aminoácido: %s' %y
            for z in range(0,len(aa2codon[y].split(':'))):
                total = 0
                for i in codon_counter[y]:
                    total += i
                #print total
                if total == 0:
                    print '  %s %.2f%%' %(aa2codon[y].split(':')[z],(codon_counter[y][z]))
                else: print '  %s %.2f%%' %(aa2codon[y].split(':')[z],(codon_counter[y][z]*100./total))
                #print '  %s %s' %(codon_actual,porcentaje)
        nueva_secuencia = secuencia[int(pos[1])-1:int(pos[2])]
        a = pos[0]


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-704233b63f26> in <module>()
      2 file = open('SP.fasta')
      3 
----> 4 secuencia = fasta2sec('SP.fasta')
      5 
      6 

NameError: name 'fasta2sec' is not defined

Alineamiento de Proteínas


In [14]:
import nwalign as nw

D = fasta2sec('Q9V429.fasta')
H = fasta2sec('P10599.fasta')

acounter = 0
bcounter = 0

a, b = nw.global_align(D, H, matrix='BLOSUM62.txt')

for aa in a:
    if aa != '-':
        acounter += 1 
        
for aa in b:
    if aa != '-':
        bcounter += 1 
        
print ((acounter*100./len(a))+(bcounter*100./len(b))) / 2


84.2307692308

In [16]:
b


Out[16]:
'M-V--KQIESKT-A-F--QEA-LD-----AAGDKLVVVDF-SATWCGPCKMIKP-FFHSLS-EKYS-NVIFLEVDVDDCQDVASE-CEVKCMPTF-QFFKKG-QKVGE-FSGAN-KEKLE-ATI--NELV'