Découverte de nouveaux transcrits

Cours

Exercice 6 : Recherche de nouveaux transcrits

1. Manipulation du GTF, se familiariser avec sa référence : A partir du fichier ITAG_pre2.3_gene_models_Ch6.gtf, compter combien il y a de transcrits.

   Penser a utiliser la commande cut sur colonne 9 pour récuperer les attributs, puis piper | ce résultat et faire un cut selon « ; » pour ne récurérer que le nom du transcript, enfin trier et compter le compter le nombre de lignes)

   Solution

   cd star-index/
   cut -f 9 ITAG_pre2.3_gene_models_Ch6.gtf | cut -d ';' -f 2 | sort -u | wc
   

   --> 2813 transcrits

   cd ..
   

   Explication de cette commande , le fichier contient les lignes suivantes:

    SL2.40ch06    ITAG_eugene    exon    3688    4407    .    +    .    gene_id "Solyc06g005000.2.1"; transcript_id "Solyc06g005000.2.1";
    SL2.40ch06    ITAG_eugene    exon    4853    5604    .    +    .    gene_id "Solyc06g005000.2.1"; transcript_id "Solyc06g005000.2.1";
    SL2.40ch06    ITAG_eugene    exon    7663    7998    .    +    .    gene_id "Solyc06g005000.2.1"; transcript_id "Solyc06g005000.2.1";
    SL2.40ch06    ITAG_eugene    exon    9148    9408    .    +    .    gene_id "Solyc06g005010.1.1"; transcript_id "Solyc06g005010.1.1";
    SL2.40ch06    ITAG_eugene    exon    13310    13330    .    +    .    gene_id "Solyc06g005020.1.1"; transcript_id "Solyc06g005020.1.1";
   

   • Récupération de la 9ieme colonne (qui contient le transcript_id)

     $ cut -f 9 star-index/ITAG_pre2.3_gene_models_Ch6.gtf | head
     gene_id "Solyc06g005000.2.1"; transcript_id "Solyc06g005000.2.1";
     gene_id "Solyc06g005000.2.1"; transcript_id "Solyc06g005000.2.1";
     gene_id "Solyc06g005000.2.1"; transcript_id "Solyc06g005000.2.1";
     gene_id "Solyc06g005010.1.1"; transcript_id "Solyc06g005010.1.1";
     gene_id "Solyc06g005020.1.1"; transcript_id "Solyc06g005020.1.1";
     gene_id "Solyc06g005020.1.1"; transcript_id "Solyc06g005020.1.1";
     

   • Récupération de la 2ieme colonne à partir de la sortie précédente en découpant sur le ';'

     $ cut -f 9 star-index/ITAG_pre2.3_gene_models_Ch6.gtf  | cut -d ';' -f 2  | head
     transcript_id "Solyc06g005000.2.1"
     transcript_id "Solyc06g005000.2.1"
     transcript_id "Solyc06g005000.2.1"
     transcript_id "Solyc06g005010.1.1"
     transcript_id "Solyc06g005020.1.1"
     transcript_id "Solyc06g005020.1.1"
     

   • Trie de facon unique (déduplique les lignes dupliqué)

     $ cut -f 9 star-index/ITAG_pre2.3_gene_models_Ch6.gtf  | cut -d ';' -f 2  | sort -u | head
     transcript_id "Solyc06g005000.2.1"
     transcript_id "Solyc06g005000.2.1"
     transcript_id "Solyc06g005010.1.1"
     transcript_id "Solyc06g005020.1.1"
     transcript_id "Solyc06g005030.1.1"
     transcript_id "Solyc06g005040.1.1"
     transcript_id "Solyc06g005050.2.1"
     

     • Compte le nombre de lignes

       cut -f 9 star-index/ITAG_pre2.3_gene_models_Ch6.gtf  | cut -d ';' -f 2  | sort -u | wc -l
       

2. Se positionner sur un noeud (si vous n'y êtes pas déjà)

   Solution

   srun --pty -c 4 -t 04:00:00 bash
   

3. Charger le module stringtie

   Solution

   $ search_module stringtie
   bioinfo/stringtie-1.3.2
   bioinfo/stringtie-2.0.5
   bioinfo/stringtie-2.1.1
   $ module load bioinfo/stringtie-2.1.1
   

4. Lancer stringtie sur chaque échantillon

   Solution

   stringtie WTAligned.sortedByCoord.out.bam -G star-index/ITAG_pre2.3_gene_models_Ch6.gtf -o transcriptWT.gtf 
   stringtie MTAligned.sortedByCoord.out.bam -G star-index/ITAG_pre2.3_gene_models_Ch6.gtf -o transcriptMT.gtf
   

5. Fusionner les annotations obtenus dans un seul fichier. Syntaxe : stringtie --merge ...

   Solution

   stringtie --merge -G star-index/ITAG_pre2.3_gene_models_Ch6.gtf -o merged.gtf transcriptWT.gtf transcriptMT.gtf
   

6. Combien de transcrits obtenez vous ? Comparer ce résultat au comptage obtenu à la première question.

   Solution

   $ awk '{$3=="exon"; print $0}' transcriptMT.gtf | grep -v '^#' | cut -f 9 | cut -f 2 -d ';' | sort -u | wc -l
   3442
   $ awk '{$3=="exon"; print $0}' transcriptWT.gtf | grep -v '^#' | cut -f 9 | cut -f 2 -d ';' | sort -u | wc -l
   3128
   $ awk '{$3=="exon"; print $0}' merged.gtf | grep -v '^#' | cut -f 9 | cut -f 2 -d ';' | sort -u | wc -l
   4668
   

7. Télécharger les nouveaux gtf sur votre ordinateur et charger le dans IGV, aller voir les zones précédement trouvées.

8. L'outil gffcompare permet d'obtenir une comparaison entre deux fichiers d'annotation.

   • Charger le module gffcompare

   • Comparer les fichiers MT et WT pour identifier des isoformes différentes dans les deux conditions

     Solution

     module load bioinfo/gffcompare-v0.11.4
     gffcompare -r transcriptMT.gtf transcriptWT.gtf
     more gffcmp.transcriptWT.gtf.tmap
     

   • Afficher le fichier tmap et recherche quelque exemple a aller visualiser dans igv

     • SL2.40ch06:1,833,134-1,847,126: (m) retention d'intron
     • SL2.40ch06:34,211,417-34,256,456 : (i) contenu dans un intron de la ref.

       Utiliser par exemple la commande suivante pour extraire les lignes avec 'i' dans la 3ieme colonne: awk -F'\t' '$3=="i"' gffcmp.transcriptWT.gtf.tmap

   • Aller dans IGV visualiser ces différences

   • Comparer le fichier mergé final avec la reférence initiale.

     Solution

     gffcompare -r ITAG_pre2.3_gene_models_Ch6.gtf merged.gtf 
     more gffcmp.merged.gtf.tmap
     awk -F'\t' '$3=="i"'  gffcmp.merged.gtf.tmap 
     awk -F'\t' '$3=="j"'  gffcmp.merged.gtf.tmap 
     awk -F'\t' '$3=="m"'  gffcmp.merged.gtf.tmap
     

   • Aller dans IGV visualiser ces différences