Lecture Notes

Collected lectures notes in computational biology

Recherche reproductible et interfaces

Le biologiste moderne, nous l'avons vu doit souvent écrire de petits programmes. Ses principales préoccupations sont (dans l'ordre):

  • Le programme doit être correct
  • Les résultats doivent être facilement reproductibles
  • Le programme doit être pratique à utiliser

Vous remarquerez que la notion de vitesse d'exécution n'est même pas dans cette liste. En effet, c'est souvent un problème secondaire.

Assurer le caractère pratique de son code

De l'importance d'une bonne interface

L'interface est la partie de votre code avec lequel l'utilisateur interagit avec votre code (dans la recherche c'est souvent vous-même, éventuellement dans quelques mois/années, mais soyez sympa avec vous-même). Il est donc important qu'elle soit claire et pratique.

Pour des petites analyses, un Notebook Jupyter bien documenté (avec des blocs de textes en markdown) et commenté (avec des docstring et autres) est satisfaisant.

Dès que l'on commence à avoir une utilisation plus intensive (utilisation sur de nombreux jeux de données, très longs calculs)... il est bénéfique de sortir l'analyse du notebook. Quitte à ne l'utiliser que pour l'exploration interactive et la présentation (avec prise de notes) des résultats.

Je vous propose de parler de deux types d'interface. L'interface en ligne de commande (CLI, command line interface) et l'interface programmatique (Advanced Programming Interface, API). On ne couvrira pas les interface graphiques (Graphical user interface, GUI) parce que ce n'est pas si utile quand on développe un programme scientifique. On peut souvent se contenter de la GUI fournie par le notebook.

L'interface console (CLI) avec argparse

L'interface en ligne de commande est celle que vous invoquez quand vous lancer votre programme depuis le terminal.

Nous avons vu dans les exercices comment utiliser input pour demander à l'utilisateur de rentrer des valeurs. Cependant c'est une très mauvaise pratique. Une bonne interface en ligne de commande pour un logiciel scientifique n'est pas interactive, et permet que tous les paramètres soient passés par des arguments.

Cela à l'avantage de permettre la compatibilité avec les autres outils unix (cat, grep, awk, count...) mais surtout cela permet l'automatisation de votre programme. Comme le lancer en parallèle sur un grand nombre de jeux de données ou avec un grand nombre de paramètre, sur votre machine ou sur une grille de calcul.

In [2]:
# Tous les arguments passés au démarrage de
# l'interpréteur se trouvent dans la liste sys.argv.
import sys
sys.argv
Out[2]:
['/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py',
 '-f',
 '/run/user/1000/jupyter/kernel-62d83444-484e-4cf0-9e89-73aa417966f4.json']

Une solution est de parser sys.argv à la main. Heureusement Python à un module dans sa bibliothèque standard qui permet de créer facilement des interfaces en ligne de commande. Il s'agit de argparse. argparse possède une API orienté objet, un bon tutoriel et est très facile d'utilisation. Voici un exemple minimal:

#! /usr/bin/env python3
""" Un exemple d'interface. 
Mettre une docstring au début d'un fichier python définit la docstring du module associé."""
import argparse
import random

# Créer un objet Argument Parser. Le champ description est utilisé pour l'aide de votre programme.
parser = argparse.ArgumentParser(description='Generate a random sequence')

# Ajouter un nouvel argument -N, qui prend un entier, qui vaut 10 par défaut
# et qui doit être présent 1 ou 0 fois (vous pouvez utiliser *,+,? comme pour les regexp, ou un nombre)
parser.add_argument('-N', help='Sequence size', type=int, default=10, nargs="?")

# Parse les arguments présents dans sys.argv
args = parser.parse_args() # => Renvoie un objet de type Namespace

# On accède aux arguments avec la syntaxe habituelle des espaces de noms
print('> Random sequence of length {}'.format(args.N))
print(''.join([random.choice('ATGC') for _ in range(args.N)]))

Si vous mettez ceci dans un fichier sequence.py, voici le résultat:

guilhem@volvox:~/lectures/python$ ./sequences.py -N 3
> Random sequence of length 3
TTC

Vous pouvez constater que python fait la prise en charge des erreurs pour vous:

guilhem@volvox:~/lectures/python$ ./sequences.py -N "Blah"
usage: sequences.py [-h] [-N [N]]
sequences.py: error: argument -N: invalid int value: 'Blah'

Et créér même l'aide de votre programme automatiquement:

guilhem@volvox:~/lectures/python$ ./sequences.py --help
usage: sequences.py [-h] [-N [N]]

Generate a random sequence

optional arguments:
  -h, --help  show this help message and exit
  -N [N]      Sequence size
In [ ]:
# Exercice :
# Écrivez une interface CLI pour votre logiciel de traitement de séquence
#- `sequence.py -N 100` doit renvoyer une séquence aléatoire de longueur 100
#- `sequence.py -N 100 --fasta` doit le faire au format fasta
#- `sequence.py -N 100 --translate` doit convertir la séquence en protéines

# Avancé: accepter une séquence d'entrée dans un fichier si on utilise `--input fichier.fa` au lieu de `-N 100`.

L'interface programmatique (API) avec un module

Un second type d'interface, à laquelle vous avez déjà tous eu à faire est l'interface programmatique (API). Typiquement les bibliothèques comme numpy, matplotlib etc. définissent une interface que vous utilisez quand vous faites import numpy. Il existe aussi des API qui sont exposées sur le web et que l'on peut interroger depuis son propre code (exemple, celles du NCBI).

De manière générale une API est un ensemble de spécifications qui régissent la manière avec laquelle l'utilisateur-developpeur va interagir avec votre code. Les API des grandes bibliothèques en sont des composantes cruciales, et ne sont pas changées à la légère (et jamais sans changer le numéro de version), car cela pourrait casser le code des utilisateurs.

Au niveau, plus modeste, du chercheur, on peut vouloir réutiliser son code dans d'autre projet, ou publier une méthode avec une interface pratique pour son utilisation. Pour cela, on écrit un paquet python et on documente son API. Voici quelques notions importantes, comme toujours je vous renvoie à la documentation officielle pour plus d'information.

Quelques mots sur import

En python tout fichier contenant des instructions python (en général portant l'extension *.py) est un module et peut être importé. Voici un exemple:

#! /usr/bin/env python3
''' fichier module.py '''
def fonction_utile(a):
    ### Cette fonction fait partie de l'espace de nom du module.
    pass

print('module chargé !') # Ce code sera exécuté au premier import de ce module.
if __name__ == "__main__":
    ### Le code dans ce bloc sera exécuté seulement si on appelle ce
    ### module avec python3 module.py ou ./module.py
    ### mais pas si on import module
    fonction_utile(4)

Dans un autre fichier...

#! /usr/bin/env python3
''' fichier main.py '''
import module
# L'espace de nom du module est maintenant accessible avec la syntaxe "module.variable" ou "module.fonction"
module.fonction_utile(10)

Notez que les modules peuvent être regroupés en un sous-module en les plaçant dans le même dossier. Si on importe le nom du dossier, python cherchera un fichier avec le nom spécial __init__.py.

Ainsi, si on a cette arborescence:

./main.py
./module/ 
./module/__init__.py 
./module/sous_module.py

et ce fichier:

"""contenu de main.py"""
import module # Importe module/__init__.py
import module.sous_module # Importe module/sous_module.py

Mais où sont les modules ?

Quand on utilise la commande import blah, python cherche un fichier blah.py dans le répertoire local (donné par os.getcwd()) puis, s'il ne trouve pas il essaye successivement les chemins disponibles dans sys.path.

Exemple:  

>>> import sys
>>> print(sys.path)
['',
 '/usr/lib/python35.zip',
 '/usr/lib/python3.5',
 '/usr/lib/python3.5/plat-x86_64-linux-gnu',
 '/usr/lib/python3.5/lib-dynload',
 '/home/guilhem/.local/lib/python3.5/site-packages',
 '/usr/local/lib/python3.5/dist-packages',
 '/usr/lib/python3/dist-packages',
 '/usr/local/lib/python3.5/dist-packages/IPython/extensions'
]

Installer un module

Installer un module en python se fait avec la commande pip (Python Package Installer). La documentation complète est la référence, mais voici quelques trucs à retenir:

# Installe le paquet nomdupaquet depuis le Python Package Index [PyPI](https://pypi.org/) (necessite d'être root)
pip install nomdupaquet

# Installe le paquet dans .local/lib/ (pas besoin d'être root)
pip --user install nomdupaquet

# Installe un paquet depuis un dossier
pip install ./dossier/paquet

# Installe un paquet depuis un dossier SANS le copier dans le path
# ce qui permet de ne pas avoir à le réinstaller à chaque modification
# on parle de `developper mode`
pip install -e ./dossier/paquet

# Installe un paquet depuis une archive de sources
pip install ./paquet.tar.gz

# Installe un paquet depuis un dépôt gitlab|github bien configuré
pip install https://gitlab.com/<depot>/repository/archive.zip?ref=master

Un path temporaire pour chaque projet: virtualenv

Par défaut pip installe dans le path global du système, ce qui nécessite les droits d'administrateur. On peut installer en mode utilisateur, mais il reste un problème: plusieurs projets peuvent avoir besoin de plusieurs versions de modules différents (certains de vos collaborateurs sont sous numpy 10 et les autre sous numpy 11), ou vous voulez être sûr qu'une mise à jour intempestive de votre système ne changera pas la version des paquets.

L'utilitaire virtualenv permet de créer un "environnement virtuel" très facilement. C'est-à-dire une installation de python indépendante de votre système hôte, avec son propre path et ses propres paquets installés. La documentation sur virtualenv est, bien sûr, la référence mais voici les quelques trucs à retenir.

# Créer un virtualenv, avec son path, python et pip.
virtualenv mon_environnement 
# Cela créer un répertoire mon_environnement.
# On peut préciser la version de python avec l'option -p (voir man virtualenv)
# Activer le virtualenv
source mon_environnement/bin/activate
# Maintenant dans cette console, tout appel a python ou pip utilisera le python de cet environnement virtuel et non celui du système.
# Il faut y réinstaller vos paquets (avec pip).

# L'environnement est ouvert dans le shell courant.
# Pour le quitter sans fermer son shell (terminal, connection SSH) il faut faire:
deactivate

Écrire un paquet Python

Écrire un paquet python n'est pas très compliqué, à condition de s'y prendre avec méthode et bien suivre la documentation. Voici un exemple minimal.

./mon_paquet/
./mon_paquet/setup.py # Configure le paquet.
./mon_paquet/mon_paquet.py # Le module à empaqueter (pourrait être un dossier comme vu plus haut)

Le contenu de setup.py donne les informations sur le paquet:

import setuptools

long_description = """# Mon paquet il est beau...

...et en plus il expose de belles fonctions."""

setuptools.setup(
    name="mon_paquet",
    version="0.0.1",
    author="Ernest Dulm",
    author_email="example@ens.fr",
    description="Un bien beau paquet",
    long_description=long_description,
    long_description_content_type="text/markdown",
    url="https://github.com/pypa/sampleproject",
    packages=setuptools.find_packages(),
)

Vous pouvez maintenant installer ce paquet avec pip ./mon_paquet.

Ou créer une archive re-distribuable avec python3 setup.py sdist (var créer ./dist/mon_paquet-0.0.1.tar.gz)

Si vous souhaitez que votre paquet expose une commande pour votre terminal c'est dans le setup.py que ça se passe.

Exercice:

  • Écrire un paquet pour votre programme de traitement des séquences.
  • L'installer n mode développeur avec pip -e <chemin>
  • L'importer depuis un notebook se trouvant ailleurs dans votre home.

Assurer le caractère reproductible de son code

Très rapidement vous (ou vos collaborateurs) serez amené à ré-utiliser du code que vous avez écrit quelques mois/ans plus tôt, pour refaire une analyse, répondre à un reviewer. Pour que ceci se passe sans douleur, il est important de bien organiser son code.

Si tout votre code se trouve dans un notebook-spaghetti, où il faut exécuter les cellules dans une ordre précis pour qu'elles fonctionnent (parfois), vous allez avoir beaucoup de mal. Voici une série de conseils généraux:

  • Vérifiez ponctuellement que vos notebook s'exécutent sans erreurs si vous faites (Kernel->Restart and Run all). Toujours faire cette vérification avant d'arrêter de travailler sur un notebook.
  • Adhérez le plus possible à un style de programmation (comme celui de la PEP8 utilisé pour la bibliothèque standard et de nombreux paquets).
  • Un notebook doit servir à l'exploration interactive d'un problème: une fois que le code logique "critique" est stabilisé dans le notebook, mettez-le dans un module à part, et importez-le dans le notebook. Il sera plus facile à réutiliser, maintenir et tester.
  • Mettez votre code dans un gestionnaire de version comme git git et tenez vous y. Software carpentry propose un bon tutoriel git pour les scientifique.
  • Quand vous sauvegardez des résultats d'analyses ou de simulations, sauvegardez toujours les paramètres utilisés. Ça peut être aussi simple qu'un write dans un fichier metadata.txt. Faites en sorte qu'il soit impossible d'exporter de nouveaux résultats sans exporter les métadonnées en même temps.
  • Quand vous recevez des données de l'extérieur (collaborateur, article). Ne modifiez jamais les données "en place". Personnellement ma méthode est d'avoir un dossier data/raw/ qui contient exactement ce qu'on m'a envoyé et un notebook data/readme.ipynb qui documente les données (date, provenance, structure), teste si elles sont correctes (nombre d'entrées, pas de 0 aux mauvais endroits etc..) et les convertit dans un format plus simple à traiter (sauvegardé dans un dossier data/processed). Comme ça, même si les données changent (nouvelle expérience, etc...), il me suffit de changer cette partie.

Exercice:

  • Réécrivez votre programme sequence.py pour qu'il adhère a la PEP8 (Vous pouvez tester votre code avec pylint)
  • Voir le notebook "Ex3" pour un exemple de bonne pratique notebook/données.

Assurer le caractère correct de son code

Peut être la pire chose qui puisse arriver à un morceau de code en recherche est de pas avoir confiance aux résultats qu'il donne. S'il est très facile de détecter quand un programme plante (pas de résultats), il y a des cas où un programme peut vous donner des résultats (parfois subtilement faux).

Cela peut vous faire perdre beaucoup de temps. Comment s'éviter cela ?

De manière générale, si un programme est faux, on veut qu'il lève une exception plutôt que de continuer silencieusement. On veut qu'il la lève le plus tôt possible et de façon la plus explicite possible.

Les assertions et la programmation défensive

La programmation défensive consiste à vérifier un certain nombre de conditions connues tout au long du programme. Une méthode typique fait appel à la commande assert.

Une commande assert est suivie d'une expression. Si cette expression est correcte, il ne se passe rien. Si non, une exception AssertError est levée.

In [13]:
assert True #Il ne se passe rien
assert False, "On peut préciser un commentaire qui s'affichera dans la trace"
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-13-7b614d909674> in <module>()
      1 assert True #Il ne se passe rien
----> 2 assert False, "On peut préciser un commentaire qui s'affichera dans la trace"

AssertionError: On peut préciser un commentaire qui s'affichera dans la trace

Les assertions permettent de vérifier que:

  • les entrées sont correctes (on parle de pré-conditions)
  • les sorties sont correctes (on parle de post-conditions)
  • les conditions qui devraient toujours être vraies le reste (on parle d'invariants) et d'arrêter le programme si ce n'est pas le cas.
In [27]:
# Exemple: 
from collections import Counter
def gc_content(seq):
    """Calcule le GC-content de la séquence seq."""
    
    if not seq:
        return None
    
    # Ceci est une pré-condition, elle doit toujours être vraie en entrée de la fonction.
    assert all([base in 'ATGC' for base in seq]), 'Sequence contains other letters than ATGC'
    
    counter = Counter(seq)
    gc =  (counter['G']+counter['G'])/len(seq)
    
     # Ceci est une post-condition, elle doit être toujours vraie en sortie de la fonction.
    assert 0<=gc<=1
    
    return gc

gc_content('ATGCTTGACTTC')
gc_content('ATGCTTGACTTCU')
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-27-4efe75a2ed48> in <module>()
     19 
     20 gc_content('ATGCTTGACTTC')
---> 21 gc_content('ATGCTTGACTTCU')

<ipython-input-27-4efe75a2ed48> in gc_content(seq)
      8 
      9     # Ceci est une pré-condition, elle doit toujours être vraie en entrée de la fonction.
---> 10     assert all([base in 'ATGC' for base in seq]), 'Sequence contains other letters than ATGC'
     11 
     12     counter = Counter(seq)

AssertionError: Sequence contains other letters than ATGC

Notez que les assertions on un coût en temps. Mais il est possible de les ignorer en utilisant le paramètre -O de l'interpréteur python.

# On créer un module python très simple.
guilhem@volvox:~/lectures/python$ echo "assert False" > test.py

#Si on l'execute, l'assertion lève une erreur
guilhem@volvox:~/lectures/python$ python3 test.py 
Traceback (most recent call last):
  File "test.py", line 1, in <module>
    assert False
AssertionError

# Mais n'est pas testée avec l'option -O
guilhem@volvox:~/lectures/python$ python3 -O test.py
#(rien)

Les Tests

Les tests sont des fonctions particulières qui ne servent qu'à faire des assertions. On peut écrire des test pour l'ensemble d'un programme, ou pour des fonctions particulières (on parle de tests unitaires).

Une méthode de programmation (test-driven development) recommande de toujours écrire les tests avant d'écrire le code qui est sensé les passer. Il existe des bibliothèques spécialisées (ex. nose) qui facilitent l'écriture et l'exécution de tests.

In [40]:
def test_gc():
    """Un expemple de test simple"""
    # Les cas à tester sont appeler des "fixtures" dans la terminologie des tests
    assert gc_content('') == 0, 'Séquence vide'
    assert gc_content('ATGC') == 0.5, 'Séquence normale'
    assert gc_content('ATAT') == 0, 'Séquence sans G'
    assert gc_content('GGCC') == 1, 'Séquence de G'
    
    # On peut aussi tester si la bonne exception est levée.
    try:
        gc_content('U')
    except AssertionError:
        pass
    else:
        raise AssertionError('Mauvais caractères')
test_gc()

Exercice:

  • Écrivez des tests et des assertions pour votre programme de traitement des séquences.
  • Faites en sorte que les test se lançent à partir de la commande sequences.py --test