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

Par­tir de quel­conques don­nées de séquen­çage haut débit brutes pour arri­ver à une ana­lyse com­plète demande au mieux, une cer­taine pra­tique de ces tech­no­lo­gies. Dans bien des cas, on va alors mettre en place un pipe­line repo­sant sur tout un tas d'outils. Il fau­dra pro­ba­ble­ment des heures pour com­prendre les para­mètres de cha­cun d’entre eux et regrou­per l'ensemble en une jolie chaine de trai­te­ment pour arri­ver au résul­tat sou­hai­té (qui res­te­ra bien sûr à inter­pré­ter hum..).

Suite aux remarques sur un autre article concer­nant les points impor­tants à ne pas lou­per dans le cas d'analyse de carte de contact chro­mo­so­mique (Hi‑C), il m'a sem­blé inté­res­sant de venir vous pro­po­ser une petite démons­tra­tion d'un pipe­line pour ce type de don­nées. L'objectif de l'article est de vous don­ner une manière par­mi d'autres d'obtenir une carte de contact à par­tir de don­nées brutes (reads). De plus, je mon­tre­rai briè­ve­ment quelques exemples de scripts pour visua­li­ser les cartes pro­duites en sor­tie du pipe­line.

Plu­tôt que de vous expo­ser une suc­ces­sion de scripts faits par mes soins, j’ai choi­si de vous pré­sen­ter HiC-Pro : un pipe­line simple à uti­li­ser et donc idéal pour les débu­tants devant ana­ly­ser des don­nées Hi‑C sans faire d'erreur. Je ne pré­tends en aucun cas que ce pipe­line est le meilleur, ni le plus effi­cace mais plu­tôt une option abor­dable pour com­men­cer. Point impor­tant, l'ensemble des codes pro­po­sés aujourd'hui sont pen­sés pour fonc­tion­ner sous GNU/​Linux, je ne connais à ce jour, aucun équi­valent sous Win­dows.

Utilisation du pipeline

Installation de la bête

Pour uti­li­ser ce pipe­line, il faut l'installer. De ce coté là rien de très dif­fi­cile, tout est dans le tuto­riel en ligne qui est bien plus com­plet que mes propres ins­truc­tions.

Petit rap­pel de ce qu'il faut avoir sur sa machine avant de com­men­cer :

ensuite il suf­fit de taper quelques lignes de com­mande et c'est par­ti !

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

Pour cette petite ses­sion tuto­riel, je vous pro­pose de faire fonc­tion­ner le pipe­line sur les don­nées de Bing Ren. J'ai choi­si d'analyser les cartes pro­ve­nant de cel­lules sur­ré­nales. Ces don­nées sont dis­po­nibles en ligne et cor­res­pondent aux réfé­rences SRA allant de SRR4271980 à SRR4271983.  Pour pou­voir obte­nir ces don­nées, il faut d’abord les télé­char­ger à par­tir de la page de la publi­ca­tion, celle-ci nous indique l'ensemble des fichiers à télé­char­ger pour construire notre carte.

Cepen­dant ces don­nées sont au for­mat SRA (comme beau­coup de don­nées de séquen­çage dis­po­nibles en ligne) et HiC-Pro ne sup­porte que le for­mat FASTQ ! Il faut donc conver­tir les fichiers au for­mat SRA vers le for­mat FASTQ, ce qui est pos­sible grâce à sra­tool­kit. Les com­mandes ci-des­sous télé­chargent et conver­tissent direc­te­ment les don­nées :

(Exemple type : ici je télé­charge les don­nées et conver­tis en for­mat FASTQ en direct)

sratoolkit.2.8.0‑centos_linux64/bin/prefetch.2.8.0 SRR4271981
sratoolkit.2.8.0‑centos_linux64/bin/fastq-dump –split-files SRR4271981

sratoolkit.2.8.0‑centos_linux64/bin/prefetch.2.8.0 SRR4271982
sratoolkit.2.8.0‑centos_linux64/bin/fastq-dump –split-files SRR4271982

sratoolkit.2.8.0‑centos_linux64/bin/prefetch.2.8.0 SRR4271983
sratoolkit.2.8.0‑centos_linux64/bin/fastq-dump –split-files SRR4271983

Utilisation : le fichier de configuration du pipeline

Main­te­nant que les don­nées sont télé­char­gées et sau­ve­gar­dées au for­mat FASTQ, il est temps de lan­cer le pipe­line !

Il nous faut alors édi­ter le fichier de confi­gu­ra­tion pro­po­sé par le pipe­line pour l'adapter à nos don­nées. À prio­ri ce n'est pas très long. Il faut avant tout véri­fier le génome de réfé­rence et l'enzyme de res­tric­tion pour que les reads soient assi­gnés aux bonnes coor­don­nées dans la carte de contact. Ici le génome à été digé­ré par Hin­dIII, nous adap­te­rons pour avoir dans la par­tie du fichier dédiée à la diges­tion les lignes sui­vantes :

Pen­sez éga­le­ment à bien indi­quer la réso­lu­tion à laquelle vous vou­lez voir vos cartes , c'est indi­qué par le para­mètre BIN_​SIZE.

Ici le pipe­line va donc géné­rer une carte de réso­lu­tion 10kb et une de 100kb. Le for­mat est ici 'upper', ce qui signi­fie que seule­ment la par­tie supé­rieure gauche de la carte sera repré­sen­tée. Cepen­dant, la carte étant symé­trique, nous pour­rons déduire l'autre par­tie faci­le­ment.

Appel du pipeline

Le pipe­line s’exécute inté­gra­le­ment en ne lan­çant qu'une seule ligne de com­mande. A par­tir des don­nées FASTQ conte­nues ici dans le réper­toire /users/­Do­cu­ments/­Hi‑C/­Bing-ren/a­dre­nal/­fast­q/­da­ta/, le pipe­line va réa­li­ser tous les trai­te­ments néces­saires jusqu'à l'obtention d'une carte qu'il sau­ve­gar­de­ra en for­mat texte. Ces résul­tats seront sto­ckés dans le réper­toire /users/­Do­cu­ments/­Hi‑C/­Bing-ren/a­dre­nal/­re­sul­tat/.

(Exemple type de ce que j'ai tapé, et qui n'est qu'une tra­duc­tion mal­adroite de la docu­men­ta­tion déjà exis­tante)

Petits détails tech­niques :

-Un fichier conte­nant la taille des chro­mo­somes est ren­sei­gné dans le fichier de confi­gu­ra­tion. L'ordre des chro­mo­somes dans ce fichier est conser­vé pour les cartes obte­nues. Faites donc atten­tion de ne pas prendre un chro­mo­some pour un autre.

-En lan­çant Hi‑C Pro, j'ai pas­sé le réper­toire /users/­Do­cu­ments/­Hi‑C/­Bing-ren/a­dre­nal/­fast­q/ comme source (option ‑i, pour input). Pour­tant, mes don­nées sont dans le sous-réper­toire /users/­Do­cu­ments/­Hi‑C/­Bing-ren/a­dre­nal/­fast­q/­da­ta/. HiC-Pro consi­dère chaque sous-réper­toire comme un échan­tillon a trai­ter indé­pen­dam­ment. Si vous n'avez qu'un seul échan­tillon, vous devez tout de même pla­cer les fichiers dans un sous-réper­toire sans quoi HiC-Pro ne fonc­tion­ne­ra pas. L'information est docu­men­tée dans le manuel de confi­gu­ra­tion de HiC-Pro, mais il sem­ble­rait qu'elle ne soit pas évi­dente à trou­ver (pour en avoir par­lé sur notre mer­veilleux chan #bioin­fo-fr, nous sommes plu­sieurs à être pas­sés sur ce point de la docu­men­ta­tion sans le remar­quer).

Analyse des résultats

Comprendre les fichiers de sorties

A sa sor­tie, Hi‑C pro retourne deux types de cartes : la carte sans trai­te­ment, et la carte nor­ma­li­sée selon la méthode de Ima­kaev, ré-implé­men­tée dans ce pipe­line (plus d'explication dans mon article pré­cé­dent). Les cartes sont retour­nées dans un for­mat texte un peu par­ti­cu­lier. Pour une carte, nous avons deux fichiers : Un fichier tabu­lé de type 'nom.matrix' don­nant pour chaque couple de coordonnés(i,j) dans la matrice son nombre de contacts et un fichier BED don­nant la cor­res­pon­dance entre une coor­don­née chro­mo­so­mique et son indice dans la matrice.

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

On com­prend 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 éga­le­ment. Dans les cases (587,2) et (2,587), il y a 6 contacts etc..

Exemple de fichier bed :

On com­prend alors que la réso­lu­tion de la matrice est de 100 000 paire de base (ou 100kb) par bin dans la matrice (sur­prise, comme nous l'avons para­mé­tré !). Nous pou­vons éga­le­ment aisé­ment trou­ver la taille totale de la matrice qui cor­res­pond tout sim­ple­ment au nombre de lignes de ce fichier mais aus­si à la der­nière valeur de la der­nière ligne du fichier, ici 30971.

Conversion de la carte

Pour pou­voir uti­li­ser cette carte, j'ai fait un petit script python3 me per­met­tant de la conver­tir en à peu près ce que je veux (en véri­té, uti­li­sant MATLAB dans l'équipe, il me sert sur­tout à avoir un fichier .hdf5). Voi­ci donc quelques lignes que je vous pro­pose pour obte­nir une carte en for­mat Num­Py array (qu'il reste à sau­ve­gar­der en for­mat hdf5 ):

J'ai main­te­nant des fonc­tions pour conver­tir ma matrice tex­tuelle en for­mat num­py qu'il me faut alors sau­ve­gar­der. Pour cela j'ai choi­si d'utiliser le for­mat hdf5, com­pa­tible avec plu­sieurs lan­gages comme python ou mat­lab. J'appelle alors mes fonc­tions décrites au des­sus avec le code sui­vant :

Quelques lignes de code pour visualiser

Je prie les puristes, ceux qui ne jurent que par les lan­gages libres de m'excuser pour la suite… mais dans cette sec­tion, je vais être ame­né à… uti­li­ser MATLAB. Mon équipe uti­lise prin­ci­pa­le­ment ce lan­gage, c'est pour­quoi mes scripts de visua­li­sa­tion sont écrits avec. D'autres articles vien­dront pro­ba­ble­ment com­plé­ter l'analyse de don­nées Hi‑C avec des lan­gages, outils et sys­tèmes dif­fé­rents de ceux pré­sen­tés ici.

À l'étape pré­cé­dente, j'ai obte­nu une carte hdf5 lisible aus­si bien en Python qu'en MATLAB. Voyons voir à quoi elle res­semble !

Pour cela rien de plus simple, les lignes sui­vantes chargent une carte au for­mat hdf5 et affichent le chro­mo­some 1 à une réso­lu­tion de 100kb à l'échelle log10. S'il faut des exemples de prise en main pour la suite des opé­ra­tions, n'hésitez pas a me l'indiquer en com­men­taire, je ferai alors un autre article ! Mal­heu­reu­se­ment, après recherche, cette par­tie ne fonc­tionne pas direc­te­ment en Sci­lab (un équi­valent libre de MATLAB), si vous sou­hai­ter uti­li­ser  Sci­lab, il faut rem­pla­cer la fonc­tion h5read par par h5open, et ima­gesc par imshow.

Et voi­là, si vous ne saviez pas trop quoi faire avec des don­nées brutes, vous avez une idée de piste sur la façon de faire. Si vous ne sou­hai­tez pas uti­li­ser HiC-Pro, il existe d'autres outils faciles à prendre en main, tels que Jui­cer. Je ne peux que vous recom­man­der éga­le­ment de voir ce que font ces outils pour construire vous-même votre propre pipe­line si vos don­né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 mer­ci au relec­teurs de ce petit article, ayant consi­dé­ra­ble­ment aidé à rendre cet article meilleur qu'il ne l'était à la base : June Sal­lou , Mat Blum et Hed­jour. Mer­ci éga­le­ment à Yoann M, admi­nis­tra­teur de la semaine !