On se propose de d'implémenter une commande count capable de compter le nombre non-redondant de transcrits dans un fichier gtf. On illustre, par ailleurs, ci-dessous l'importance de bien choisir les structures stockant les données, le choix d'une liste ou d'un dictionnaire ayant un effet remarquable sur le temps d'éxécution.
# -*- coding: utf8 -*-Sauvegardez le fichier sous le nom "count.py".
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
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";
Nous devrons à un moment ou à un autre pouvoir conserver quelque part la liste nous redondante des transcripts. On propose de stocker l'information en utilisant trois objets differents et, donc, d'implémenter trois solutions:
- L'une utilisant une liste dans laquelle on stockera la liste non-redondante des transcripts.
- L'autre utilisant un dictionnaire dont les clefs seront les transcripts, et les valeurs '1' si le transcript a été vu.
- L'une utilisant un ensemble (set) qui constitue un ensemble non ordonné d'éléments uniques.
Nous verrons par la suite que ce choix a un impact très important sur la vitesse d'exécution de notre programme. La première version (liste) prenant plusieurs minutes contre quelques secondes pour la deuxième et troisième version (dictionnaire et ensemble).
Le pseudo -code ci-dessous illustre la stratégie du programme:
- Créer une liste vide, list_of_tx_name.
- 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 nom du transcript de la ligne courante n'est pas dans la liste list_of_tx_name:
- l'ajouter à list_of_tx_name.
- Imprimer la taille du tableau.
L'implémentation de cette solution en Python se traduit par le code suivant:
# -*- coding: utf-8 -*-
# La ligne précédente permet d'utiliser des caractères accentués
# dans les commentaires
# On importe le module re
# afin de pouvoir rechercher
# une expression régulière dans une
# chaîne de caractères.
import re
# Pensez à adapter le chemin vers le fichier.
# Ici, il est dans le dossier data présent dans le
# répertoire supérieur
file_handler = open("../data/hg38_5k.gtf", "r")
# Cette liste contiendra l'ensemble des transcrits (tx)
list_of_tx_name = list()
# On utilise une boucle for pour
# lire des lignes renvoyées par l'objet
# file_handler
for line in file_handler:
# On substitue dans 'line' 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 "', "", line)
# 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 transcrit n'est pas dans la liste, on l'ajoute
if tx_name not in list_of_tx_name:
list_of_tx_name = list_of_tx_name + [tx_name] # ou tx_list.append(tx_name)
# Après être sorti de la boucle
# on imprime le nombre de transcrits.
# len(tx_list) est un entier (int). Pour pouvoir le concatener
# avec une chaîne de caractères (le message), il faut le transformer
# en chaîne de caractères (str()).
print("Le nombre de transcrits est: " + str(len(list_of_tx_name))) # n = 5000
On voit de manière évidente que le temps d'exécution est long (8.592s) pour un ensemble d'opérations qui semble relativement simple. Quel est le problème ? Ici nous utilisons une liste et quand celle ci s'étend, les temps de recherche (tx_name not in list_of_tx_name) deviennent relativement importants. En effet, pour ce test, il faut comparer une chaîne de caractères à toutes les autres. La solution est ici d'utiliser plutôt un dictionnaire.
Quand on utilise un dictionnaire, les cases du tableau (les clefs) sont associées à des valeurs. La recherche du clef dans la liste des clefs est optimisée afin de ne pas avoir à comparer lettre à lettre une chaîne de caractère à toutes les autres. Le pseudo -code ci-dessous illustre la stratégie du programme:
- Créer un dictionnaire vide, dict_of_tx_name.
- 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 nom du transcript de la ligne courante n'est pas défini dans dict_of_tx_name:
- l'ajouter à dict_of_tx_name.
- Imprimer la taille de dict_of_tx_name (i.e le nombre de clefs).
L'implémentation de cette solution en Python se traduit par le code suivant (les changements apportés sont notés par ***):
# -*- coding: utf-8 -*-
# La ligne précédente permet d'utiliser des caractères accentués
# dans les commentaires
# On importe le module re
# afin de pouvoir rechercher
# une expression régulière dans une
# chaîne de caractères.
import re
# Pensez à adapter le chemin vers le fichier.
# Ici, il est dans le dossier data présent dans le
# répertoire supérieur
file_handler = open("../data/hg38_5k.gtf", "r")
# ***
# Ce dictionnaire contiendra l'ensemble des transcrits (tx)
dict_of_tx_name = dict()
# On utilise une boucle for pour
# lire des lignes renvoyées par l'objet
# file_handler
for line in file_handler:
# On substitue dans 'line' 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 "', "", line)
# 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 transcrit n'est pas dans la liste, on l'ajoute
# ***
if tx_name not in dict_of_tx_name:
# ***
# On déclare une nouvelle clef dans le dictionnaire
# et on place sa valeur à 1.
dict_of_tx_name[tx_name] = 1
# Après être sorti de la boucle
# on imprime le nombre de transcrits.
# len(tx_list) est un entier (int) qui correspond au nombre de clef du dictionnaire. Pour pouvoir le concatener
# avec une chaîne de caractères (le message), il faut le transformer
# en chaîne de caractères (str()).
print("Le nombre de transcrits est: " + str(len(dict_of_tx_name))) # n = 5000
Alors que le résultat est identique, le temps d'exécution tombe à 0.627s. Soit un gain de temps important.
Le set est lui aussi un objet très optimisé pour les recherche. Il ne contient que des éléments uniques qui sont non ordonnés. Dans notre cas nous sommes intéressés par la liste non-redondante et l'ordre dans cette liste nous importe peu. Il apparait comme une solution très intéressante. Il occupera par ailleurs moins d'espace mémoire que le dictionnaire.
- Créer un ensemble vide, set_of_tx_name.
- 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 nom du transcript de la ligne courante n'est pas défini dans set_of_tx_name:
- l'ajouter à set_of_tx_name (méthode add).
- Imprimer la taille de set_of_tx_name (i.e le nombre de clefs).
L'implémentation de cette solution en Python se traduit par le code suivant (les changements apportés sont notés par ***):
# -*- coding: utf-8 -*-
# La ligne précédente permet d'utiliser des caractères accentués
# dans les commentaires
# On importe le module re
# afin de pouvoir rechercher
# une expression régulière dans une
# chaîne de caractères.
import re
# Pensez à adapter le chemin vers le fichier.
# Ici, il est dans le dossier data présent dans le
# répertoire supérieur
file_handler = open("../data/hg38_5k.gtf", "r")
# ***
# Cet ensemble contiendra l'ensemble des transcrits (tx)
set_of_tx_name = set()
# On utilise une boucle for pour
# lire des lignes renvoyées par l'objet
# file_handler
for line in file_handler:
# On substitue dans 'line' 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 "', "", line)
# 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)
set_of_tx_name.add(tx_name) # on utilise la méthode add.
# Après être sorti de la boucle
# on imprime le nombre de transcrits.
# len(tx_list) est un entier (int) qui correspond au nombre de clef du dictionnaire. Pour pouvoir le concatener
# avec une chaîne de caractères (le message), il faut le transformer
# en chaîne de caractères (str()).
print("Le nombre de transcrits est: " + str(len(set_of_tx_name))) # n = 5000
La vitesse ici est très proche du dictionnaire cependant cette solution est sans doute la meilleure car le volume occupé en mémoire par un ensemble est réduit par rapport à un dictionnaire.