Ho mappato le letture dell'RNA sul genoma di riferimento, utilizzando LAST in modalità split e ho convertito l'allineamento MAF in SAM con maf-convert.
Il mio problema è che le trascrizioni non vengono riportate in modo spliced, il che significa che un transcript_ID viene segnalato più volte nel file SAM di allineamento con un identico flag bit per bit in $ 2 . Da quello che ho capito, ciò è dovuto al fatto che solo gli esoni sono mappati (un esone per riga) e riportati come allineamenti locali, poiché il software non può gestire la combinazione di schemi di giunzione esone-introne (ancora), il cui comportamento è chiaro da il filo SIGARO.
Per un esempio più concreto e visivo, si consideri la mappatura della trascrizione FBtr0344900 al genoma di riferimento come fornita da LAST:
$ cat last_aln .sam | grep FBtr0344900FBtr0344900 0 4 42774 100 384 = 144H * 0 0 TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATttttattttgatATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCAT * NM: i: 0 AS: i: 2304FBtr0344900 0 4 43231 100 384H144 = * 0 0 CTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT * NM: i: 0 AS: i: 864 Ed ecco la mappatura della stessa trascrizione FBtr0344900 fornita da STAR - software che riporta l'allineamento proprio nel modo in cui ne ho bisogno:
cat star_aln.sam | grep FBtr0344900
FBtr0344900 0 4 42774 255 384M73N144M * 0 * NH: i: 1 HI: i: 1 NM: i: 0 MD: Z: 528 Dalla discussione con l'autore, sembra che non riesco a ottenere ciò di cui ho bisogno direttamente dall'attuale versione LAST, e non è un problema tecnico. Quindi, devo modificare l'output da solo. L'obiettivo sarebbe quello di ottenere almeno una linea SIGAR che rappresenti una trascrizione completa.
La mia domanda è: conosci qualche software per farlo? Ciò di cui avrei bisogno è un file contenente una riga per trascrizioni mappate in modo univoco, con tre campi: transcript_ID , start position e CIGAR string .
Da parte mia, ho proceduto in questo modo:
1) estrai i campi di cui ho bisogno per l'interessato dal file SAM:
$ cut -f1,4,6FBtr0344900 42774 384 = 144HFBtr0344900 43231 384H144 = 2) dividi la riga CIGAR per rimuovere gli elementi che non mi interessano - Semplifico il comando qui, supponendo solo di averlo solo corrispondenze perfette (a cui sono interessato) e hard-clipping (a cui non sono interessato):
$ cut -f3 | sed 's / H / _ / g' | sed 's / = / = / g' | sed 's / \ w * _ \ s * //' | sed 's / = // g'384144 3) incolla il CIGAR modificato nel file originale, con paste , ottenendo:
$ incolla (1) (2) | cut -f1,2,4FBtr0344900 42774 384FBtr0344900 43231 144
4) unisci righe che iniziano con lo stesso transcript_id :
$ awk -F '\ t' -v OFS = '\ t' '{x = $ 1; $ 1 = ""; a [x] = a [x] $ 0} END {for (x in a) print x, a [x]} '| FBtr0344900 42774384 43231 144 5) calcola il nuovo sigaro, calcolando la lunghezza dell'introne come la formula aritmetica intron_length = (next_exon_start_coordinate - exon_length - previous_exon_start_coordinate) , in questo caso semplice sopra: intron_length = 43231-384-42774
$ awk '{printf ("% s", $ 1)}; {for (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 73 144 6) idealmente, con un metodo che io capisci, modifico la stringa CIGAR (aggiungendo M, N alternati alla fine di ogni campo eccetto il primo), questo è l'aspetto del file finale:
FBtr0344900 42774 384M73N144M Il problema del mio approccio di base è che:
- Non sono sicuro di come tenere conto del SAM basato su 1: devo aggiungere
+ 1ad ogniexon_start_coordinate? Non sembra, poiché l'uscita STAR ha la stessa identica stringa di sigari che ho calcolato sull'uscita STAR. - GRANDE problema: questo funzionerà solo per le letture mappate nel filamento diretto: come renderlo fattibile con le letture mappate nel filamento inverso? Se mantengo il mio approccio attuale, avrò dimensioni di introni negative ...
Qualsiasi suggerimento è benvenuto!