J'ai mappé les lectures d'ARN sur le génome de référence, en utilisant LAST en mode fractionné, et j'ai converti l'alignement MAF en SAM avec maf-convert.
Mon problème est que les transcriptions ne sont pas rapportées de manière épissée, ce qui signifie qu'un transcript_ID est signalé plusieurs fois dans le fichier SAM d'alignement avec un drapeau bit à bit identique dans 2 $ . D'après ce que je comprends, cela est dû au fait que seuls les exons sont mappés (un exon par ligne) et signalés comme des alignements locaux, car le logiciel ne peut pas gérer (encore) la combinaison de modèles d'épissage exon-intron, dont le comportement est clair la ficelle CIGAR.
Pour un exemple plus concret et visuel, considérons le mappage du transcript FBtr0344900 au génome de référence tel que donné par 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 Et voici le mappage du même transcript FBtr0344900 que celui donné par STAR - logiciel qui rapporte l'alignement comme j'en ai besoin:
chat star_aln.sam | grep FBtr0344900
FBtr0344900 0 4 42774 255 384M73N144M * 0 * NH: i: 1 HI: i: 1 NM: i: 0 MD: Z: 528 De la discussion avec l'auteur, il semble que je ne puisse pas obtenir ce dont j'ai besoin directement à partir de la dernière version actuelle, et ce n'est pas un problème technique. Je dois donc modifier moi-même la sortie. Le but serait d'obtenir au moins une ligne CIGAR représentant une transcription complète.
Ma question est la suivante: connaissez-vous un logiciel pour faire cela? Ce dont j'aurais besoin est un fichier contenant une ligne par transcription mappée de manière unique, comportant trois champs: transcript_ID , position de départ et chaîne CIGAR . p>
De mon côté, j'ai procédé comme ceci:
1) extraire les champs dont j'ai besoin pour l'intéressé du fichier SAM:
$ cut -f1,4,6FBtr0344900 42774384 = 144HFBtr0344900 43231 384H144 = 2) diviser la ligne CIGAR pour supprimer les éléments qui ne m'intéressent pas - Je simplifie la commande ici, en supposant simplement que je n'ai correspondances parfaites (qui m'intéressent) et hard-clipping (qui ne m'intéressent pas):
$ cut -f3 | sed 's / H / _ / g' | sed 's / = / = / g' | sed 's / \ w * _ \ s * //' | sed 's / = // g'384144 3) collez le CIGAR modifié dans le fichier d'origine, avec collez , ce qui donne:
$ coller (1) (2) | coupe -f1,2,4FBtr0344900 42774 384FBtr0344900 43231144
4) fusionner les lignes commençant par le même 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 42774 384 43231 144 5) calculer le nouveau cigare, en calculant la longueur d'intron comme la formule arithmétique intron_length = (next_exon_start_coordinate - exon_length - previous_exon_start_coordinate) , dans ce cas simple ci-dessus: intron_length = 43231-384-42774
$ awk '{printf ("% s", $ 1)}; {pour (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) idéalement, avec une méthode que je vais comprendre, je vais modifier la chaîne CIGAR (en ajoutant un autre M, N à la fin de chaque champ sauf le premier), voici à quoi le fichier final devrait ressembler:
FBtr0344900 42774 384M73N144M Le problème de mon approche de base est que:
- Je ne sais pas comment rendre compte du SAM basé sur 1: dois-je ajouter
+ 1à chaqueexon_start_coordinate? Cela ne ressemble pas, car la sortie STAR a exactement la même chaîne de cigare que j'ai calculée sur la sortie STAR. - GRAND problème: Cela ne fonctionnera que pour les lectures mappées dans le brin avant: comment le rendre réalisable avec des lectures mappées dans le brin inverse? Si je garde mon approche actuelle, j'aurai des tailles d'introns négatives ...
Toute suggestion est la bienvenue!