GNU GENERAL PUBLIC LICENSE
Version 3, 29 June 2007Ejecutar secuencialmente cada celda del notebook, cambiando segun corresponda los nombres de archivo.
In [1]:
from Bio import SeqIO
# para leer fasta
# si da error "No module named Bio",
# install Biopython, ej: pip install Biopython o sudo apt-get install python-biopython
In [2]:
#lectura de secuencias objetivo
records=list(SeqIO.parse("DQB 9F.fa", "fasta"))
#Aqui deberia prever la lectura de los leidos del cromatograma en R
In [3]:
records[0]
Out[3]:
In [5]:
#convierto las secuencias a strings y descripciones
seq=[str(r.seq) for r in records]
desc=[str(r.description).strip() for r in records]
print desc[0]
print seq[0]
#convierto las secuencias en strings - OJO QUE DEJO EL ID
In [6]:
# Funciones de conversión de formato utilizadas en el código principal.
def seq_a_lista(s):
#recibe un string y devuelve una lista
# proceso las [A/G] eliminando los [] y uniendo los dos valores en torno a / como listas
# TODO: ver que pasa con los [./.]. De movida debe quedar .. ya que no estoy usando sets, que eliminarian la duplicación
s=s.replace('[','').replace(']','')
sl= [x for x in s]
for i,x in enumerate(sl):
if x=='/':
sl[i-1]+=sl[i+1]
del sl[i:i+2]
return sl
def lista_a_seq(l):
'''recibe una lista y devuelve un string con corchetes y barra en los polimorfismos
'''
ll=[x if len(x)<2 else '[' + x[0] + '/' + x[1] + ']' for x in l] #'/'.join(x)
return ''.join(ll)
def seq_a_string(s):
'''recibe un string con corchetes, lo transforma en lista y lo traduce a IUPAC
'''
return lista_a_iupac(seq_a_lista(s))
def guion_a_vacio(l): # creo que ya no se usa, ver. recibe un string o una lista y convierte guiones en vacíos
return [ x if (x<>'-') else '' for x in la]
def vacio_a_guion(l):
'''recibe una lista y convierte los vacíos en guiones -
si recibe un string lo devuelve tal cual
'''
return [ x if x else '-' for x in l]
def unir_haplotipos(la,lb):
'''recibe dos listas de haplotipos de referencia
maneja secuencias con guiones y las convierte a listas con vacíos
devuelve una lista de la union de las dos listas elemento a elemento como conjunto (sin duplicados)
'''
la=[ x if (x<>'-') else '' for x in la]
lb=[ x if (x<>'-') else '' for x in lb]
return [''.join(list(set(la[i]+lb[i]))) for i in range(len(la))]
def ref_haplo(la):
'''recibe un haplotipo de referencia, con los guiones convertidos en vacíos
completa los guiones - con los datos del haplotipo de referencia
'''
return [ refseq[i] if not x or x=='-' else x for i,x in enumerate(la)]
#return [ x if (x in 'ACGT') else refseq[i] for i,x in enumerate(la)] ## no son equivalentes, y este funciona siempre? Ver
In [7]:
print len(seq)
In [8]:
for s in seq: #imprimo separado por espacios para ver
print ' '.join(seq_a_lista(s))
In [9]:
#lectura del alineamiento
f=open('alineamiento DQB1.txt')
lineas=f.readlines()
# Queda en el diccionario d hasheado por id
d={}
ref=''
for i,l in enumerate(lineas):
if (len(l)>20 and l[0]==' ' and l[1]!=' '): # uso un patron para identificar las lineas con datos
id= l[:20].strip() #las corto por ancho fijo del id
#print i, id
sseq= l[21:].replace(' ','').rstrip() # elimino espacis
#print id
#print seq
if id in d.keys(): # y las guardo (concatenadas) en un diccionario indexado por id
d[id]+=sseq
else:
d[id]=sseq
if not ref: ref=id # me estoy quedando con la primera secuencia asumiendo que es la de referencia
d[id]=d[id].replace('*','-')
f.close()
In [10]:
# no hace falta, pero pongo la secuencia de referencia en esta variable (despues reorganizo)
ref
refseq = d[ref]
print ref, refseq
In [12]:
# Convierto las secuencias en listas y las guardo en seqs
seqs=[]
for s in seq:
seqs.append(seq_a_lista(s))
print seqs
In [13]:
print seqs[0]
In [14]:
#Verifico la cantidad de secuencias de referencia y sus identificadores
print len(d), d.keys()
In [15]:
#Generación de todos los haplotipos posibles, unidos de a pares
#Guardo en el log
haplos=[]
hh=d.keys()
f=open('haplos.log','w')
def log(l): #aca defino la funcion para generar los logs
for x in l:
f.write(str(x))
f.write('\n')
for x in hh:
for y in hh:
log((x,y))#print x,y
log(('x :', d[x]))#print 'x:', d[x]
log(('y :', d[y]))#print 'y:', d[y]
log(('xr:',ref_haplo(d[x])))
log(('yr:',ref_haplo(d[y])))
dd= unir_haplotipos(ref_haplo(d[x]), ref_haplo(d[y]))
log(('z0:', lista_a_seq((dd))))
log(('z1:', lista_a_seq(vacio_a_guion(dd))))
log(('z2:', lista_a_seq(dd)))
log(('z3:', dd))
haplos.append([[x,y],dd])
f.close()
* lectura, par de haplos, ranking, diferencias
In [16]:
#Visualizo cuantas combinaciones se generaron, cuantos haplotipos de referencia hay,
#y cuantos deberian haberse generado
print len(haplos), len(hh), len(hh)*len(hh)
#Y muestro un ejemplo
print hh[11], haplos[11]
In [17]:
print len(seqs)
In [18]:
# Visualizo las listas como secuencias
for ss in seqs:
print lista_a_seq(ss)
In [ ]:
# Visualizo la primer secuencia como lista y como secuencia
s=seqs[0]
print s
print ':', lista_a_seq(s)
#print s, len(s)
In [19]:
# Visualizo un haplotipo como lista y como secuencia, idem con la combinación de haplotipos
print 'haplos:', haplos[4]
h=ref_haplo(haplos[4][1])
print 'h:',lista_a_seq(h)
#print h, len(h)
In [ ]:
print refseq
print h
In [20]:
# Diccionario de conversión a IUPAC
IUPAC={
'A':'A',
'C':'C',
'G':'G',
'T':'T',
'U':'U',
'AG':'R',
'CT':'Y',
'CG':'S',
'AT':'W',
'GT':'K',
'AC':'M',
'CGT':'B',
'AGT':'D',
'ACT':'H',
'ACG':'V',
'ACGT':'N',
'.':'.',
'-':'-'
}
In [24]:
#Ejemplo
IUPAC['AGT']
Out[24]:
In [22]:
# Cálculo del diccionario inverso de IUPAC
IUPAC_1=dict([ [x[1],x[0]] for x in list(IUPAC.items())])
In [25]:
#Ejemplo inverso
IUPAC_1['D']
Out[25]:
In [26]:
# Función de conversion de iupac a secuencia como lista, y ejemplo
s='ACGTRYSWMBHVN'
s='AHRV'
def iupac_a_seq(s):
t=[]
for x in s:
t.append(IUPAC_1[x])
return t
iupac_a_seq(s)
Out[26]:
In [27]:
# Función de conversión de secuencia como lista a IUPAC, con ejemplo
s=['A', 'C', 'G', 'T', 'GA', 'CT', 'CG', 'AT', 'AC', 'CGT', 'ACT', 'ACG', 'ACGT']
def lista_a_iupac(s):
''' Recibe una lista del tipo del ejemplo, la pasa a iupac, devuelve string
'''
t=[]
for x in s:
t.append(IUPAC[''.join(sorted(list(x)))])
return ''.join(t)
lista_a_iupac(s)
Out[27]:
In [28]:
# Funciones de comparación de secuencias y haplotipoa
def seq_compare(s,t): #se supone que estan alineados y tienen la misma longitud
match=0
nomatch=0
compared=0
nocompared=0
total=len(s)
#print 'S T i M N C'
for i, x in enumerate(s):
y=t[i]
if x=='-' or y=='-':
nocompared+=1
elif x==y:
match +=1
compared +=1
else:
nomatch +=1
compared +=1
#print x, t[i], i, match, nomatch, compared
return (match, nomatch, compared, nocompared, total)
def haplo_compare(s,t): #se supone que estan alineados y tienen la misma longitud
comp=[]
match=0
nomatch=0
compared=0
nocompared=0
total=len(s)
#print 'S T i M N C'
for i, x in enumerate(s):
yi=set(t[i])
xi=set(x)
if xi==set('-') or yi==set('-'):
nocompared+=1
comp +='-'
elif xi.issubset(yi):
match +=1
compared +=1
comp +='|'
else:
nomatch +=1
compared +=1
comp +='X'
#print x, t[i], i, match, nomatch, compared
return [match, nomatch, compared, nocompared, total, ''.join(comp)]
seq_compare('ASDFGH','ASEF-H')
Out[28]:
In [29]:
# Ejemplo usando la secuencia de referencia
seq_compare(refseq, refseq)
Out[29]:
In [30]:
# prueba del algoritmo de alineamiento
# para adaptarse a las diferencias sin introducir gap
# -> local y no global, penalizando altamente el gap (habria que buscar un mínimo)
f=open('alineamiento.txt','w')
from Bio import pairwise2
from Bio.pairwise2 import format_alignment
s=lista_a_seq(s)
h1=lista_a_seq(h[1])
alignments = pairwise2.align.localms(refseq,h1[:200],2,-1,-5,-5)
log(('#alignements=',len(alignments)))
for a in alignments:
log(('--------------'))
log((format_alignment(*a)))
log(('a=',a))
log(('seq_compare',seq_compare(a[0], a[1])))
# print '--------------'
# print format_alignment(*a)
# print a
# print seq_compare(a[0], a[1])
def alinear(s,t):
# print 'align:'
# print 's:', s
# print 't:', t
#log(('align:'))
#log(('s:', s))
#log(('t:', t))
al = pairwise2.align.localms(s,t,2,-1,-5,-5)
# for a in al:
# log(('--------------'))
# log((format_alignment(*a)))
# log(('a=',a))
# log(('seq_compare',seq_compare(a[0], a[1])))
#log((a[0][0], a[0][1]))
#log(('pairwise'))
#log((a))
return al[0][0], al[0][1]
f.close()
In [31]:
s
Out[31]:
In [32]:
print h[1]
In [ ]:
q='GTCATYTGDGTDGCDGGYTCDCC'
q='TTGAGBARRCCRTACATRAGC'
q='GGGCTRACRAARGCYCKRTCCC'
q='CCHTACAAYGCYGATTTYGACG'
q='GCTYATGTAYGGYYTVCTCAA'
q='GCTYTDCCKGGYATGAAYGTYCG'
r=iupac_a_seq(q)
print r
In [33]:
# Función que genera todas las secuencias posibles que machean con una secuencia degenerada
#esto es lo que hice para Suani!!!
def gen(s,a=''):
if len(s)>1: #print len(s)
for x in s[0]:
gen(s[1:],a+x)
else: #es el ultimo, cerrar
for x in s[0]:
print a + x
s=r
gen(s)
In [34]:
# Genero el archivo con los haplotipos generados
hap=open('haplos.txt','w')
for h in haplos:
hap.write(str(h[0]) + ' ' + lista_a_seq(h[1]) + '\n')
hap.close()
bak=''' for i,s in enumerate(seq): f=open(desc[i]+'comparacion.txt','w') print i,':',s for j,h in enumerate(haplos): print j, h1=h[1] s1=seq_a_lista(s) h1=h[1] log(('------------------\n')) #print '-------------------' x, y= alinear(lista_a_iupac(s1),lista_a_iupac(h1)) # res=haplo_compare(x,y) log((h[0])) log((x)) log((str(''.join(res[-1])))) log((y)) log(('match, nomatch, compared, nocompared, total, comp')) log((str(res[:-1])+'\n')) f.close() '''
In [35]:
# Procedimiento principal
# Main Loop
# Alinea y compara cada secuencia con cada par de haplotipos, y ordena los resultados por ranking
for i,s in enumerate(seq):
r=open(desc[i]+'comparacion.txt','w')
rank=[]
print i,':',s
for j,h in enumerate(haplos):
print j,
h1=h[1]
#print 'h1:',h1
#print lista_a_iupac(h1)
s1=seq_a_lista(s)
#print 's1:',s1
#print lista_a_iupac(s1)
#log(('------------------\n')) #print '-------------------'
x, y= alinear(lista_a_iupac(s1),lista_a_iupac(h1)) #
res=haplo_compare(x,y)
res.extend([h[0],x,y])
rank.append(res)
debug='''
log((h[0])) # los id
log((x))
log((str(''.join(res[-1]))))
log((y))
log(('match, nomatch, compared, nocompared, total, comp'))
log((str(res[:-1])+'\n'))
print((h[0])) # los id
print((x))
print((str(''.join(res[-1]))))
print((y))
print(('match, nomatch, compared, nocompared, total, comp'))
print((str(res[:-1])+'\n'))
'''
#print 'sorting...'
rank.sort(reverse=True)
#r=open('rank.txt','w')
for x in rank:
r.write('Haplotipo: '+ ' '.join(x[6])+'\n')
r.write(', '.join([str(xx) for xx in x[:5]])+'\n')
r.write(str(x[7])+'\n')
r.write(str(x[5])+'\n')
r.write(str(x[8])+'\n\n')
r.close()
In [36]:
# Escribo el ranking en el archivo rank.txt
rank.sort(reverse=True)
r=open('rank.txt','w')
for x in rank:
r.write('Haplotipo: '+ ' '.join(x[6])+'\n')
r.write(', '.join([str(xx) for xx in x[:5]])+'\n')
r.write(str(x[7])+'\n')
r.write(str(x[5])+'\n')
r.write(str(x[8])+'\n\n')
r.close()
In [37]:
f.close()
De aqui en adelante solo son pruebas y ejemplos de conversión
In [40]:
s=seq_a_lista(seq[0])
h=haplos[22]
t=h[1]
print 'S:', s
print ' '.join(s)
print 'S_iupac:'
si= ''.join(lista_a_iupac(s))
print 'si:', si
print 'haplo:',h[0]
print t
ti=lista_a_iupac(t)
print ti
ti =''.join(ti)
print 'ti:',ti
In [42]:
alignments = pairwise2.align.localms(si,ti,2,-1,-5,-5)
In [43]:
for a in alignments:
print '--------------'
print format_alignment(*a)
print a
print seq_compare(a[0], a[1])
In [44]:
print a[0]
print a[1]
print seq_compare(a[0],a[1])
In [45]:
print iupac_a_seq(a[0])
print iupac_a_seq(a[1])
print haplo_compare(iupac_a_seq(a[0]),iupac_a_seq(a[1]))
In [46]:
print haplo_compare(iupac_a_seq(a[1]),iupac_a_seq(a[0]))
In [47]:
p= iupac_a_seq(a[0])
q= iupac_a_seq(a[1])
r= haplo_compare(iupac_a_seq(a[0]),iupac_a_seq(a[1]))
print ''.join(p)
print ''.join(q)
print ''.join(r[-1])
In [48]:
f=open('prisec.txt')
ps=f.readlines()
f.close()
print len(ps)
a=ps[0]
b=ps[1]
print a
print b
In [ ]:
ab=unir_haplotipos(a,b)
In [ ]:
abseq=lista_a_seq(ab)
print abseq
In [ ]:
seq_a_lista(abseq)
In [ ]:
seq[0]
In [ ]:
seqs[0]