Langage : shell
Commandes présentées : grep, split (succintement)
Niveau : débutant
Présentation de la commande grep
La commande grep est disponible nativement sur la plupart des systèmes d'exploitation GNU/Linux. La plupart des utilisateurs utilisent cette commande pour rechercher un mot ou un groupe de mots, que nous appellerons motif (pattern en anglais), dans un fichier texte. Cependant cette commande ne se limite pas à du simple cas par cas.
Grep recherche le motif en parcourant tout le fichier texte du début jusqu'à la fin. Ainsi, autant pour un fichier de quelques lignes, le résultat sera quasi immédiat, autant pour un fichier de plusieurs milliers de lignes le résultat peut être plus ou moins long.
Dans ce billet je vous présenterai les différentes façons dont je me sers régulièrement de grep, que ce soit de la simple recherche d'un mot clé à la recherche, plus ou moins complexe, de plusieurs motifs.
Pour les exemples qui suivront, je me base sur la liste des informations sur les gènes fournies par Entrez Gene du NCBI. Si vous êtes sous GNU/Linux et que vous souhaitez reproduire les exemples présentés, vous pouvez saisir les commandes suivantes dans un terminal :
1 2 3 4 5 6 |
# wget : télécharge le fichier depuis une source externe<br> wget ftp ://ftp.ncbi.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz<br> # gunzip : décompresse le fichier au format gzip<br> gunzip Homo_sapiens.gene_info.gz<br> # mv : renomme le fichier pour une extension plus explicite<br> mv Homo_sapiens.gene_info Homo_sapiens.tsv |
Méthode basique
La commande grep vous apporte un grand soutien dans la recherche de motif(s), vous l'utilisez déjà certainement de façon très basique comme dans cet exemple :
1 2 3 |
grep HLA-A Homo_sapiens.tsv<br> 9606 3105 HLA-A DAQB-90C11.16-002 HLAA HGNC :4931|MIM :142800|Ensembl :ENSG00000206503|HPRD :00826|Vega :OTTHUMG00000130501 6 6p21.3 major histocompatibility complex, class I, A protein-coding HLA-A major histocompatibility complex, class I, A O HLA class I histocompatibility antigen, A-1 alpha chain|MHC class I antigen HLA-A heavy chain|antigen presenting molecule|leukocyte antigen class I-A 20121110<br> 9606 100873924 HLA-AS1 - - HGNC :42509 1 - HLA antisense RNA 1 miscRNA HLA-AS1 HLA antisense RNA 1O- 20121021 |
Premier constat : si le motif est compris dans un motif plus grand alors que vous ne souhaitez récupérer que les lignes qui correspondent exactement à ce motif, alors vous récupérerez toutes les lignes comprenant ce motif, y compris celles ne vous intéressant pas. Le bruit alors généré peut s'avérer problématique en fonction du nombre d’occurrences retournées !
Dans l'exemple présenté ci-dessus, vous pouvez constater que grep retourne le gène HLA‑A et le gène HLA-AS1.
Rechercher un motif exact
Pour trouver le motif exact, vous pouvez jouer avec les options et plus particulièrement l'option -w (ou - ‑word-regexp). Exemple :
1 2 |
grep -w HLA-A Homo_sapiens.tsv<br> 9606 3105 HLA-A DAQB-90C11.16-002 HLAA HGNC :4931|MIM :142800|Ensembl :ENSG00000206503|HPRD :00826|Vega :OTTHUMG00000130501 6 6p21.3 major histocompatibility complex, class I, A protein-coding HLA-A major histocompatibility complex, class I, A O HLA class I histocompatibility antigen, A-1 alpha chain|MHC class I antigen HLA-A heavy chain|antigen presenting molecule|leukocyte antigen class I-A 20121110 |
Rechercher plusieurs motifs
Pour chercher plusieurs motifs, trois solutions s'offrent à vous :
- faire une boucle sur la liste des motifs et faire des grep successifs ;
- jouer avec les expressions régulières, chose que je vous recommande pour quelques motifs ;
- passer en argument un fichier de motif, qui est plus recommandé si vous avec beaucoup de motifs à chercher.
Méthode de la boucle
Cette méthode est assez basique, je m'en suis servie dans mes tout premiers scripts shell, mais elle présente un gros inconvénient : le temps de calcul !
Bien que grep soit une commande rapide pour un motif, plus vous ferez de grep successifs et plus cela mettra de temps. Ainsi, si il ne faut qu'un dixième de seconde pour afficher le résultat pour un motif, faites le calcul pour la recherche de 100 motifs. Notez toutefois que, plus vous aurez de motifs à rechercher, plus il faudra du temps avant d'afficher les résultats.
Pour les plus débutants d'entre vous je vous présente une syntaxe à reproduire pour cette méthode, comme ça vous pourrez vous faire votre propre idée du temps de calcul en comparant avec les deux autres méthodes.
1 2 3 4 5 6 |
motifs=(HLA-A HLA-DRA ADD2 BCR1)<br> for motif in ${motifs[*]}; do grep -w $motif Homo_sapiens.tsv ; done<br> 9606 3105 HLA-A DAQB-90C11.16-002 HLAA HGNC :4931|MIM :142800|Ensembl :ENSG00000206503|HPRD :00826|Vega :OTTHUMG00000130501 6 6p21.3 major histocompatibility complex, class I, A protein-coding HLA-A major histocompatibility complex, class I, A O HLA class I histocompatibility antigen, A-1 alpha chain|MHC class I antigen HLA-A heavy chain|antigen presenting molecule|leukocyte antigen class I-A 20121110<br> 9606 3122 HLA-DRA DASS-397D15.1 HLA-DRA1|MLRW HGNC :4947|MIM :142860|Ensembl :ENSG00000204287|HPRD :00833|Vega :OTTHUMG00000031269 6 6p21.3 major histocompatibility complex, class II, DR alpha protein-coding HLA-DRA major histocompatibility complex, class II, DR alphaO HLA class II histocompatibility antigen, DR alpha chain|MHC cell surface glycoprotein|MHC class II antigen DRA|histocompatibility antigen HLA-DR alpha 20121110<br> 9606 119 ADD2 - ADDB HGNC :244|MIM :102681|Ensembl :ENSG00000075340|HPRD :00037|Vega :OTTHUMG00000129710 2 2p13.3 adducin 2 (beta) protein-coding ADD2 adducin 2 (beta) O beta-adducin|erythrocyte adducin subunit beta 20121110<br> 9606 613 BCR - ALL|BCR1|CML|D22S11|D22S662|PHL HGNC :1014|MIM :151410|Ensembl :ENSG00000186716|HPRD :01044|Vega :OTTHUMG00000150655 22 22q11.23 breakpoint cluster region protein-coding BCR breakpoint cluster region O BCR/FGFR1 chimera protein|FGFR1/BCR chimera protein|breakpoint cluster region protein|renal carcinoma antigen NY-REN-26 20121110 |
Méthode des expressions régulières
Pour chercher plusieurs motifs à l'aide des expressions régulières, vous allez devoir ajouter l'option -E (ou - ‑extended-regexp). Exemple :
1 2 3 |
grep -wE "HLA‑A|BCR1" Homo_sapiens.tsv<br> 9606 613 BCR - ALL|BCR1|CML|D22S11|D22S662|PHL HGNC :1014|MIM :151410|Ensembl :ENSG00000186716|HPRD :01044|Vega :OTTHUMG00000150655 22 22q11.23 breakpoint cluster region protein-coding BCR breakpoint cluster region O BCR/FGFR1 chimera protein|FGFR1/BCR chimera protein|breakpoint cluster region protein|renal carcinoma antigen NY-REN-26 20121110<br> 9606 3105 HLA-A DAQB-90C11.16-002 HLAA HGNC :4931|MIM :142800|Ensembl :ENSG00000206503|HPRD :00826|Vega :OTTHUMG00000130501 6 6p21.3 major histocompatibility complex, class I, A protein-coding HLA-A major histocompatibility complex, class I, A O HLA class I histocompatibility antigen, A-1 alpha chain|MHC class I antigen HLA-A heavy chain|antigen presenting molecule|leukocyte antigen class I-A 20121110 |
Explication(s) de la ligne de commande : le caractère | est un caractère spécial utilisé dans les expressions régulières. Ce caractère indique à la commande grep qu'il doit trouver au moins un des deux motifs. Si on devait le formuler à un être humain, nous pourrions le traduire par : "dans le fichier Homo_sapiens.tsv, cherche moi toutes les lignes qui ont exactement le terme HLA‑A ou le terme BCR1 et indique-les moi".
Méthode du fichier de motifs
Pour passer un fichier de motifs à la commande grep, vous devez passer par l'option -f (ou - ‑file=FILE).
Le fichier de motifs doit contenir un motif par ligne, ici, un extrait des gènes intervenant dans la glycolyse/glycogénogenèse chez l'homme (source KEGG) :
1 2 3 4 5 6 7 8 9 |
HK3<br> HK1<br> HK2<br> HKDC1<br> [..]<br> ADPGK<br> BPGM<br> PCK1<br> PCK2 |
Une fois votre fichier de motifs créé, vous pouvez le fournir à la commande grep qui gérera alors l'ordre dans lequel les motifs seront recherchés et trouvés :
1 |
grep -w -f "genes.txt" Homo_sapiens.tsv |
Notez toutefois que plus votre fichier de motifs contient des motifs, plus la commande grep mettra du temps avant de vous retourner le résultat. Pour cela je vous conseille de vous orienter du côté de la commande split afin de découper le fichier de motifs en plusieurs fichiers comme dans l'exemple suivant :
1 2 3 4 5 6 |
split -l 30 genes.txt genes.split.<br> for gene in $( ls genes.split.* )<br> do<br> grep -f $gene -w fichier.txt<br> rm $gene<br> done |
Cette méthode, dont j'ai trouvé la source sur ce blog en anglais, présente un bien meilleur avantage que celui de faire une boucle sur la liste des motifs avant de faire des grep successifs. Avec un fichier de motifs votre recherche sera plus rapide que si vous deviez faire 30 grep les uns à la suite des autres, et votre machine n'en sera que plus heureuse !
Il est vrai que sur l'exemple que je vous ai fourni ce n'est pas très parlant. Sur 65 gènes, grep ‑f ne met pas plus de temps qu'avec plusieurs fichiers. Cependant j'ai pu constater une importante différence, en terme de temps, sur une liste de 9 312 motifs, et ça, c'est sans contexte un must to know pour des motifs très nombreux.
Savoir sur quelle ligne apparait le motif
L'un des avantages que grep nous apporte, c'est la possibilité de connaître la ligne du fichier sur laquelle le motif a été trouvé. Pour cela vous devrez utiliser l'option -n comme dans l'exemple ci-dessous :
1 |
grep -w BCR1 Homo_sapiens.tsv -n |
Sur mon système, une simple Crunchbang installée sous VirtualBox, voici le résultat que j'obtiens :
1 2 |
grep -w BCR1 Homo_sapiens.tsv -n<br> 506 :9606 613 BCR - ALL|BCR1|CML|D22S11|D22S662|PHL HGNC :1014|MIM :151410|Ensembl :ENSG00000186716|HPRD :01044|Vega :OTTHUMG00000150655 22 22q11.23 breakpoint cluster region protein-coding BCR breakpoint cluster region O BCR/FGFR1 chimera protein|FGFR1/BCR chimera protein|breakpoint cluster region protein|renal carcinoma antigen NY-REN-26 20121110 |
Le gène BCR1 a été trouvé sur la ligne 506. Et pour les amoureux des couleurs qui sont parmi vous, essayez l'option - ‑color si votre motif n'apparaît pas explicitement 😉 !
Afficher 2 lignes avant et 2 lignes après le motif trouvé
Grep est également capable de vous permettre d'afficher un nombre de lignes avant ou après le motif une fois celui-ci trouvé, en voici un exemple :
1 2 3 4 5 6 |
grep -w BCR1 -B 2 -A 2<br> 9606 610 HCN2 - BCNG-2|BCNG2|HAC-1 HGNC :4846|MIM :602781|Ensembl :ENSG00000099822|HPRD :07213|Vega :OTTHUMG00000180590 19 19p13.3 hyperpolarization activated cyclic nucleotide-gated potassium channel 2 protein-coding HCN2 hyperpolarization activated cyclic nucleotide-gated potassium channel 2 O brain cyclic nucleotide gated channel 2|brain cyclic nucleotide-gated channel 2|potassium/sodium hyperpolarization-activated cyclic nucleotide-gated channel 2 20121110<br> 9606 611 OPN1SW - BCP|BOP|CBT HGNC :1012|MIM :613522|Ensembl :ENSG00000128617|HPRD :01836|Vega :OTTHUMG00000158311 7 7q32.1 opsin 1 (cone pigments), short-wave-sensitive protein-coding OPN1SW opsin 1 (cone pigments), short-wave-sensitive O blue cone photoreceptor pigment|blue-sensitive opsin|short-wave-sensitive opsin 1 20121110<br> 9606 613 BCR - ALL|BCR1|CML|D22S11|D22S662|PHL HGNC :1014|MIM :151410|Ensembl :ENSG00000186716|HPRD :01044|Vega :OTTHUMG00000150655 22 22q11.23 breakpoint cluster region protein-coding BCR breakpoint cluster region O BCR/FGFR1 chimera protein|FGFR1/BCR chimera protein|breakpoint cluster region protein|renal carcinoma antigen NY-REN-26 20121110<br> 9606 616 BCRP4 - BCR-4|BCR4|BCRL4 HGNC :1017|MIM :113660 22 22q11.22 breakpoint cluster region pseudogene 4 pseudo BCRP4 breakpoint cluster region pseudogene 4 O - 20121110<br> 9606 617 BCS1L - BCS|BCS1|BJS|FLNMS|GRACILE|Hs.6719|PTD|h-BCS HGNC :1020|MIM :603647|Ensembl :ENSG00000074582|HPRD :04708|Vega :OTTHUMG00000133114 2 2q33 BC1 (ubiquinol-cytochrome c reductase) synthesis-like protein-coding BCS1L BC1 (ubiquinol-cytochrome c reductase) synthesis-like O BCS1-like protein|h-BCS1|mitochondrial chaperone BCS1|mitochondrial complex III assembly 20121110 |
Le mot de la fin
La commande grep, bien qu'étant une commande très informatique et non destinée au domaine de la bioinformatique, est un outil majeur dont nous nous servons quotidiennement. Elle s'avère être un outil puissant pour peu que l'on sache l'utiliser. De plus, si vous avez des notions sur les pipelines, vous pouvez constater qu'elle peut être utilisée à d'autres fins que la simple recherche de motifs, en combinaison avec d'autres commandes et/ou d'autres programmes.
Merci à Guillaume Collet et Yoann M. pour leur relecture, et nos discussions enrichissantes sur cette commande.
Laisser un commentaire