Chercher des motifs dans un fichier

Lan­gage : shell
Com­mandes pré­sen­tées : grep, split (suc­cin­te­ment)
Niveau : débu­tant

Présentation de la commande grep

La com­mande grep est dis­po­nible nati­ve­ment sur la plu­part des sys­tèmes d'exploitation GNU/​Linux. La plu­part des uti­li­sa­teurs uti­lisent cette com­mande pour recher­cher un mot ou un groupe de mots, que nous appel­le­rons motif (pat­tern en anglais), dans un fichier texte. Cepen­dant cette com­mande ne se limite pas à du simple cas par cas.
Grep recherche le motif en par­cou­rant tout le fichier texte du début jusqu'à la fin. Ain­si, autant pour un fichier de quelques lignes, le résul­tat sera qua­si immé­diat, autant pour un fichier de plu­sieurs mil­liers de lignes le résul­tat peut être plus ou moins long.
Dans ce billet je vous pré­sen­te­rai les dif­fé­rentes façons dont je me sers régu­liè­re­ment de grep, que ce soit de la simple recherche d'un mot clé à la recherche, plus ou moins com­plexe, de plu­sieurs motifs.

Pour les exemples qui sui­vront, je me base sur la liste des infor­ma­tions sur les gènes four­nies par Entrez Gene du NCBI. Si vous êtes sous GNU/​Linux et que vous sou­hai­tez repro­duire les exemples pré­sen­tés, vous pou­vez sai­sir les com­mandes sui­vantes dans un ter­mi­nal :

# wget : télécharge le fichier depuis une source externe

wget ftp://ftp.ncbi.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz

# gunzip : décompresse le fichier au format gzip

gunzip Homo_sapiens.gene_info.gz

# mv : renomme le fichier pour une extension plus explicite

mv Homo_sapiens.gene_info Homo_sapiens.tsv

Méthode basique

La com­mande grep vous apporte un grand sou­tien dans la recherche de motif(s), vous l'utilisez déjà cer­tai­ne­ment de façon très basique comme dans cet exemple :

grep HLA-A Homo_sapiens.tsv

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

9606    100873924   HLA-AS1 -   -   HGNC:42509  1   -   HLA antisense RNA 1 miscRNA HLA-AS1 HLA antisense RNA 1O-   20121021

Pre­mier constat : si le motif est com­pris dans un motif plus grand alors que vous ne sou­hai­tez récu­pé­rer que les lignes qui cor­res­pondent exac­te­ment à ce motif, alors vous récu­pé­re­rez toutes les lignes com­pre­nant ce motif, y com­pris celles ne vous inté­res­sant pas. Le bruit alors géné­ré peut s'avérer pro­blé­ma­tique en fonc­tion du nombre d’occurrences retour­nées !
Dans l'exemple pré­sen­té ci-des­sus, vous pou­vez consta­ter que grep retourne le gène HLA-A et le gène HLA-AS1.

Rechercher un motif exact

Pour trou­ver le motif exact, vous pou­vez jouer avec les options et plus par­ti­cu­liè­re­ment l'option -w (ou - -word-regexp). Exemple :

grep -w HLA-A Homo_sapiens.tsv

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 cher­cher plu­sieurs motifs, trois solu­tions s'offrent à vous :

• faire une boucle sur la liste des motifs et faire des grep suc­ces­sifs ;
• jouer avec les expres­sions régu­lières, chose que je vous recom­mande pour quelques motifs ;
• pas­ser en argu­ment un fichier de motif, qui est plus recom­man­dé si vous avec beau­coup de motifs à cher­cher.

Méthode de la boucle

Cette méthode est assez basique, je m'en suis ser­vie dans mes tout pre­miers scripts shell, mais elle pré­sente un gros incon­vé­nient : le temps de cal­cul !
Bien que grep soit une com­mande rapide pour un motif, plus vous ferez de grep suc­ces­sifs et plus cela met­tra de temps. Ain­si, si il ne faut qu'un dixième de seconde pour affi­cher le résul­tat pour un motif, faites le cal­cul pour la recherche de 100 motifs. Notez tou­te­fois que, plus vous aurez de motifs à recher­cher, plus il fau­dra du temps avant d'afficher les résul­tats.
Pour les plus débu­tants d'entre vous je vous pré­sente une syn­taxe à repro­duire pour cette méthode, comme ça vous pour­rez vous faire votre propre idée du temps de cal­cul en com­pa­rant avec les deux autres méthodes.

motifs=(HLA-A HLA-DRA ADD2 BCR1)

for motif in ${motifs[*]}; do grep -w $motif Homo_sapiens.tsv; done

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

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

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

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 cher­cher plu­sieurs motifs à l'aide des expres­sions régu­lières, vous allez devoir ajou­ter l'option -E (ou - -exten­ded-regexp). Exemple :

grep -wE "HLA-A|BCR1" Homo_sapiens.tsv

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

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 com­mande : le carac­tère | est un carac­tère spé­cial uti­li­sé dans les expres­sions régu­lières. Ce carac­tère indique à la com­mande grep qu'il doit trou­ver au moins un des deux motifs. Si on devait le for­mu­ler à un être humain, nous pour­rions le tra­duire par : "dans le fichier Homo_sapiens.tsv, cherche moi toutes les lignes qui ont exac­te­ment le terme HLA-A ou le terme BCR1 et indique-les moi".

Méthode du fichier de motifs

Pour pas­ser un fichier de motifs à la com­mande grep, vous devez pas­ser par l'option -f (ou - -file=FILE).
Le fichier de motifs doit conte­nir un motif par ligne, ici, un extrait des gènes inter­ve­nant dans la glycolyse/​glycogénogenèse chez l'homme (source KEGG) :

HK3

HK1

HK2

HKDC1

[..]

ADPGK

BPGM

PCK1

PCK2

Une fois votre fichier de motifs créé, vous pou­vez le four­nir à la com­mande grep qui gére­ra alors l'ordre dans lequel les motifs seront recher­chés et trou­vés :

grep -w -f "genes.txt" Homo_sapiens.tsv

Notez tou­te­fois que plus votre fichier de motifs contient des motifs, plus la com­mande grep met­tra du temps avant de vous retour­ner le résul­tat. Pour cela je vous conseille de vous orien­ter du côté de la com­mande split afin de décou­per le fichier de motifs en plu­sieurs fichiers comme dans l'exemple sui­vant :

split -l 30 genes.txt genes.split.

for gene in $( ls genes.split.* )

do

    grep -f $gene -w fichier.txt

    rm $gene

done

Cette méthode, dont j'ai trou­vé la source sur ce blog en anglais, pré­sente un bien meilleur avan­tage que celui de faire une boucle sur la liste des motifs avant de faire des grep suc­ces­sifs. 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 heu­reuse !
Il est vrai que sur l'exemple que je vous ai four­ni ce n'est pas très par­lant. Sur 65 gènes, grep -f ne met pas plus de temps qu'avec plu­sieurs fichiers. Cepen­dant j'ai pu consta­ter une impor­tante dif­fé­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 nom­breux.

Savoir sur quelle ligne apparait le motif

L'un des avan­tages que grep nous apporte, c'est la pos­si­bi­li­té de connaître la ligne du fichier sur laquelle le motif a été trou­vé. Pour cela vous devrez uti­li­ser l'option -n comme dans l'exemple ci-des­sous :

grep -w BCR1 Homo_sapiens.tsv -n

Sur mon sys­tème, une simple Crunch­bang ins­tal­lée sous Vir­tual­Box, voi­ci le résul­tat que j'obtiens :

grep -w BCR1 Homo_sapiens.tsv -n

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é trou­vé sur la ligne 506. Et pour les amou­reux des cou­leurs qui sont par­mi vous, essayez l'option - -color si votre motif n'apparaît pas expli­ci­te­ment 😉 !

Afficher 2 lignes avant et 2 lignes après le motif trouvé

Grep est éga­le­ment capable de vous per­mettre d'afficher un nombre de lignes avant ou après le motif une fois celui-ci trou­vé, en voi­ci un exemple :

grep -w BCR1 -B 2 -A 2

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

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

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

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

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 com­mande grep, bien qu'étant une com­mande très infor­ma­tique et non des­ti­née au domaine de la bio­in­for­ma­tique, est un outil majeur dont nous nous ser­vons quo­ti­dien­ne­ment. Elle s'avère être un outil puis­sant pour peu que l'on sache l'utiliser. De plus, si vous avez des notions sur les pipe­lines, vous pou­vez consta­ter qu'elle peut être uti­li­sée à d'autres fins que la simple recherche de motifs, en com­bi­nai­son avec d'autres com­mandes et/​ou d'autres pro­grammes.

---

Mer­ci à Guillaume Col­let et Yoann M. pour leur relec­ture, et nos dis­cus­sions enri­chis­santes sur cette com­mande.