Objectif

On se propose de d'implémenter une commande nb_exons capable de renvoyer pour chaque transcript le nombre d'exons rencontrés.

Préparation de l'éditeur

Lancez le programme geany. Créez un nouveau fichier. A la première ligne de ce fichier, indiquez:

# -*- coding: utf8 -*-

Sauvegardez le fichier sous le nom "nb_exons.py".

Téléchargement

A titre d'exemple, nous utiliserons un fichier gtf. Ce fichier qui contient les coordonnées des exons de 5000 transcrits dans le génome humain peut être téléchargé ici: hg38_5k.gtf

Format du fichier

Le format GTF est un fichier tabulé (i.e. dont les colonnes sont séparées par des tabulations) qui contient des informations concernant des éléments génomiques (souvent des transcripts et leurs exons). Les informations portées par les colonnes sont les suivantes:

  • Colonne 1 - Nom du chromosome.
  • Colonne 2 - Source de l'annotation.
  • Colonne 3 - Type d'élément (exon, CDS, 5'UTR, 3'UTR,...).
  • Colonne 4 - Position de départ (Start).
  • Colonne 5 - Position de fin (End).
  • Colonne 6 - Score.
  • Colonne 7 - Le brin (+ ou -).
  • Colonne 8 - Frame (0,1,2). '0' indique que la première base de l'élément est la première base d'un codon, '1' que la deuxième base est la première base d'un codon, etc...
  • Colonne 9 - Attribute. Une liste de couples clefs/valeurs séparées par un point virgule.
1       ensembl exon    182393  182746  .       +       .       gene_id "ENSG00000279928"; transcript_id "ENST00000624431";
1       ensembl exon    183114  183240  .       +       .       gene_id "ENSG00000279928"; transcript_id "ENST00000624431";
1       ensembl exon    183922  184158  .       +       .       gene_id "ENSG00000279928"; transcript_id "ENST00000624431";
1       ensembl exon    195263  195411  .       -       .       gene_id "ENSG00000279457"; transcript_id "ENST00000623083";
1       ensembl exon    925741  925800  .       +       .       gene_id "ENSG00000187634"; transcript_id "ENST00000616016";
1       ensembl exon    925741  925800  .       +       .       gene_id "ENSG00000187634"; transcript_id "ENST00000618181";
1       ensembl exon    925741  925800  .       +       .       gene_id "ENSG00000187634"; transcript_id "ENST00000618323";
1       ensembl exon    925741  925800  .       +       .       gene_id "ENSG00000187634"; transcript_id "ENST00000620200";
1       ensembl exon    925922  926013  .       +       .       gene_id "ENSG00000187634"; transcript_id "ENST00000616016";
1       ensembl exon    925922  926013  .       +       .       gene_id "ENSG00000187634"; transcript_id "ENST00000618181";
1       ensembl exon    925922  926013  .       +       .       gene_id "ENSG00000187634"; transcript_id "ENST00000618323";

Solution

Le pseudo -code ci-dessous illustre la stratégie du programme:

  • Créer un dictionnaire permettant de mémoriser le nom du transcrit (clef) et le nombre d'exons associés.
  • Ouvrir le fichier.
  • Parcourir le fichier ligne à ligne.
  • Substituer les caractères précédant 'transcript_id "' par une chaîne de caractères vides.
  • Substituer les caractères succédant l'identifiant du transcript par une chaîne de caractères vides.
  • Si le transcript_id courant n'est pas dans l'un des dictionnaire :
    • Initialiser la clef du dictionnaires avec le nom du transcrit et la valeur 1.
  • Sinon:
    • Ajouter 1 à la valeur correspondant à la clef du transcrit dans le dictionnaire.
  • Imprimer transcrit et la valeur correspondant à sa clef dans le dictionnaire
In []:
# -*- coding: utf-8 -*-
# La ligne précédente permet d'utiliser des caractères accentués
# dans les commentaires

"""
Author D. Puthier
date: Wed Oct 28 16:44:29 CET 2015
program: nb_exons.py
This program is intended to compute the number of exons in the gtf file.
"""
# On importe le module re
# afin de pouvoir rechercher
# une expression régulière dans une
# chaîne de caractères.
import re

# Ce dictionnaire contiendra le nombre d'exon de chaque transcrit

nb_exons = dict()


# Pensez à adapter le chemin vers le fichier.

file_handler = open("../data/hg38_5k.gtf", "r")


# On utilise une boucle for pour
# lire des lignes renvoyées par l'objet
# file_handler

for line in file_handler:

    # On découpe la ligne sur le séparateur "\t"

    columns = line.split("\t")

    if columns[2] == "exon":

        # On substitue dans la colonne 9 n'importe quel caractère (.)
        # répété 0 fois ou plus (*) et
        # suivit de transcript_id, d'un espace et d'un double guillemet par
        # une chaine de caractères vide ("").

        tx_name = re.sub('.* transcript_id "', "", columns[8])

        # On substitue tout ce qui se trouve à droite du nom du transcrit: un guillemet suivit de
        # n'importe quel caractère (.) répété 0 fois ou plus (*) et d'un caractère retour à la ligne.
        tx_name = re.sub('".*\n', "", tx_name)


        # Si le transcript n'est pas connu dans le dictionnaire
        # On l'initialise à 1 pour la clef tx_name.
        if not nb_exons.get(tx_name):

            nb_exons[tx_name] = 1


        else:
            # Sinon on ajoute 1 à la valeur associée à la clef
            # tx_name

                nb_exons[tx_name] += 1

# Pour toutes les clefs du dictionnaire nb_exons

for tx_name in nb_exons:

    print(tx_name + "\t" + str(nb_exons[tx_name]))