29.《Bioinformatics Data Skills》之

很多命令行的工具速度非常快速，包括专门用于匹配的工具grep(图1)，追求速度的时候可以使用。grep非常强大，这里介绍一下它的常见用法。

图1 grep与其它工具速度对比

grep基本输入为需要匹配的字符串与文本文件（如下，这里加引号只是为了安全，可以选择不加）：

$ grep "Olfr1413" Mus_musculus.GRCm38.75_chr1_genes.txt
ENSMUSG00000058904      Olfr1413

grep是通过部分匹配的方式，任何包含字符串的行都会被返回：

$ grep "Olfr" Mus_musculus.GRCm38.75_chr1_genes.txt | head -n3
ENSMUSG00000067064      Olfr1416
ENSMUSG00000057464      Olfr1415
ENSMUSG00000042849      Olfr1414

如果我们想要匹配Olfr开头但是不包括Olfr1413的基因，可以使用grep的组合(如下，-v取否)：

$ grep "Olfr" Mus_musculus.GRCm38.75_chr1_genes.txt | grep -v "Olfr1413" | head -n3
ENSMUSG00000067064      Olfr1416
ENSMUSG00000057464      Olfr1415
ENSMUSG00000042849      Olfr1414

不过以上的做法是存在问题的，因为结果不仅会排除掉基因Olfr1413，还会排除掉基因Olfr1314a等（如果有的话）。为了避免这种情况，grep有一个参数-w，代表完全匹配（即匹配的字符串前后为空格）。通过下面的例子比较结果的区别：

$ cat example.txt
bio
bioinfo
bioinformatics
computational biology

# 未加-w返回部分匹配的结果
$ grep "bioinfo" example.txt
bioinfo
bioinformatics

# 加-w后只返回完全匹配的结果
$ grep -w "bioinfo" example.txt
bioinfo

grep的返回结果只包括匹配的行（例如下）:

$ grep "AGATCGG" contam.fastq | head -n 3
TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG
GCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAAC

但有时我们想要查看上下文，这个时候可以使用以下参数：

1. -A{n}代表展示匹配行后（After）n行的信息，例如查看后1行的信息（结果间以"--"分隔）：

   $ grep -A1 "AGATCGG" contam.fastq | head -n 6
   TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
   +
   --
   CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG
   +
   --
   

2. -B{n}则展示前(Back)n行的信息，例如：

   $ grep -B1 "AGATCGG" contam.fastq | head -n 6
   @DJB775P1:248:D0MDGACXX:7:1202:12362:49613
   TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
   --
   @DJB775P1:248:D0MDGACXX:7:1202:12782:49716
   CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG
   --
   

3. -C{n}代表同时展示前后n行的信息，例如：

   $ grep -C1 "AGATCGG" contam.fastq | head -n 6
   @DJB775P1:248:D0MDGACXX:7:1202:12362:49613
   TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
   +
   --
   @DJB775P1:248:D0MDGACXX:7:1202:12782:49716
   CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG
   

grep是支持正则表达式的，比如说我们想要匹配基因Olfr1413或Olfr1411，我们可以采用如下方式：

$ grep "Olfr141[13]" Mus_musculus.GRCm38.75_chr1_genes.txt
ENSMUSG00000058904      Olfr1413
ENSMUSG00000062497      Olfr1411

不过这样只能使用基本的正则表达式，可以指定-E参数扩展到perl级别。例如我们想要匹配基因Olfr1413或者Olfr1408：

$ grep -E "(Olfr1413|Olfr1408)" Mus_musculus.GRCm38.75_chr1_genes.txt
ENSMUSG00000058904      Olfr1413
ENSMUSG00000062527      Olfr1408

当我们只关注匹配的结果的数目的时候，可以采用-c参数，比如说我们想要看以Olfr开头的基因的数目：

$ grep -c $'\tOlfr' Mus_musculus.GRCm38.75_chr1_genes.txt
27

grep一旦匹配到结果后，会返回整行的信息。有时我们只关心匹配到的信息，可以采用-o参数只保留匹配信息：

$ grep -o -E 'gene_id "(\w)+"' Mus_musculus.GRCm38.75_chr1.gtf | head -n3
gene_id "ENSMUSG00000090025"
gene_id "ENSMUSG00000090025"
gene_id "ENSMUSG00000090025"

我们通过以下方式将所有的gene id都保存下来：

grep -E -o 'gene_id "(\w)+"' Mus_musculus.GRCm38.75_chr1.gtf |\
cut -f2 -d" " | \
sed 's/"//g' | \
sort |\
uniq > mm_gene_id.txt

$ wc -l mm_gene_id.txt
2027 mm_gene_id.txt

以上就是grep常见的用法介绍，其余用法可以通过man grep了解。