Vim puede ser una herramienta dura al principio, pero a poco que se maneja un poco parece evidente la potencia de este programa para el manejo de archivos de texto plano.
Ls forma de borrar múltiples líneas con Vim no puede ser más sencilla. Para acceder a los comando de Vim basta con pulsar ESC para salir del modo INSERTAR, si es que estamos en él, tecleamos la orden correspondiente. Veamos unos cuantos ejemplos:
Borrar desde la línea donde está el cursor (.) hasta el final del archivo ($)
:.,$d
Borrar desde la línea 1 a la 39 (incluidas ambas)
:1,39d
Borrar desde la línea donde está el cursor (.) hasta el final, menos las 10 últimas líneas
:.,$-10d
Borrar desde la línea 1 hasta la línea donde está el cursor (.) incluida
:1,.d
Fácil, ¿verdad?
NGS Bioinformática
Bioinformatica en español. Aquí almaceno los trucos o piezas de código que me han servido, cosas interesantes que he visto o leído o inquietudes personales en este mundillo.
lunes, 30 de marzo de 2015
viernes, 12 de abril de 2013
Hashes y diccionarios anidados 'al vuelo'
Las estructuras de datos de Perl conocidas como 'hashes' tienen su equivalente en Python con el nombre de 'diccionarios'. En ellos se accede a un valor (p.e. un número de teléfono) mediante una clave (p.e. un nombre). Los hashes/diccionarios planos son fáciles de implementar en ambos lenguajes:
PERL
my %hash = (); # Definimos un hash vacío
$hash{'Fulano'} = '123456789'; # Introducimos el valor '123456789' con la clave de acceso 'Fulano'
PYTHON
diccionario = {} # Definimos un diccionario vacío
diccionario['Fulano'] = '123456789' # Introducimos el valor '123456789' con la clave de acceso 'Fulano'
Perl va perdiendo gradualmente posiciones frente a Python. Sin embargo, cuando se trata de rellenar valores en diccionarios anidados 'al vuelo', Perl tiene una flexibilidad que no tiene Python. En Perl sólo tenemos que definir un 'hash' vacío y luego podemos ir rellenándolo al nivel que nosotros queramos, haciendo tantos niveles anidados como sean necesarios 'on the fly'. Sin embargo Python necesita que le definamos a priori los niveles de los que va a constar el diccionario. La forma inmediata de hacerlo cada vez que metemos un entrada en el diccionario sería la siguiente:
>>> diccionario = {'clave1':{'clave2':{'clave3':'valor'}}}
>>> diccionario['clave1']['clave2']['clave3']
'valor'
Sin embargo si tenemos que iterar en una gran cantidad de datos con distintos niveles de anidamiento, hay una forma más fácil mediante el módulo defaultdict y las funciones lambda que merecerán un próximo post:
>>> from collections import defaultdict
>>> diccionario = defaultdict(lambda: defaultdict())
De esta forma hemos definido un diccionario con 2 niveles interpuestos entre la raíz y el valor final. Por ejemplo:
>>> diccionario['a']['b'] = 'ab'
>>> dicccionario['a']['d']
'ad'
Sin embargo, si intentamos introducir por ejemplo el siguiente registro en el diccionario, se producirá un mensaje de error:
>>> diccionario['a']['d']['m'] = 'acm'
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: 'str' object does not support item assignment
Si quisiésemos un diccionario con 3 niveles interpuestos entre la raiz y el valor final deberíamos hacer lo siguiente:
>>> diccionario = defaultdict(lambda: defaultdict(lambda: defaultdict()))
Veamos un ejemplo a partir de un archivo SAM/BAM de un alineamiento de lecturas single-end. Supongamos por queremos contar las lecturas del alineamiento que caen en cada cromosoma y, dentro de cada cromosoma, en que hebra. Hay por tanto un total de 2 niveles interpuestos entre la raiz y el valor final:
diccionario/hash -> cromosoma -> hebra -> valor
La forma de informarnos de en que hebra estamos ya la hemos visto en un post anterior.
PERL
my %hash = {};
while (my $line = <SAM>) { # Lee el archivo SAM línea a línea
chomp $line;
if ($line =~ /^@/) {next;} # Descarta las líneas de cabecera
else {
my @chunks = split ("\t",$line);
my $bitwise = $chunks[1]; # Selecciona la columna del 'bitwise flag'
my $chr = $chunks[2]; # Selecciona el cromosoma
my $byte = sprintf("%012b", $bitwise); # Convierte el número en una cadena de 12 bits
my $strand = substr($byte, -5, 1); # Selecciona el bit que indica la hebra de la lectura
my $strand_name = '';
if ($strand) {$strand_name = 'reverse';}
else {$strand_name = 'forward';}
$hash{$chr}{$strand_name}++; # Sumamos una lectura a la entrada del hash: facilísimo !!!!!!
}
}
PYTHON
diccionario = defaultdict(lambda: defaultdict()) # Definimos un diccionario vacío con 2 niveles de anidamiento
for line in sam_file:
if line[0] == '@': # Descarta las líneas de cabecera
continue
else:
chunks = line.split()
bitwise = int(chunks[1]) # Selecciona la columna del 'bitwise flag'
chr = str(chunks[2]) # Selecciona el cromosoma
bits = bin(bitwise) # Transforma a binario, pero con un prefijo '0b'
bits = bits[2:] # Quita el prefijo '0b'
byte = "{0:0>12}".format(bits) # Añade ceros hasta obtener una cadena de 12 bits
strand = int(byte[-5]) # Selecciona el bit que indica la hebra de la lectura
if strand:
strand_name = 'reverse'
else:
strand_name = 'forward'
diccionario[chr][strand_name]++ # Sumamos una lectura a la entrada del diccionario
Aunque hemos podido hacer la implementación de diccionarios anidados en Python algo más llevadero, seguimos constreñidos por la necesidad de definir previamente el número de niveles de los que constará y, por tanto, escribir un número equivalente de 'defaultdict' en nuestra definición del diccionario. Si definimos varios diccionarios anidados, este trabajo se multiplica. Sin embargo, podemos definir una función recursiva que permite definir un diccionario con el número máximo de niveles de anidamiento sin necesidad de escribir tantos 'defaultdict' y que valdrá para cualquier diccionario con el número de niveles que queramos:
>>> from collections import defaultdict
>>> def diccionarioAnidado():
... return defaultdict(diccionarioAnidado)
Hagamos una prueba:
>>> diccionario2 = diccionarioAnidado()
>>> diccionario2['a']['b']['c']['d'] = 'abcd'
>>> diccionario2['a']['b']['c']['d']
'abcd'
>>> diccionario3 = diccionarioAnidado()
>>> diccionario3['x']['y']['x'] = 'xyz'
>>> diccionario3['x']['y']['x']
'xyz'
¡Esto es otra cosa!
PERL
my %hash = (); # Definimos un hash vacío
$hash{'Fulano'} = '123456789'; # Introducimos el valor '123456789' con la clave de acceso 'Fulano'
PYTHON
diccionario = {} # Definimos un diccionario vacío
diccionario['Fulano'] = '123456789' # Introducimos el valor '123456789' con la clave de acceso 'Fulano'
Perl va perdiendo gradualmente posiciones frente a Python. Sin embargo, cuando se trata de rellenar valores en diccionarios anidados 'al vuelo', Perl tiene una flexibilidad que no tiene Python. En Perl sólo tenemos que definir un 'hash' vacío y luego podemos ir rellenándolo al nivel que nosotros queramos, haciendo tantos niveles anidados como sean necesarios 'on the fly'. Sin embargo Python necesita que le definamos a priori los niveles de los que va a constar el diccionario. La forma inmediata de hacerlo cada vez que metemos un entrada en el diccionario sería la siguiente:
>>> diccionario = {'clave1':{'clave2':{'clave3':'valor'}}}
>>> diccionario['clave1']['clave2']['clave3']
'valor'
Sin embargo si tenemos que iterar en una gran cantidad de datos con distintos niveles de anidamiento, hay una forma más fácil mediante el módulo defaultdict y las funciones lambda que merecerán un próximo post:
>>> from collections import defaultdict
>>> diccionario = defaultdict(lambda: defaultdict())
De esta forma hemos definido un diccionario con 2 niveles interpuestos entre la raíz y el valor final. Por ejemplo:
>>> diccionario['a']['b'] = 'ab'
>>> dicccionario['a']['d']
'ad'
Sin embargo, si intentamos introducir por ejemplo el siguiente registro en el diccionario, se producirá un mensaje de error:
>>> diccionario['a']['d']['m'] = 'acm'
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: 'str' object does not support item assignment
Si quisiésemos un diccionario con 3 niveles interpuestos entre la raiz y el valor final deberíamos hacer lo siguiente:
>>> diccionario = defaultdict(lambda: defaultdict(lambda: defaultdict()))
Veamos un ejemplo a partir de un archivo SAM/BAM de un alineamiento de lecturas single-end. Supongamos por queremos contar las lecturas del alineamiento que caen en cada cromosoma y, dentro de cada cromosoma, en que hebra. Hay por tanto un total de 2 niveles interpuestos entre la raiz y el valor final:
diccionario/hash -> cromosoma -> hebra -> valor
La forma de informarnos de en que hebra estamos ya la hemos visto en un post anterior.
PERL
my %hash = {};
while (my $line = <SAM>) { # Lee el archivo SAM línea a línea
chomp $line;
if ($line =~ /^@/) {next;} # Descarta las líneas de cabecera
else {
my @chunks = split ("\t",$line);
my $bitwise = $chunks[1]; # Selecciona la columna del 'bitwise flag'
my $chr = $chunks[2]; # Selecciona el cromosoma
my $byte = sprintf("%012b", $bitwise); # Convierte el número en una cadena de 12 bits
my $strand = substr($byte, -5, 1); # Selecciona el bit que indica la hebra de la lectura
my $strand_name = '';
if ($strand) {$strand_name = 'reverse';}
else {$strand_name = 'forward';}
$hash{$chr}{$strand_name}++; # Sumamos una lectura a la entrada del hash: facilísimo !!!!!!
}
}
PYTHON
diccionario = defaultdict(lambda: defaultdict()) # Definimos un diccionario vacío con 2 niveles de anidamiento
for line in sam_file:
if line[0] == '@': # Descarta las líneas de cabecera
continue
else:
chunks = line.split()
bitwise = int(chunks[1]) # Selecciona la columna del 'bitwise flag'
chr = str(chunks[2]) # Selecciona el cromosoma
bits = bin(bitwise) # Transforma a binario, pero con un prefijo '0b'
bits = bits[2:] # Quita el prefijo '0b'
byte = "{0:0>12}".format(bits) # Añade ceros hasta obtener una cadena de 12 bits
strand = int(byte[-5]) # Selecciona el bit que indica la hebra de la lectura
if strand:
strand_name = 'reverse'
else:
strand_name = 'forward'
diccionario[chr][strand_name]++ # Sumamos una lectura a la entrada del diccionario
Aunque hemos podido hacer la implementación de diccionarios anidados en Python algo más llevadero, seguimos constreñidos por la necesidad de definir previamente el número de niveles de los que constará y, por tanto, escribir un número equivalente de 'defaultdict' en nuestra definición del diccionario. Si definimos varios diccionarios anidados, este trabajo se multiplica. Sin embargo, podemos definir una función recursiva que permite definir un diccionario con el número máximo de niveles de anidamiento sin necesidad de escribir tantos 'defaultdict' y que valdrá para cualquier diccionario con el número de niveles que queramos:
>>> from collections import defaultdict
>>> def diccionarioAnidado():
... return defaultdict(diccionarioAnidado)
Hagamos una prueba:
>>> diccionario2 = diccionarioAnidado()
>>> diccionario2['a']['b']['c']['d'] = 'abcd'
>>> diccionario2['a']['b']['c']['d']
'abcd'
>>> diccionario3 = diccionarioAnidado()
>>> diccionario3['x']['y']['x'] = 'xyz'
>>> diccionario3['x']['y']['x']
'xyz'
¡Esto es otra cosa!
domingo, 31 de marzo de 2013
Bitwise flag del formato SAM
La segunda columna de un archivo en formato SAM está reservado al denominado bitwise flag. Este es un número en formato decimal, que debe ser traducido a formato binario de 12 bits. Una vez traducido, cada uno de los bits, de derecha a izquierda, confirma (bit = 1) o rechaza (bit = 0) cada una de las declaraciones siguientes que se hacen sobre una lectura:
bit 1: La pareja ha tenido múltiples lecturas en la secuenciación
bit 2: Cada lectura del par ha sido correctamente alineada según el alineador
bit 3: Esta lectura no ha sido mapeada
bit 4: La pareja de esta lectura no ha sido mapeada
bit 5: La lectura mapea en la hebra reversa
bit 6: La pareja de la lectura mapea en la hebra reversa
bit 7: Esta lectura es la primera del par
bit 8: Esta lectura es la última del par
bit 9: Esta lectura tiene un alineamiento secundario
bit 10: Esta lectura no ha pasado los controles de calidad
bit 11: Esta lectura es un duplicado de PCR u óptico
bit 12: Esta lectura tiene un alineamiento suplementario
Este 'flag' puede ser muy útil. Fundamentalmente me ha servido hasta el momento para diferenciar la hebra de mapeo (bit #5) o para filtrar los duplicados de PCR (bit #11).
¿Como podemos acceder cómodamente a la información de este flag? Incluyo dos ejemplos de como seleccionar una lectura según su hebra de mapeo en Perl y Python que espero que sean útiles:
PERL
while (my $line = <SAM>) { # Lee el archivo SAM línea a línea
chomp $line;
if ($line =~ /^@/) {next;} # Descarta las líneas de cabecera
else {
my @chunks = split ("\t",$line);
my $bitwise = $chunks[1]; # Selecciona la columna del 'bitwise flag'
my $byte = sprintf("%012b", $bitwise); # Convierte el número en una cadena de 12 bits
my $strand = substr($byte, -5, 1); # Selecciona el bit que indica la hebra de la lectura
if ($strand) {print "La lectura mapea en la hebra reversa\n";}
else {print "La lectura mapea en la hebra directa\n";}
}
}
PYTHON
for line in sam_file:
if line[0] == '@': # Descarta las líneas de cabecera
continue
else:
chunks = line.split()
bitwise = int(chunks[1]) # Selecciona la columna del 'bitwise flag'
bits = bin(bitwise) # Transforma a binario, pero con un prefijo '0b'
bits = bits[2:] # Quita el prefijo '0b'
byte = "{0:0>12}".format(bits) # Añade ceros hasta obtener una cadena de 12 bits
strand = int(byte[-5]) # Selecciona el bit que indica la hebra de la lectura
if strand:
print "La lectura mapea en la hebra reversa"
else:
print "La lectura mapea en la hebra directa"
bit 1: La pareja ha tenido múltiples lecturas en la secuenciación
bit 2: Cada lectura del par ha sido correctamente alineada según el alineador
bit 3: Esta lectura no ha sido mapeada
bit 4: La pareja de esta lectura no ha sido mapeada
bit 5: La lectura mapea en la hebra reversa
bit 6: La pareja de la lectura mapea en la hebra reversa
bit 7: Esta lectura es la primera del par
bit 8: Esta lectura es la última del par
bit 9: Esta lectura tiene un alineamiento secundario
bit 10: Esta lectura no ha pasado los controles de calidad
bit 11: Esta lectura es un duplicado de PCR u óptico
bit 12: Esta lectura tiene un alineamiento suplementario
Este 'flag' puede ser muy útil. Fundamentalmente me ha servido hasta el momento para diferenciar la hebra de mapeo (bit #5) o para filtrar los duplicados de PCR (bit #11).
¿Como podemos acceder cómodamente a la información de este flag? Incluyo dos ejemplos de como seleccionar una lectura según su hebra de mapeo en Perl y Python que espero que sean útiles:
PERL
while (my $line = <SAM>) { # Lee el archivo SAM línea a línea
chomp $line;
if ($line =~ /^@/) {next;} # Descarta las líneas de cabecera
else {
my @chunks = split ("\t",$line);
my $bitwise = $chunks[1]; # Selecciona la columna del 'bitwise flag'
my $byte = sprintf("%012b", $bitwise); # Convierte el número en una cadena de 12 bits
my $strand = substr($byte, -5, 1); # Selecciona el bit que indica la hebra de la lectura
if ($strand) {print "La lectura mapea en la hebra reversa\n";}
else {print "La lectura mapea en la hebra directa\n";}
}
}
PYTHON
for line in sam_file:
if line[0] == '@': # Descarta las líneas de cabecera
continue
else:
chunks = line.split()
bitwise = int(chunks[1]) # Selecciona la columna del 'bitwise flag'
bits = bin(bitwise) # Transforma a binario, pero con un prefijo '0b'
bits = bits[2:] # Quita el prefijo '0b'
byte = "{0:0>12}".format(bits) # Añade ceros hasta obtener una cadena de 12 bits
strand = int(byte[-5]) # Selecciona el bit que indica la hebra de la lectura
if strand:
print "La lectura mapea en la hebra reversa"
else:
print "La lectura mapea en la hebra directa"
Suscribirse a:
Comentarios (Atom)