生物信息 awk 用法进阶

作者: 黄树嘉 | 来源:发表于2019-04-20 12:24 被阅读41次
    图源:王颜公子

    全文6,829字(含代码),阅读18分钟。

    在掌握了上一篇文章中 awk 基础用法的之后,这一篇文章我们来进一步深入地理解和应用 awk。

    理解AWK的工作原理

    首先,第一个应该加深理解的地方就是 awk 的工作原理(或者说是执行流程)。理解了其工作原理本身,也有助于我们写出更好的 awk 。下面这个图来自 runoob.com 上一篇关于 awk 的文章,它非常清楚明白地描述出了 awk 的工作原理和执行流程,可以说理解 awk 的原理看这一张图几乎就足够了(下图)。

    图源:runoob.com

    总的来说,awk 的执行流程可以分成三个大的部分:

    • 读输入文件之前需要执行的代码段,由 BEGIN 关键字所标识;

    • BODY块,这里是自动循环并处理输入文件的代码段,也是我们处理数据的核心之处,默认情况下,我们编写的 awk 其实都是BODY块;

    • 读取并处理了全部输入文件的内容之后才执行的代码段,由 END 关键字所标识。

    命令的结构如下:

    $ awk  'BEGIN{动作} pattern{动作} END{动作}'
    

    这里的 pattern 属于BODY块,你可以写上一些正则表达式或者条件判断语句,虽然这些语句也可以在 大括号{} 里正式的BODY块中完成,但是写在外面可以使整个命令看起来更加清爽。如:

    $ awk  'BEGIN{OFS="\t";print "#CHROM\tPOS\tINFO"} $1!~/^#/ && $6>40 {print $1,$2,$8}' demo.vcf #CHROM POS INFO chr22 17662679 CMDB_AF=0.030044,CMDB_AC=420,CMDB_AN=13442
    chr22 17662699 CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
    chr22 17662699 CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
    chr22 17662793 CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
    chr22 17662793 CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
    chr22 17663076 CMDB_AF=0.053564,CMDB_AC=534,CMDB_AN=9525
    

    上面的语句就是这样的一个例子,BEGIN 中设定了输出内容的表头和输出分隔符,然后是 pattern,接着是BODY块的主程序。

    所以,awk 的工作原理和执行流程是这样的:

      1. 在所有处理操作之前,先读取 BEGIN 关键字标识起来的代码段,并执行之,给一些预设变量赋值或者输出表头信息;
      1. 然后执行 BODY 块,一行一行往下完成文本的处理;
      1. 在 BODY 执行过程中,对每一行,按照指定的分隔符,把当前整行的内容进行切分,并填充到 awk 内置的数据域中,如 0 标示所有数据域(也就是原来的行内容),1 表示第一个域,$n 表示第 n 个域;
      1. 如果 BODY 前有 pattern 匹配和条件判断语句,那么在依次执行时,只有符合 pattern 条件的才会执行 BODY 中的动作;
      1. 循环读取到整个文件结束之后,就完成了 BODY 块的执行;
      1. 执行 END 代码段,在 END 块中完成最终结果的输出。

    自定义变量

    在看过上一篇文章之后,我想大家一定还多少还记得 awk 的内置变量(比如 NF,FS,OFS等),它们可以帮助我们完成很多的事情。但是内置的变量毕竟是固定的,缺乏灵活性,有些操作它们就不能够胜任了,特别是当我们需要从外部传入参数的时候,它们就通通都不好使了。这个时候我们就需要有一个能够自定义变量的方式,-v 参数在 awk 中就是用于补足这一个需求的,它是这样使用的:

    $ awk -v 变量名字和赋值 '{动作}' 文件名 
    

    来一个实际的例子:

    $ awk -v qual=40 '$1!~/^#/ && $6>qual {print $1,$2,$8}' demo.vcf
    chr22 17662679 CMDB_AF=0.030044,CMDB_AC=420,CMDB_AN=13442
    chr22 17662699 CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
    chr22 17662699 CMDB_AF=0.031047,CMDB_AC=441,CMDB_AN=13553
    chr22 17662793 CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
    chr22 17662793 CMDB_AF=0.050419,CMDB_AC=842,CMDB_AN=16135
    chr22 17663076 CMDB_AF=0.053564,CMDB_AC=534,CMDB_AN=9525
    

    在上面这个例子里,我们通过 -v 参数设置一个自定义变量 qual 并给它赋值为 40, 然后在BODY主程序中 qual 被用于一个条件判断语句,把符合这个条件的 demo.vcf 内容输出出来,非常方便。而且对于自定义变量来说,最大的一个好处是,让 awk 可以和外部进行充分交互,通过接受外部参数,完成内部动作。

    而且 -v 还可以多重设置,把多个变量输入到 awk 执行代码段之中,这真的是一个很有用功能。如:

    $ awk -v qual=40 -v pos=17662793 '$1!~/^#/ && $6>qual && $2>pos {print $1,$2,$8}' demo.vcf
    chr22 17663076 CMDB_AF=0.053564,CMDB_AC=534,CMDB_AN=9525
    

    在上面这个命令里面,我不但通过自定义参数要求 6 > qual,还同时要求只输出那些2 > pos 的结果。你如果有更多的需要,可以不断往后加上 -v 设置变量。

    数组

    awk 中也有数组的概念和数据组织形式,不过与其说是数组,不如说更像是哈希表,原因是它的数组索引可以不必像通常我们所知的那样。
    首先,它的数组语法格式这样的:

    array_name[index] = value
    

    其中:

    • array_name 是数组的名称;
    • index是数组的索引,这个索引可以是数字下标也可以是字符下标;
    • value是数组中元素的值

    接下来,我们先看一下应该如何创建和访问数组中的元素:

    $ awk 'BEGIN{sites["chrom"]="chr22"; sites["pos"]=17662679; print sites["chrom"], sites["pos"]}'
    

    这个命令执行之后,print出来的结果是:

    chr22 17662679
    

    在上面代码中,我定义了一个名字为 sites 的数组,这个数组的索引下标我不是用通常的数字,而是字符——后面再举例子讲数字下标,这个做法与哈希表如出一辙(或者说,就是哈希)。用字符索引代替数字索引的好处是,可以用名称来获得对应的 value,建立起索引和 value 之间的一个映射关系,甚至可以像哈希表那样通过 index 进行信息查找。

    这个方式还可以 “人为地” 制造出多维数组。只需要你把索引的命名按照多维数组那样的形式来进行就可以。比如,以一个二维数组为例,我们可以用 array_name["0,0"]、array_name["0,1"]、 array_name["1,0"]、array_name["1,1"]分别代表一个 2×2 数组中的各个元素,这里就不额外举例子了。

    以上是字符下标的数组,接着我举一个数字下标的数组例子:

    $ echo  "this is a variant in vcf file"  |  awk  '{split($0, array, " "); for(i=1; i<=length(array); i++){print array[i]} }' this
    is
    a
    variant in vcf file
    

    在这个例子里面,我想你也可以看出来,数字下标的数组一般都是通过文本处理而产生的,比如这里我就是通过 split 函数,把 “this is a variant in vcf file” 这一个字符串,按照空格,将它切分为一个数组,数组中的元素为这字符串中的每一个单词。然后,再写一个循环语句将其输出(循环语句中 length函数,可以获取到该数组的长度),值得注意的一个地方是,awk 数组的第一个元素下标是 1 而不是 0。

    另外,如果要删除掉数组中的某个元素,只需要通过 delete 语句就可以实现,语法:

    delete array_name[index]
    

    这样就可以随意把任意一个 index 索引的元素删除掉。

    其实,awk 的数组功能,我们在生物信息数据分析的场景中用的不多,就算真要用到,这个分析任务的复杂性也往往不是在 awk 仅用数组就可以解决的,这个时候可能也是需要写成脚本的时候了。但不管如何,数组的创建和使用方法还是值得在这里描述清楚的。特别是在数组上也可以有更多的操作,比如,还可以用 asort 对数据元素进行排序,或者使用 asorti 对数组索引进行排序。

    再谈条件判断与循环语句

    awk 虽然是一个 文本文件处理程序,但其实它也像是一个编程语言,所以在常见编程语言中该有的功能和语法表达形式,其实它也照样有。比如,之前提到的 if - else 语句,这里我还要再说上一说,同时也把循环语句补充上来。

    先说 if 的语法:

    if  (条件)  { 动作 }
    

    中间的执行动作,都括在大括号里。由于之前(见上一篇文章)已经给过不少例子了,所以这里我想偷个懒,只要大家能够看明白的,就不多举例子了。

    除了 if 语句,紧接着的就是 if-else 语句,它的语法结构是:

    if (条件) {
      动作
    } else {
      动作
    }
    

    if 中的判断条件符合了,就执行 if 中的动作,否则执行 else 中的动作,这是一个比较常用的语句功能。除了上面两种之外,其实 awk 也有 if-else-if 语句,我们可以用它来创建多个 if-else 组合,实现多条件判断。

    if  (条件1){ 动作 }  else  if  (条件2)  { 动作 }  else  if  (条件3)  { 动作 }  else  { 动作 }
    

    关于 awk 的 if 语句就在这里都补充完成了。接下来说一说,awk 中的另一个重要语句:循环。

    循环也是常规编程语言用有的核心语法,在 awk 中也不例外。虽然,awk 在处理文本数据的时候,BODY 语句会自动循环执行的,但是它的循环是在文本文件中一行行往下进行的循环。如果我们需要在每一行文本处理中都做出一些其他的循环操作,那么就需要使用 awk 提供出来的循环语句。

    awk 的循环语句有两种:for 和 while 。

    对于 for 循环来说,它的语法是这样的:

    for  (起始条件初始化; 终止条件; 迭代起始条件)  { 动作 }
    

    对于有过编程基础的朋友来说,应该对这种结构非常熟悉,几乎所有常见的编程语言,都是类似的for循环结构。它在执行的时候,先初始化起始条件,然后与终止条件比较,如果条件为真,那么执行 for 循环中的动作——也就是执行循环体,然后执行第三部分“迭代起始条件”——这个迭代一般是递增或者递减操作,然后再继续和终止条件进行比较,只要比较结果为真,就一直循环下去;直到条件为假,才终止 for 循环并退出这个执行语句。下面就是一个简单的循环输出数字的 awk 语句:

    $ awk  'BEGIN{ for(i=0; i<4; i++){print i} }' 
    1
    2
    3
    

    之所以把这个语句中用在 BEGIN 里,目的其实就是想省下对具体文件的处理,方便作为例子。至于在具体的项目中,还应该按照具体的文件处理需求来执行。

    对于 while 循环来说,它的语法结构为:

    while  (终止条件)  { 动作 }
    

    相比于 for 循环语句,while 语句要简单得多。它只检查 while 后面的条件是否为真,如果是真,那么执行,如果为假,那么结束循环。这里用数字输出作为例子:

    $ awk  'BEGIN{i=1; while(i<4){print i; ++i;} }' 
    1
    2
    3
    

    在 for 或者 while 循环中,并不是只有等到终止条件为假的时候,才可以退出循环。有时在执行的过程中,我们也可以强制中断循环体或者跳过某一次循环。能够完成这两个功能的是 awk 循环中提供的 break 和 continue 语句,而且这两个都是只在循环体(执行动作的语句)中使用的语句。

    break 语句可以让我们在碰到某个条件的时候就强制退出循环,而 continue 语句则可以让在碰到某个条件之后,直接忽略在 continue 之下的执行动作,直接回到循环头进入下一次循环迭代。比如,我们用 continue 举个例子,输出所有 1-10 之间的奇数:

    $ awk  'BEGIN{ for(i=1; i<=10; i++){ if(i % 2 > 0){print i;} else { continue; }} }'
    1
    3
    5
    7
    9 
    

    自定义函数

    awk 中自定义函数的语句是 function ,使用这个语句,就越来越像是在编程了,虽然能够做的事情更多了,但代价是整个 awk 也会因此变得更加复杂。

    函数的好处,除了功能模块化之外,就是提高代码的复用性。在 awk 中我们自定义函数的语法是:

    function function_name(参数1,参数2,参数3,...){ 动作 }
    

    其实跟前面的语句有类似之处,都是关键字+名称+参数(或者判断条件)+动作的模式。这里函数前面的 function 关键字是必须,它规定了这是一个自定义的函数。其中:

    • function_name 是函数名字;
    • 大括号括起来的一系列执行动作是该函数所要完成的具体功能

    另外,函数的定义一般要在其它 awk 操作之前完成。我自己没有合适的例子,就借用网上的一个 awk 函数来举例吧。下面代码定义了两个功能很简单的函数,它们分别用于数字比较之后,返回数据中的最小值和最大值,然后还定义了一个 main 函数作为主函数来调用它们。而且,一般来说,当需要自定义函数时,代码都会比较长,已经不适合在一行命令中写下,所以会写成一份真正的 awk 脚本文件,这个文件的后缀用 .awk,比如这里我们就可以将其命名为 function_demo.awk ,其中的所有 awk 代码如下:

    # 返回最小值
    function find_min(num1, num2){
      if (num1 < num2)
        return num1
      return num2
    }
    # 返回最大值
    function find_max(num1, num2){
      if (num1 > num2) {
        return num1
      } else {
        return num2
      }
    }
    
    # 主函数
    function main(num1, num2){
      # 查找最小值
      result = find_min(num1, num2)
      print "Minimum =", result
    
      # 查找最大值
      result = find_max(num1, num2)
      print "Maximum =", result
    }
    
    # 整个脚本还是从这里开始执行
    BEGIN {
      main(30, 20)
    }
    

    这时,通过 awk -f 执行这个脚本,我们就可以得到如下结果:

    $ awk -f function_demo.awk
    Minimum = 20
    Maximum = 30
    

    要再提醒大家的是,这个脚本里只定义了 BEGIN 代码段,这是为了可以在不用有任何文件输入时也能执行。但在实际使用的时候,我们是需要定义 BODY 代码段的,甚至还有 END 代码段的,并且在最后还要有一份待处理的文件作为输入。

    还能同时处理多个文件?

    其实从 awk 本来的设计理念来看,它最适合的场景是一次只处理一份文件。但如果在某些情况下,我们非要同时处理多个文件,awk 也能做到,只是这个情况用的很少,而且也相对费劲一些。我自己从未如此使用过,它也不是本文的重点,所以这里我也不打算进一步展开,只是想告诉大家 awk 是有能力这样做的,大家真有需要了,再从网上或者它的手册中找到它的具体用法吧。

    小结

    这篇文章就在这里结束吧。如无意外这应该也是最近两篇 awk 文章中的最后一篇,四千五百多字(不含代码)。看完这一篇,再加上上一篇的 awk 基础用法,我们其实已经可以用 awk 来实现很多工作了,包括很复杂的文本处理,都完全可以通过 awk 实现。但是,我觉得要提醒一下大家,awk 是动态语言,执行效率并不是很高,处理一些比较小的文件,确实没有什么问题。但,如果要处理大型的文件,比如 BAM 之类的,那么不建议用 awk 。而且,awk 的功能毕竟还是比较单一,在处理多文件处理方面也不是很灵活,也不能很好地与其他代码进行交互,更加没有什么基于 awk 开发的包来支持更多的分析,它本身是一把精致的匕首,我们就不要过多地将其它当大刀来使。任何工具或者编程语言都应该是用在它最合适的地方上才好,用不着因为手里拿着一个锤子,所以就要把世界都当成了钉子。对我来说,使用 awk 主要还是图它在基本文本处理方面的简单、方便和快捷,可以只用一行命令就搞定很多事情,如果复杂了我也不一定要用 awk 了。

    参考链接

    http://www.runoob.com/w3cnote/awk-work-principle.html http://www.runoob.com/w3cnote/awk-user-defined-functions.html


    如果喜欢更多的生物信息和组学文章,搜索并关注我的微信公众号“碱基矿工”(ID: helixminer)

    碱基矿工

    你还可以读

    这是我的知识星球:『达尔文生信星球』(原名:解螺旋技术交流圈),是一个我与读者朋友们的私人朋友圈,如今已有超过400人在星球中一起学习和交流。我有9年前沿而完整的生物信息学、NGS领域的科研经历,在该领域发有多篇Nature、Cell级别的科学文章,我希望借助这个知识星球可以与更多的志同道合者沟通和交流,同时也把自己的一些微薄经验分享给更多对组学感兴趣的伙伴们。
    这是知识星球上第一个与基因组学和生物信息学强相关的圈子,也是官方评定的优秀星球。如今已经累计超过1100个主题,希望能够借此营造一个高质量的组学知识圈和人脉圈,通过提问、彼此分享、交流经验、心得等,促进彼此更好地学习生信知识,共同提升基因组数据分析和解读的能力。

    在这里你可以结识到全国优秀的基因组学和生物信息学专家,同时可以分享你的经验、见解和思考,有问题也可以向我提问和星球里的星友们提问。

    达尔文生信星球

    相关文章

      网友评论

        本文标题:生物信息 awk 用法进阶

        本文链接:https://www.haomeiwen.com/subject/rlbkgqtx.html