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"

No hay comentarios:

Publicar un comentario