生物信息 awk 简明教程和基本用法

全文6,271字，阅读16分钟。配图来源：Julia Evens

----/ START /----

awk 是处理文本文件的一个应用程序，几乎所有的Linux以及MacOS都自带这个程序。

在这篇文章中，我想给大家介绍如何用这个程序来解决一些基本的生物信息数据处理和文本处理的问题，特别适合对此不熟悉的同学和读者朋友。

简述

我们日常项目中很多的数据分析和处理工作并非都需要编写复杂的程序才能完成，很多小修小改、查找、替换、简单的数据计算等工作，其实可以用一些Linux/MacOS中自带的命令行工具来完成。awk 就是这一类工具中的一个，它依次处理文件中的每一行，并读取里面的每一个字段，对于我们在生信中很多每行格式都相同的文本文件来说，awk 可能是最方便的一个工具，不但可以省去很多不必要的脚本和程序，还可以通过对它的灵活应用，快速理解手中的数据。

其实，把 awk 说成是一个程序工具并不十分准确。实际上，它还是一种解释型的编程语言（类似于Perl），即写即用，响应快，错了重改也方便，也有人习惯称这一类编程语言为脚本语言。不过在这里我只介绍它的命令行用法，对于很多生物信息的数据分析场景，应该是足够的，与之类似的还有 sed。

事先说明一下，awk 毕竟是命令行工具，所以我在这篇文章中所用到的例子都只能在Linux或者MacOS命令行中才能执行。如果你用的是windows电脑，那么需要安装一个命令行工具才可以（比如git bash——https://segmentfault.com/a/1190000006683607，而不能是win自带的命令行）。

基本用法 接下来就是正文了，awk 其实十分简单，它在命令行中的基本用法就是下面这个形式：

# 标准用法的形式
$ awk 处理动作 文件名
# 例子
$ awk '{print $0}' demo.vcf

在这个例子中，demo.vcf 是 awk 要处理的文本文件——注意我这里反复强调必须是文本文件，而不是BAM或者.gz这一类非文本文件，如果想用 awk 处理这类文件，那么需要先转换为文本文件才行，假如文件不大，那么可以不做单独转换，直接用管道操作来完成即可。回到刚刚的例子，demo.vcf 前面的单引号内有一个大括号（注意，这个单引号是必须的，而在包含判断、输出等复杂语句的时候大括号也是必须的），里面是对文件中每一行内容的处理动作，比如这里是：print $0，其中 print 是打印命令，而 $0 代表当前完整的一行，所以上面这个命令的执行结果就是把 demo.vcf 每一行都原样打印出来。

由于文章这里不便打开 demo.vcf 文件，我们就用标准输入来演示这个例子吧。

$ echo "this is a variant in vcf file" | awk '{print $0}'
this is a variant in vcf file

上面代码中，echo 也是 linux 一个命令行工具，作用是输出标准输入的文本内容，这里不展开，print $0 就是把标准输入的 "this is a variant in vcf file" 这一句话，重新输出到屏幕。

大家应该也注意到上面的命令里有一个 “|” 竖线，这个就是 Linux/Unix/MacOS 的管道操作符，一个非常有用的符号。它可以把前一个命令的结果作为标准输入传输到后一个命令中去处理，而且还可以多重串联下去，就像成语接龙一样，前一个管道处理完再传给下一个管道去处理，然后再下一个，如果你愿意的话，甚至可以一直接下去，这样做的好处是减少系统 IO，加快处理速度。我前面说到 awk 只能处理文本文件，那当我们的文件不是文本格式时，比如是 gz 压缩文件或者BAM文件的时候，要用 awk 处理的话，就需要先做转换然后通过管道把数据传过给 awk 来分析，比如：

$ samtools view demo.bam | awk '{print $0}'

这里就是先通过 samtools view 将 demo.bam 转为可读的文本，然后用管道("|")把数据传到后面的 awk 中去处理——我们这里为例子的方便只是做原样输出。不过，通过这种形式进行数据分析的时候，应该注意的地方是，被处理的 demo.bam 文件不能太大，否则，管道前一个命令（samtools view）转换出来的文本信息会一直累积到计算机内存中，最后很可能把机器内存撑爆。

默认情况下，awk 将根据空格和制表符(tab)，把每一行自动切分成若干个字段，并在系统里依次用 $1,$2,$3,... 代表第一个字段、第二个字段、第三个字段等等。

$ echo "this is a variant in vcf file" | awk '{print $4}'
variant

上述代码中，$4 代表了 "this is a variant in vcf file" 这一句话中的第四个字段 "variant"。如果把这一段话换为一份文件，那么这个命令就会把那份文件中各行的第四列都打印输出出来。

除此之外，对于某些不是以空格和tab作为分隔符存储的文件，或者在文件中的某一列的信息中是以其它分隔符串接起来的，比如 VCF 的 INFO 那一列，它是 VCF 的第八列，该列中的信息往往比较丰富，并且各个字段之间是通过逗号(,)连接在一起的，如下：

CMDB_AF=0.030044,CMDB_AC=420,CMDB_AN=13442
CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
CMDB_AF=0.053564,CMDB_AC=534,CMDB_AN=9525

此时，如果我们想获得该 INFO 中的某一个字段信息，那么这个时候除了要提取出这一列之外，还需要通过自定义输入分隔符才能将其进行切割。自定义输入分隔符，在 awk 中用的是 -F 参数，例子：

$ awk '{if($1!~"^#"){print $8}}' demo.vcf | awk -F ',' '{print $2}'
CMDB_AC=420
CMDB_AC=441
CMDB_AC=441
CMDB_AC=842
CMDB_AC=842
CMDB_AC=534

这个命令完成了 demo.vcf 中 INFO 这一列信息中第2个字段信息的提取。其中 通过 -F 参数重新设置了输入分隔符为逗号，从而完成了对INFO的切分，然后再提取出字段。该操作命令中前半部分的语句 "if($1!~"^#"){print $8}" 是把VCF 的header信息过滤掉，由于 VCF 的 Header 中每一行都是以 # 开头的，所以 $1!~"^#" 就可以忽略掉这些 # 开头的行。

BEGIN 语句

另外在上面的例子中，除了使用 -F 参数之外，还有另一个方法也可以完成这个操作，就是通过 BEGIN 语句，在执行实际命令之前初始化输入分隔符：

$ awk '{if($1!~"^#"){print $8}}' demo.vcf | awk 'BEGIN{FS=","}{print $2}'

上面这个命令也能够实现同样的功能，只不过，在 BEGIN 中就不是通过 -F 参数了，代替之的是一个内置变量 FS。同时，如果需要的话，我们还可以在其中设置多重分隔符，如 FS="[:,]"(或者 -F '[:,]')，代表同时用冒号和逗号作为输入分隔符切分数据，这种方式在比较复杂的文本环境中应用起来会更加方便。

此外，既然可以设置输入分隔符，自然也可以定义输出分隔符。我这里还是用 BEGIN 定义作为例子：

$ awk '{if($1!~"^#"){print $8}}' demo.vcf | awk 'BEGIN{FS=","; OFS="###"}{print $1,$2,$3}'
CMDB_AF=0.030044###CMDB_AC=420###CMDB_AN=13442
CMDB_AF=0.031047###CMDB_AC=441###CMDB_AN=13553
CMDB_AF=0.031047###CMDB_AC=441###CMDB_AN=13553
CMDB_AF=0.050419###CMDB_AC=842###CMDB_AN=16135
CMDB_AF=0.050419###CMDB_AC=842###CMDB_AN=16135
CMDB_AF=0.053564###CMDB_AC=534###CMDB_AN=9525

awk 默认的输出分隔符是空格，这个例子在 BEGIN 语句中则通过 OFS 参数将输出分隔符修改为 "###"，当然，最后想用什么输出分隔符，完全取决于我们的实际需要。

有BEGIN就有END

与 BEGIN 语句对应的是 END 语句。awk 在默认情况，是每处理完一行数据，就可以输出一次。但是有时候，我们不希望那么快地把信息输出，特别是当我们进行累加统计的时候，比如，计算基因组的平均覆盖深度，这个时候，我们希望能够在读取完整份文件之后，再统一输出结果。这个时候 END 就派上用场了，顾名思义，只有直到文件结束了，它才把最后结果输出出来。

$ awk '{total_depth += ($3-$2)*$4; total_length += $(3-$2);}END{print total_depth/total_length}' seq_depth.bed
53.4

上面的代码里面 seq_depth.bed 是通过 bedtools2 生成的一份基因组覆盖深度文件，为 bed 格式，第一列是染色体ID，第二列是起始位置，第三列是终止位置，第四列是该区域各个位点的覆盖深度，其中每一个bed区域里各个位点的深度都是一样的，所以只留下一个值，这也是为什么我在上面累加深度的时候需要用 ($3-$2)*$4 的原因。在整个命令中，直到最后读完整份 seq_depth.bed 才print 出最终的平均深度，比如这里的 53.4。

内置变量

其实，除了上述通过 $+数字 的形式表示某个字段之外，awk 本身还有一些默认变量。其中包括，变量 NF 表示当前行按照输入分隔符切分之后一共有多少列（或者说多少字段），所以 $NF就表示最后一个字段，在一些列数非常多的文件中 NF 是很有用的，我们不用数数 数到眼花，也能立刻获得最后一个字段，或者立刻知道每一行都有多少字段。如：

$ echo "this is a variant in vcf file" | awk '{print $NF}'
file

$(NF-1)则代表倒数第二个字段，$(NF-2)表示倒数第三，依此类推。

有表示列数，自然也就有表示行数的。awk 中的变量 NR 就是表示当前所处理的是第几行。

$ awk '{if($1!~"^#"){print $8}}' demo.vcf | awk -F ',' '{print NR")", $2}'
1) CMDB_AC=420
2) CMDB_AC=441
3) CMDB_AC=441
4) CMDB_AC=842
5) CMDB_AC=842
6) CMDB_AC=534

在这个例子中唯一需要注意的是，print 输出的字段中，如果各个字段之间没通过逗号隔开，那么输出时，中间也不会加入任何分隔符，比如这里 NR 后面直接跟了 ")"，输出的时候 ")" 就紧贴着行数出来。

awk 内置的变量还有这些，其实有不少我们在上面已经用过了，这里再做汇总：

FILENAME：当前文件名
FS：字段分隔符，默认是空格和制表符
RS：行分隔符，用于分割每一行，默认是换行符
OFS：输出字段的分隔符，用于打印时分隔字段，默认为空格
ORS：输出记录的分隔符，用于打印时分隔记录，默认为换行符
OFMT：数字输出的格式，默认为％.6g

内置函数

awk 除了有好用的内置变量之外，也提供了不少好用的内置函数。这些函数可以让我们很方便地对原始数据进行一些基本的处理。比如，tolower() 用于将字符转换为小写。

$ awk '{if($1!~"^#"){print $8}}' demo.vcf | awk -F ',' 'tolower($2)}'
cmdb_ac=420
cmdb_ac=441
cmdb_ac=441
cmdb_ac=842
cmdb_ac=842
cmdb_ac=534

上述代码中，tolower把对应字段都转换为了小写。其他的常用函数还有如下这些：

tolower()：字符转为小写。
length()：返回字符串长度。
substr()：返回子字符串。
sin()：正弦。
cos()：余弦。
sqrt()：平方根。
rand()：随机数。

但是真正完整的内置函数列表，需要查看awk手册(https://www.gnu.org/software/gawk/manual/html_node/Built_002din.html#Built_002din)

条件判断

awk还可以自定条件判断语句，只把符合条件预期的结果输出。命令模式如：

$ awk '条件 动作' 文件名

需要注意的是，条件判断要写在动作之前。请看下面一个例子：

$ awk '$6 > 40' demo.vcf
这里只把 demo.vcf 中第六列大于40（也就是质量值>40）的行输出出来。

我们也可以写一个正则表达式，把符合匹配条件的行输出，比如上述例子也出现过，把VCF的Header过滤掉：

$ awk '$1!~/^#/' demmo.vcf

条件判断是很自由的，我们可以依据自己的需要任意设计条件，包括大于、小于、等于、匹配、与或非、异或等等逻辑判断条件都可以设置。

if 语句

除了上面提到的条件判断之外，awk也有 if 语句，可以用来编写更加灵活复杂的条件，上述例子也已经有所涉及。

$ awk '{if($1!~"^#" && $6>40){print $8}}' demo.vcf
CMDB_AF=0.030044,CMDB_AC=420,CMDB_AN=13442
CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
CMDB_AF=0.053564,CMDB_AC=534,CMDB_AN=9525

这个代码的意思就是，只把VCF文件中，质量值大于40的变异的INFO信息 print 出来。

有 if 结构自然也会有与它配对的 else 部分。

$ awk '{if(/^#/){print $0}else{if($6>40){print $0}}}' demo.vcf

(内容比较多，附图)

小结

有关 awk 的基本用法，这篇文章就介绍到这里吧，下一篇是 awk 的进阶（进阶篇已经优先在我的知识星球中给所有星友分享了）。虽然还不是十分全面，但我觉得能够掌握好上面的使用方法，并灵活应用，其实已经可以用一行命令处理很多基本的分析需求了，不必为了一个小功能费劲去写一个程序。更重要的是，像 awk 这一类的处理程序，可以很方便且快速地帮我们了解一个刚接触到的数据，对于加深对数据的理解都是很有好处的，并不用总是要等到你学会写程序之后才行。更多的 awk 教程推荐这本书《Effective awk programming》，貌似只有英文版。

下面这个是《sed 和 awk》中文版的：

参考链接

http://www.ruanyifeng.com/blog/2018/11/awk.html http://www.runoob.com/linux/linux-comm-awk.html

----/ END /----

※ ※ ※

你还可以读

• 如何理解GWAS中Manhattan plot和QQ plot所传递的信息
• 如何有效使用CMDB基因频率数据库
• 我应该如何正确设置GATK VQSR的模型训练参数？
• 一篇文章说清楚基因组结构性变异检测的方法
• 一篇文章说清楚什么是“插入片段”？