À la plateforme, je fais régulièrement des analyses de données de séquençage de nouvelle génération (Next Generation Sequencing ou NGS). L’une des questions qui revient fréquemment chez mes clients est : comment ouvrir les fichiers de séquences générés? Considérant l’énorme taille de ces fichiers (souvent plusieurs millions de lignes) et, par conséquent, l’espace qu’ils requièrent en mémoire, ils ne devraient pas être ouverts d’une quelconque façon, ils devraient plutôt être processés.

La plupart des programmes conçus pour traiter les données de NGS procèdent séquentiellement (stream), c’est-à-dire qu’ils chargent seulement la quantité de données requise du disque en mémoire, la traitent, produisent en sortie un résultat et libèrent la mémoire utilisée.

Historiquement, les ordinateurs n’ont jamais eu autant de mémoire vive (RAM). Ainsi, dans les années 70 (!), une multitude d’outils ont été développés sous UNIX pour effectuer ce genre d’opérations directement dans un terminal. Ces outils sont maintenant couramment utilisés et sont disponibles par défaut sous Linux et MacOS.

Je présenterai ici quelques-uns des outils les plus utiles pour manipuler et traiter de larges fichiers. Les outils présentés sont de plus des incontournables pour construire de petits pipelines pour l’analyse de données de NGS.

sed

sed est apparu en 1974 comme un éditeur de flux (abréviation de Stream EDitor, « éditeur de flux ») pour UNIX capable d’effectuer toutes sortes de manipulations sur des fichiers textes. Son utilisation la plus fréquente est probablement comme outil de remplacement de texte même si plusieurs autres fonctionnalités sont disponibles. Par exemple, supposons que vous téléchargez une annotation de gènes de Ensembl et que vous réalisez après coup que les chromosomes sont numérotés de 1 à MT alors que dans votre génome de référence, ils étaient numérotés de chr1 à chrM. Deux appels consécutifs à sed vous permettront de convertir le format de l’un vers l’autre: 

sed 's/^/chr/g' Homo_sapiens.GRCh37.75.gtf | sed 's/chrMT/chrM/g' > Homo_sapiens.GRCh37.75.modified.gtf

sed peut aussi être utile pour extraire des parties d’un fichier. Pour un simple fichier delimité par des tabulations (tab-delimited), un champ donné, disons le champ x, peut être extrait facilement avec la commande cut -f x . Par contre, pour un cas moins bien structuré, par exemple un fichier GTF où les attributs sont conservés sous forme de liste sans ordre particulier, vous aurez besoin de la puissance de sed et des expressions régulières pour extraire le champ désiré. Par exemple, pour extraire le nom des gènes dans un fichier GTF, en sachant que le huitième champ contient une chaîne de caractères similaire à gene_name « GAPDH », vous pouvez utiliser la commande sed suivante pour récupérer cette portion de chacune des lignes : 

sed 's/.*gene_name "\([^"]*\)".*;/\1/g' Homo_sapiens.GRCh37.75.gtf

Vous pouvez ensuite diriger (pipe) le résultat vers deux autres commandes très utiles, sort et uniq, pour obtenir une liste non redondante de tous les gènes.

sed 's/.*gene_name "\([^"]*\)".*;/\1/g' Homo_sapiens.GRCh37.75.gtf | sort | uniq > gene_list.txt

awk
awk a été développé quelques années plus tard, en 1977, avec un but similaire en tête, c’est-à-dire faciliter les modifications de fichiers texte une ligne à la fois. awk offre par contre un plus grand nombre de fonctionnalités. En particulier, il permet l’utilisation de conditions et de variables. Prenons un exemple. Supposons que vous voulez calculer la couverture moyenne de régions spécifiques dans une expérience de capture d’exome. Vous pouvez utiliser les excellents outils de bedtools pour calculer, dans un premier temps, la profondeur à chacune des positions des éléments d’un fichier d’annotation et diriger ce résultat vers awk pour obtenir la valeur moyenne de couverture. Dans ce cas-ci, le dernier champ de chaque ligne ($NF) est additionné pour faire une somme, mais seule la moyenne est retournée en sortie.

bedtools coverage -d  -a targets.bed -b sample.bam -split \
| awk '($1 == "chr1") { sum += $NF } END { print sum/NR }'

Redirection de flux (Stream redirection)

Dans l’exemple précédent, nous avons sauvé de l’espace disque en n’écrivant pas la sortie de bedtools sur le disque. Toutefois, si nous voulons maintenant calculer le nombre de bases ayant une couverture supérieure ou égale à 20X, il faudrait attendre quelques minutes pendant que nous recalculons la couverture ou bien il faudrait sauver le résultats intermédiaire sur le disque (typiquement quelques GB). Nous pouvons optimiser ceci en utilisant la commande tee.

Cet outil est utilisé pour séparer la sortie standard (standard output) d’un programme pour qu’elle puisse être à la fois affichée et sauvegardée dans un fichier. En utilisant quelques trucs additionnels, nous pouvons diriger la sortie vers deux processus à la fois : 

bedtools coverage -d -a targets.bed -b sample.bam -split \
 | tee >(awk '{ sum += $NF } END { print sum/NR }' > coverage_mean.txt) \
 | awk '$NF>=20' | wc -l > coverage_20x.txt

Un dernier exemple concernant les données à passer en entrée (input) à des programmes, plus particulièrement quand vous voulez utiliser la sortie d’un programme comme entrée d’un autre programme qui n’a ps été conçu pour accepter stdin comme entrée ou si vous voulez diriger deux stdout en même temps vers le même programme. Reprenons l’exemple utilisant bedtools encore une fois. Nous pouvons calculer la couverture sur le chromosome 1 en utilisant l’unique commande suivante : 

bedtools coverage -d  -a <(grep ^chr1 targets.bed) -b <(samtools view sample.bam chr1) -split

Ici, nous avons remplacé les paramètres de fichiers (-a et -b) par des appels à d'autres commandes. En arrière-plan, des descripteurs de fichiers seront créés comme tampons pour diriger les données entre les programmes, mais encore une fois, rien ne sera écrit sur disque.

Évidemment, nous n'avons survolé ici que la surface de toutes les fonctionnalités offertes par les commandes UNIX dans le terminal. Plusieurs autres outils existent et méritent d'être mentionnés notamment grep, cat, join, cut, wc, etc. Je vous encourage à aller y jeter un coup d'oeil pour améliorer vos compétences en matière de scripts (shell scripting)!