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 :

bowtie2 servira à l'alignement des reads
Python (>= 2.7) avec pysam (>=0.8.3), bx(>=0.5.0), numpy(>=1.8.2) et scipy(>=0.15.1)
R avec RColorBrewer, ggplot2
g++ : pour compiler le code
Samtools (>=0.1.19)

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

tar -zxvf HiC-Pro-master.tar.gz
cd HiC-Pro-master
## Edit config-install.txt file if necessary
make configure
make install

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)

sratoolkit.2.8.0-centos_linux64/bin/prefetch.2.8.0 SRR4271980
sratoolkit.2.8.0-centos_linux64/bin/fastq-dump --split-files SRR4271980

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

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 :

#######################################################################
## Digestion Hi-C
#######################################################################

GENOME_FRAGMENT = HindIII_resfrag_hg19.bed
LIGATION_SITE = AAGCTAGCTT
MIN_FRAG_SIZE = #si laissé vide, des valeurs par défaut seront utilisées
MAX_FRAG_SIZE =
MIN_INSERT_SIZE =
MAX_INSERT_SIZE =

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

#######################################################################
## Contact Maps
#######################################################################

BIN_SIZE = 10000 100000
MATRIX_FORMAT = upper

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) :

1    473    1
1    700    1
1    804    1
1    859    1
1    1194    2
1    1889    1
1    1956    1
1    2263    1
1    2443    1
1    2497    1
...
587 2 6
...

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 :

chr1    0    100000    1
chr1    100000    200000    2
chr1    200000    300000    3
chr1    300000    400000    4
chr1    400000    500000    5
chr1    500000    600000    6
chr1    600000    700000    7
chr1    700000    800000    8
chr1    800000    900000    9
chr1    900000    1000000    10

...

chrY    59300000    59373566    30970
chrM    0    16571    30971

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 ):

#in : _abs.bedfile
#out : dict chr=[begin,end] pour chaque chromosome de la matrice, contient egalement la taille totale de la carte
def loadabsdatafile(filein):
    fin=open(filein,"r")
    d={}
    d["Total"]=0
    l=fin.readline()
    while l:
        ls=l.split()
        if d.__contains__(ls[0]):
            d[ls[0]][1]=float(ls[3])
        else:
            d[ls[0]]=[float(ls[3]),float(ls[3])]
        d["Total"]+=1
        l=fin.readline()
    fin.close()
    return d

#in : un fichier de matrice, sa taille
#out : la matrice en numpy array
def loadmatrix(filein,sizemat):
    print(sizemat)
    mat=np.zeros((sizemat,sizemat))
    fin=open(filein,"r")
    l=fin.readline()
    while l:
        ls=l.split()
        i=float(ls[0])-1
        j=float(ls[1])-1
        #print(i,j)
        v=float(ls[2])
        mat[i,j]=v
        mat[j,i]=v
        l=fin.readline()
    fin.close()
    return mat

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 :

d=loadabsdatafile("adrenal_100000_abs.bed")
print(d) #histoire de voir l'ensemble des chromosomes disponibles dans la carte
mat=loadmatrix("adrenal_100000.matrix",d["Total"])
print(mat.shape) # taille de la matrice complete
fh5 = h5py.File('adrenal_full.hdf5', "w")
fh5['data'] = mat
fh5.close()

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.

Hicmat = h5read('adrenal_full.hdf5','/data');
%chagement du chromosome 1 pour l'afficher.
Hicmatchr1=Hicmat(1:2493,1:2493);
figure,imagesc(Hicmatchr1);
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!



Pour continuer la lecture :


Commentaires

2 réponses à “Vous ne savez pas comment analyser vos données Hi‑C ? Exemple d'utilisation de HiC-pro”

  1. Si je peux me per­mettre quelques conseils de style python pour la fonc­tion

    loa­dabs­da­ta­file

    :

    Les com­men­taires expli­quant les argu­ments de la fonc­tion auraient tout à fait leur place dans une doc­tring

    Il est conseillé de gérer l'ouverture et la fer­me­ture des fichiers avec le "context mana­ger" with

    Il est conseillé d'aérer le code, par exemple en met­tant des espaces de part et d'autre des signes "=", et de choi­sir des noms de variables faci­li­tant la com­pré­hen­sion du code.

    On peut faire une boucle

    for

    sur les lignes d'un fichier ouvert en lec­ture "direc­te­ment", sans uti­li­ser

    read­line

    et

    while

    :

    for ligne in fichier :
    chrom, start, end, indice = ligne.split()

    Je n'avais jamais vu

    _​_​contains_​_​

    . Typi­que­ment, on teste si une clé existe dans un dic­tion­naire avec

    in

    :

    if chrom in d :
    d[chrom][1] = float(indice) # Pour­quoi pas int(indice) ?
    else :
    d[chrom] = [float(indice), float(indice)]

    Voi­là, j'espère que le code ne sera pas for­ma­té trop n'importe com­ment.

    Mer­ci pour l'article, en tout cas.

    1. Avatar de Mathurin
      Mathurin

      Mer­ci Bli ! Effec­ti­ve­ment le code pro­po­sé est loin d'être par­fait, je vais pro­fi­ter de tes conseils pour mettre à jour le code :).

Laisser un commentaire