# -*- coding: utf-8 -*- """ Created on Wed Apr 13 12:23:49 2016 @author: rodri """ #---------------------------------------------------------------- #%%---------------- JUGANDO CON GENOMAS-------------------------- #---------------------------------------------------------------- #Este trocito de código nos permite leer un genoma de una URL #En particular, el genoma de E coli import urllib2 f=urllib2.urlopen("http://vis.usal.es/rodrigo/documentos/bie/E-coli.txt", "r") ec=f.readlines() #toma todas las líneas del archivo y las mete en una lista ec=ec[0] #el fichero de genoma sólo tiene una línea import string ec=string.replace(ec, "\n", "") #para reemplazar retornos de línea (\n) f.close() #%% PREGUNTA #ec es una cadena con el genoma de Escherichia coli #¿Cuántos nucleótidos exactamente tiene el genoma de E coli? #¿Cuáles son sun 10 primeros nucleótidos? #%%---------------------- COMPLEMENTARIA ---------------------- #Retorna la secuencia complementaria de seq #Adenina (A) y timina (T) son complementarias, #al igual que guanina (G) y citosina (C) def complement(seq): comp=[] us=seq.upper() for s in us: if(s=='A'): comp.append('T') elif(s=='T'): comp.append('A') elif(s=='G'): comp.append('C') elif(s=='C'): comp.append('G') return('%s' % ''.join(map(str, comp))) #%% print complement("ACGTTCGT") print complement(ec[:10]) #%%---------------------- INVERSA ---------------------- #Retorna la secuencia inversa de seq def reverse(seq): rev=[] for i in range(len(seq)-1,-1,-1): rev.append(seq[i]) return('%s' % ''.join(map(str, rev))) reverse("ACGTTCGT") #%%----------- INVERSA COMPLEMENTARIA ----------- #Retorna la secuencia inversa y complementaria de seq def revcomp(seq): return complement(reverse(seq)) print revcomp("ACGTTCGT") #Una cadena de ADN es doble. Generalmente sólo se indica una cadena #porque a partir de ella podemos obtener la otra, pues contiene los nucleótidos #complementarios y se lee en el orden inverso #%% PREGUNTA #Utilizando la función revcomp, ¿cuál sería la cadena complementaria de los # 20 nucleótidos que están justo en el centro del genoma de E coli? #%% ---------------- FRECUENCIA --------------- #Retorna un diccionario con la frecuencia de aparición de cada nucleótido en seq def freq(seq): f={'A':0,'C':0,'G':0,'T':0} for x in seq: if(x in f.keys()): f[x]=f[x]+1 for k in f.keys(): f[k]=(float)(f[k])/len(seq) return f print freq("TAAAGACTGCCGAGAGGCCAACACGAGTGCTAGAACGAGGGGCGTAAACGCGGGTCCGAT") print freq("ACGTTCGT") #%%Calculamos las frecuencias en E coli freqEC=freq(ec) print freqEC #%% PREGUNTA #¿Tienen lógica los valores de frecuencia de aparición? #------------------------------------------- #%%------------ CONTENIDO GC --------------- #------------------------------------------- #Si la distribución de los nucleótidos fuera aleatoria, cada uno tendría una #frecuencia de alrededor del 25%. Esto es justo lo que pasaba con E. coli: #Sin embargo, algunas especies tienen un contenido más alto de unos nucleótidos # que de otros. #Este es el caso de las bacterias, que se clasifican en grandes grupos según su #contenido en G y C (las frecuencias sumadas de ambos nucleótidos) #%% --------------- STREPTOMYCES COELICOLOR ---------------- #La familia de actinobacterias tienen genomas con alto contenido en GC. #El género Streptomyces es un ejemplo, y en particular Streptomyces coelicolor, #utilizada en la producción de antibióticos. import urllib2 f=urllib2.urlopen("http://vis.usal.es/rodrigo/documentos/bie/GCF_000203835.1_ASM20383v1_genomic.fna", "r") sc=f.readlines() #toma todas las líneas del archivo y las mete en una lista f.close() #el fichero está en formato FASTA, la primera línea contiene una descripción sc=sc[1:] sc="".join(sc) #en sc tenemos el genoma completo de S coelicolor #%% PREGUNTA: Cuál es el contenido en GC de S. coelicolor? # (en porcentaje, hasta el 5º decimal) #%% --------------- PLASMODIUM FALCIPARUM ---------------- #Otro ejemplo es el del parásito Plasmodium falciparum (causante de la malaria) import urllib2 f=urllib2.urlopen("http://vis.usal.es/rodrigo/documentos/bie/Plasmodium_falciparum.fa", "r") pf=f.readlines() #toma todas las líneas del archivo y las mete en una lista #el fichero está en formato FASTA, la primera línea contiene una descripción pf=pf[1:] pf="".join(pf) import string pf=string.replace(pf, "\n", "") #pequeña modificación para reemplazar retornos de línea (\n) f.close() #%% PREGUNTA: Cuál es el contenido en GC del P. falciparum? #Coincide con lo que dice la wikipedia en el artículo sobre su contenido GC? #https://en.wikipedia.org/wiki/GC-content #(apartado "Application in systematics") #------------------------------------------- #%% ----- ORIGENES DE REPLICACIÓN ---------- #------------------------------------------- #En ciertas bacterias, el origen de replicación estará donde haya #un cambio en la frecuencia de G y C #%% -------------------- SKEW --------------------- #Retornar una cadena numérica que vaya acumulando las G y restando las C def skew(seq): count=0 sk=[] sk.append(count) for i in range(len(seq)): if(seq[i]=='G'): count=count+1 elif(seq[i]=='C'): count=count-1 sk.append(count) return(sk) #%% Pruebas seq="CATGGGCATCGGCCATACGCC" sk=skew(seq) sk #%% Min skew #Retornar la posición o posiciones donde el skew es mínimo def minimumSkew(seq): sk=skew(seq) mins=[] msk=min(sk) for i in range(len(sk)): if(sk[i]==msk): mins.append(i) return(mins) #%% Pruebas seq="TAAAGACTGCCGAGAGGCCAACACGAGTGCTAGAACGAGGGGCGTAAACGCGGGTCCGAT" sk=skew(seq) minimumSkew(seq) #%% Gráficas #Mostrando de manera gráfica el skew import matplotlib.pyplot as plt plt.plot(sk, label='skew') plt.ylabel('skew') plt.xlabel('position') plt.show() #%% Con E coli sk=skew(ec) ms=minimumSkew(ec) plt.plot(sk, label='skew') plt.ylabel('skew') plt.xlabel('position') plt.title="Escherichia coli" plt.show() oricEC=ec[ms[0]:(ms[0]+500)] #%% Con Vibrio cholerae #Cargamos su genoma import urllib2 f=urllib2.urlopen("http://vis.usal.es/rodrigo/documentos/bie/Vibrio_cholerae.txt", "r") vc=f.readlines() #toma todas las líneas del archivo y las mete en una lista vc=vc[0] #el fichero de genoma sólo tiene una línea import string vc=string.replace(vc, "\n", "") #para reemplazar retornos de línea (\n) f.close() #%%Calculamos el skew sk=skew(vc) print minimumSkew(vc) plt.plot(sk, label='skew') plt.ylabel('skew') plt.xlabel('position') plt.title="Vibrio cholerae" plt.show() #%% Ahora con Thermotoga! import urllib2 f=urllib2.urlopen("http://vis.usal.es/rodrigo/documentos/bioinfo/avanzada/genomas/Thermotoga-petrophila.txt", "r") tg=f.readlines() #toma todas las líneas del archivo y las mete en una lista tg=tg[0] #el fichero de genoma sólo tiene una línea #%% sk=skew(tg) print minimumSkew(tg) plt.plot(sk, label='skew') plt.title="Thermotoga petrophila" plt.ylabel('skew') plt.xlabel('position') plt.show()