Objectif

On se propose de d'implémenter une commande tx_len capable de renvoyer pour chaque transcript d'un fichier gtf la couverture génomique (i.e la taille exons et introns compris). Le fichier donné en entrée par l'utilisateur n'est pas forcément trié par coordonnées génomiques. Etant donné que les coordonnées génomiques sont toujours données en fonction du brin plus du génome, cela revient à rechercher dans le fichier GTF le START le plus petit et le END le pus grand pour chaque transcrit.

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 "tx_len.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

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.

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 son 'start' (valeur).

Créer un dictionnaire permettant de mémoriser le nom du transcrit (clef) et son 'end' (valeur).

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 (choisir l'un ou l'autre):

Initialiser tous les dictionnaires

Sinon:

Si la valeur de 'start' courante pour le transcrit est inférieure à la valeur connu pour ce transcrit:

Assigner cette valeur au transcrit courant.

Si la valeur de 'end' courante pour le transcrit est supérieure à la valeur connu pour ce transcrit:

Assigner cette valeur au transcrit courant.

Parcourir les clefs de l'un des dictionnaire:

Imprimer transcrit et sa valeur end-start + 1

# -*- 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: tx_len.py
This program is intended to compute the genomic coverage of transcripts (exon
 + intron length).
"""
# 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 les coordonnées de
# départ (start) des transcrits

transcript_starts = dict()

# Ce dictionnaire contiendra les coordonnées de
# fin (end) des transcrits
transcript_ends = dict()


# 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")


# 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")

    # On récupère le start et le end du transcrit courant

    start_current = int(columns[3]) # attention au type
    end_current = int(columns[4]) # attention au type



    # 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 dans transcript_starts
    # il n'est pas non plus dans les autres dictionnaires.
    # On les initialise.
    if not transcript_starts.get(tx_name):

        transcript_starts[tx_name] = start_current
        transcript_ends[tx_name] = end_current

    else:
        # Sinon le start est plus petit que le start connu,
        # on le conserve.
        if start_current < transcript_starts[tx_name]:

            transcript_starts[tx_name] = start_current

        # Si le end est plus grand que celui connu pour ce transcrit,
        # On le conserve.
        if end_current > transcript_ends[tx_name]:

            transcript_ends[tx_name] = end_current

# Pour toutes les clefs du dictionnaire transcript_starts

for tx_name in transcript_starts:

    print(tx_name + "\t" + str(transcript_ends[tx_name] - transcript_starts[tx_name] + 1))

L'ordinateur nous donne une note

Bonne nouvelle, plus besoin d'examens, l'ordinateur va pourvoir noter la syntaxe utilisée dans notre programme (et il sera on ne peut plus impartial !). Ce n'est pas encore parfait mais ça devrait le devenir dans les cours suivants.

Installez pylint avec la commande suivante:

pip3 install pylint --user

Ensuite demander à pylint de noter votre programme. Pour se faire taper la commande suivante:

pylint tex_len.py # adaptez le chemin/nom du fichier éventuellement.

Notez que nous ne pourrons pas avoir la note maximale car nous n'utilisons pas de fonction à ce stade. Une prochaine fois espérons le.