Mapeé lecturas de ARN para referenciar el genoma, usando LAST en modo dividido, y convertí la alineación MAF a SAM con maf-convert.

Mi problema es que las transcripciones no se informan de manera empalmada, lo que significa que un transcript_ID se informa varias veces en el archivo SAM de alineación con un indicador bit a bit idéntico en $ 2 . Por lo que entiendo, esto se debe al hecho de que solo se mapean los exones (un exón por fila) y se informan como alineaciones locales, ya que el software no puede lidiar con la combinación de patrones de empalme exón-intrón (todavía), cuyo comportamiento está claro a partir de la cadena CIGAR.

Para un ejemplo más concreto y visual, consideremos el mapeo de la transcripción FBtr0344900 al genoma de referencia dado por LAST:

  $ cat last_aln .sam | grep FBtr0344900FBtr0344900 0 4 42 774 100 384 = 144H * 0 0 TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATttttattttgatATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCAT * NM: i: 0 AS: i: 2304FBtr0344900 0 4 43 231 100 384H144 = * 0 0 CTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT * NM: i: 0 AS: i: 864  

Y aquí está el mapeo de la misma transcripción FBtr0344900 dada por STAR - software que informa la alineación de la manera correcta:

  cat star_aln.sam | grep FBtr0344900
FBtr0344900 0 4 42774255 384M73N144M * 0 * NH: i: 1 HI: i: 1 NM: i: 0 MD: Z: 528  

De discusión con el autor, parece que no puedo obtener lo que necesito directamente de la ÚLTIMA versión actual, y no es un problema técnico. Entonces, tengo que modificar la salida yo mismo. El objetivo sería obtener al menos una línea CIGAR que represente una transcripción completa.

Mi pregunta es, ¿conoce algún software para hacer esto? Lo que necesitaría es un archivo que contenga una línea por cada transcripción mapeada de forma única, con tres campos: transcript_ID , posición de inicio y CIGAR string .

Por mi parte, procedí así:

1) extraer los campos que necesito para los interesados ​​del archivo SAM:

  $ cut -f1,4,6FBtr0344900 42774 384 = 144HFBtr0344900 43231 384H144 =  

2) dividir la línea CIGAR para eliminar elementos que no me interesan: simplifico el comando aquí, suponiendo que solo tengo coincidencias perfectas (que me interesan) y recortes duros (que no me interesan):

  $ cut -f3 | sed 's / H / _ / g' | sed 's / = / = / g' | sed 's / \ w * _ \ s * //' | sed 's / = // g'384144  

3) pegue el CIGAR modificado en el archivo original, con paste , resultando en:

  $ paste (1) (2) | cortar -f1,2,4FBtr0344900 42774 384FBtr0344900 43231 144  

4) fusionar filas que comienzan con el mismo transcript_id :

  $ awk -F '\ t' -v OFS = '\ t' '{x = $ 1; $ 1 = ""; a [x] = a [x] $ 0} FIN {para (x en a) imprimir x, a [x]} '| FBtr0344900 42774384 43231144 

5) Calcule el nuevo puro, calculando la longitud del intrón como la fórmula aritmética intron_length = (next_exon_start_coordinate - exon_length - previous_exon_start_coordinate) , en este caso simple arriba: intron_length = 43231-384-42774

  $ awk '{printf ("% s", $ 1)}; {para (i = 4; i< = NF; i + = 2) {printf ("% s% d% s% d", OFS, $ (i-1), OFS, $ i - $ (i-1) - $ (i-2))}}; {printf ("% s% d% s% s", OFS, $ NF, OFS, RS)} 'FBtr0344900 384 73144  

6) idealmente, con algún método que averiguar, modificaré la cadena CIGAR (agregando M, N alternativos al final de cada campo excepto el primero), así es como debería verse el archivo final:

  FBtr0344900 42774 384M73N144M  

El problema de mi enfoque básico es que:

1. No estoy seguro de cómo contabilizar el SAM basado en 1: ¿debo agregar + 1 en cada exon_start_coordinate ? No lo parece, ya que la salida de STAR tiene exactamente la misma cuerda de cigarro que calculé en la salida de STAR.
2. GRAN problema: esto solo funcionará para las lecturas mapeadas en la hebra delantera: ¿cómo hacerlo factible con lecturas mapeadas en la hebra inversa? Si mantengo mi enfoque actual, tendré tamaños de intrones negativos ...

¡Cualquier sugerencia es bienvenida!