Didacticiel :
Vous ne savez pas comment analyser vos données Hi-C? Exemple d'utilisation de HiC-pro

Partir de quelconques données de séquençage haut débit brutes pour arriver à une analyse complète demande au mieux, une certaine pratique de ces technologies. Dans bien des cas, on va alors mettre en place un pipeline reposant sur tout un tas d'outils. Il faudra probablement des heures pour comprendre les paramètres de chacun d’entre eux et regrouper l'ensemble en une jolie chaine de traitement pour arriver au résultat souhaité (qui restera bien sûr à interpréter hum..).

Suite aux remarques sur un autre article concernant les points importants à ne pas louper dans le cas d'analyse de carte de contact chromosomique (Hi-C), il m'a semblé intéressant de venir vous proposer une petite démonstration d'un pipeline pour ce type de données. L'objectif de l'article est de vous donner une manière parmi d'autres d'obtenir une carte de contact à partir de données brutes (reads). De plus, je montrerai brièvement quelques exemples de scripts pour visualiser les cartes produites en sortie du pipeline.

Plutôt que de vous exposer une succession de scripts faits par mes soins, j’ai choisi de vous présenter HiC-Pro : un pipeline simple à utiliser et donc idéal pour les débutants devant analyser des données Hi-C sans faire d'erreur. Je ne prétends en aucun cas que ce pipeline est le meilleur, ni le plus efficace mais plutôt une option abordable pour commencer. Point important, l'ensemble des codes proposés aujourd'hui sont pensés pour fonctionner sous Linux, je ne connais à ce jour, aucun équivalent sous windows.

Utilisation du pipeline

Installation de la bête

Pour utiliser ce pipeline, il faut l'installer. De ce coté là rien de très difficile, tout est dans le tutoriel en ligne qui est bien plus complet que mes propres instructions.

Petit rappel de ce qu'il faut avoir sur sa machine avant de commencer :

ensuite il suffit de taper quelques lignes de commande et c'est parti !

Téléchargement d'un exemple de jeu de données

Pour cette petite session tutoriel, je vous propose de faire fonctionner le pipeline sur les données de Bing Ren. J'ai choisi d'analyser les cartes provenant de cellules surrénales. Ces données sont disponibles en ligne et correspondent aux références SRA allant de SRR4271980 à SRR4271983.  Pour pouvoir obtenir ces données, il faut d’abord les télécharger à partir de la page de la publication, celle-ci nous indique l'ensemble des fichiers à télécharger pour construire notre carte.

Cependant ces données sont au format SRA (comme beaucoup de données de séquençage disponibles en ligne) et HiC-Pro ne supporte que le format FASTQ ! Il faut donc convertir les fichiers au format SRA vers le format FASTQ, ce qui est possible grâce à sratoolkit. Les commandes ci-dessous téléchargent et convertissent directement les données :

(Exemple type : ici je télécharge les données et convertis en format FASTQ en direct)

Utilisation : le fichier de configuration du pipeline

Maintenant que les données sont téléchargées et sauvegardées au format FASTQ, il est temps de lancer le pipeline !

Il nous faut alors éditer le fichier de configuration proposé par le pipeline pour l'adapter à nos données. À priori ce n'est pas très long. Il faut avant tout vérifier le génome de référence et l'enzyme de restriction pour que les reads soient assignés aux bonnes coordonnées dans la carte de contact. Ici le génome à été digéré par HindIII, nous adapterons pour avoir dans la partie du fichier dédiée à la digestion les lignes suivantes :

Pensez également à bien indiquer la résolution à laquelle vous voulez voir vos cartes , c'est indiqué par le paramètre BIN_SIZE.

Ici le pipeline va donc générer une carte de résolution 10kb et une de 100kb. Le format est ici 'upper', ce qui signifie que seulement la partie supérieure gauche de la carte sera représentée. Cependant, la carte étant symétrique, nous pourrons déduire l'autre partie facilement.

Appel du pipeline

Le pipeline s’exécute intégralement en ne lançant qu'une seule ligne de commande. A partir des données FASTQ contenues ici dans le répertoire /users/Documents/Hi-C/Bing-ren/adrenal/fastq/data/, le pipeline va réaliser tous les traitements nécessaires jusqu'à l'obtention d'une carte qu'il sauvegardera en format texte. Ces résultats seront stockés dans le répertoire /users/Documents/Hi-C/Bing-ren/adrenal/resultat/.

./HiC-Pro_2.7.9/bin/HiC-Pro -i ~/Documents/Hi-C/Bing-ren/adrenal/fastq/ -o ~/Documents/Hi-C/Bing-ren/adrenal/resultat/ -c config-hicpro-bingren-human.txt

(Exemple type de ce que j'ai tapé, et qui n'est qu'une traduction maladroite de la documentation déjà existante)

Petits détails techniques :

-Un fichier contenant la taille des chromosomes est renseigné dans le fichier de configuration. L'ordre des chromosomes dans ce fichier est conservé pour les cartes obtenues. Faites donc attention de ne pas prendre un chromosome pour un autre.

-En lançant Hi-C Pro, j'ai passé le répertoire /users/Documents/Hi-C/Bing-ren/adrenal/fastq/ comme source (option -i, pour input). Pourtant, mes données sont dans le sous-répertoire /users/Documents/Hi-C/Bing-ren/adrenal/fastq/data/. HiC-Pro considère chaque sous-répertoire comme un échantillon a traiter indépendamment. Si vous n'avez qu'un seul échantillon, vous devez tout de même placer les fichiers dans un sous-répertoire sans quoi HiC-Pro ne fonctionnera pas. L'information est documentée dans le manuel de configuration de HiC-Pro, mais il semblerait qu'elle ne soit pas évidente à trouver (pour en avoir parlé sur notre merveilleux chan #bioinfo-fr, nous sommes plusieurs à être passés sur ce point de la documentation sans le remarquer).

Analyse des résultats

Comprendre les fichiers de sorties

A sa sortie, Hi-C pro retourne deux types de cartes : la carte sans traitement, et la carte normalisée selon la méthode de Imakaev, ré-implémentée dans ce pipeline (plus d'explication dans mon article précédent). Les cartes sont retournées dans un format texte un peu particulier. Pour une carte, nous avons deux fichiers : Un fichier tabulé de type 'nom.matrix' donnant pour chaque couple de coordonnés(i,j) dans la matrice son nombre de contacts et un fichier BED donnant la correspondance entre une coordonnée chromosomique et son indice dans la matrice.

Exemple de fichier matrix (avec i, j et le nombre de contact, par ligne) :

On comprend donc que dans la case (1,473) et la case (473,1) de la matrice, il y a 1 contact, dans la case (1,700) et (700,1) 1 également. Dans les cases (587,2) et (2,587), il y a 6 contacts etc..

Exemple de fichier bed :

On comprend alors que la résolution de la matrice est de 100 000 paire de base (ou 100kb) par bin dans la matrice (surprise, comme nous l'avons paramétré !). Nous pouvons également aisément trouver la taille totale de la matrice qui correspond tout simplement au nombre de lignes de ce fichier mais aussi à la dernière valeur de la dernière ligne du fichier, ici 30971.

Conversion de la carte

Pour pouvoir utiliser cette carte, j'ai fait un petit script python3 me permettant de la convertir en à peu près ce que je veux (en vérité, utilisant MATLAB dans l'équipe, il me sert surtout à avoir un fichier .hdf5). Voici donc quelques lignes que je vous propose pour obtenir une carte en format NumPy array (qu'il reste à sauvegarder en format hdf5 ):

J'ai maintenant des fonctions pour convertir ma matrice textuelle en format numpy qu'il me faut alors sauvegarder. Pour cela j'ai choisi d'utiliser le format hdf5, compatible avec plusieurs langages comme python ou matlab. J'appelle alors mes fonctions décrites au dessus avec le code suivant :

Quelques lignes de code pour visualiser

Je prie les puristes, ceux qui ne jurent que par les langages libres de m'excuser pour la suite... mais dans cette section, je vais être amené à... utiliser MATLAB. Mon équipe utilise principalement ce langage, c'est pourquoi mes scripts de visualisation sont écrits avec. D'autres articles viendront probablement compléter l'analyse de données Hi-C avec des langages, outils et systèmes différents de ceux présentés ici.

À l'étape précédente, j'ai obtenu une carte hdf5 lisible aussi bien en Python qu'en MATLAB. Voyons voir à quoi elle ressemble !

Pour cela rien de plus simple, les lignes suivantes chargent une carte au format hdf5 et affichent le chromosome 1 à une résolution de 100kb à l'échelle log10. S'il faut des exemples de prise en main pour la suite des opérations, n'hésitez pas a me l'indiquer en commentaire, je ferai alors un autre article! Malheureusement, après recherche, cette partie ne fonctionne pas directement en Scilab (un équivalent libre de MATLAB), si vous souhaiter utiliser  Scilab, il faut remplacer la fonction h5read par par h5open, et imagesc par imshow.

Résultat de l'affichage de la carte du chromosome 1.

 

Et voilà, si vous ne saviez pas trop quoi faire avec des données brutes, vous avez une idée de piste sur la façon de faire. Si vous ne souhaitez pas utiliser HiC-Pro, il existe d'autres outils faciles à prendre en main, tels que Juicer. Je ne peux que vous recommander également de voir ce que font ces outils pour construire vous-même votre propre pipeline si vos données sortent un peu de la norme. Avec la série d'articles déjà écrits sur le site en quelques années, vous avez toutes les pistes pour !

Un grand merci au relecteurs de ce petit article, ayant considérablement aidé à rendre cet article meilleur qu'il ne l'était à la base : June Sallou , Mat Blum et Hedjour. Merci également à Yoann M, administrateur de la semaine!

  • À propos de
  • Bioinformaticien de formation ayant complété son cursus par un master de machine learning appliqué aux langues à Nantes. Actuellement en Doctorat au laboratoire physique théorique de la matière condensé à Paris, ma thèse a pour sujet les éléments répétés et l'organisation des génomes eucaryotes. J'aime python, le machine learning et les techniques d'études des chromosomes.

Catégorie: Didacticiel | Tags: , , , ,

Laisser un commentaire