著名“错误使用”案例——Excel的单元格格式

之前看过一篇文章，讲述的是2016年发表在《Genome Biology》上的一项研究报道——Microsoft Excel中的单元格格式的自动转换导致研究的期刊中将近20%的论文存在基因名转换错误。比如会把RIKEN identifier “2310009E13”作为科学计数法转换成“2.31E+19”，或者将基因名“SEPT2”作为日期转换成“2018/09/02”。其实早在2004年该问题就已经在《BMC Bioinformatics》上有相关报道，但似乎没引起警觉。不过Excel作为一款在各种领域广泛使用的电子表格软件，微软不太可能为了这类问题而牺牲其“智能”、“易用”的特色。并且这也只能算是用户的“错误使用”，而不是Excel的bug。下面是两篇论文的原文，有兴趣可以一看：

• http://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-1044-7
• https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-5-80

不得不小心的awk之坑，Caveat Emptor！

以上是联想到的题外话，本文主要介绍的是“GNU awk”（即gawk）中的一个坑。和Excel的故事非常类似，都可以归咎为理解不够深入而导致的“错误使用”。而且相比Excel格式转换这种错误，这个问题隐蔽性更高，即使对编程语言的精度问题有所了解，但是还是会掉到坑里，甚至让我一度以为awk出现了bug。

具体表现

在做生物信息学分析时，我们经常要和各种p-value打交道，比如富集、差异表达、正选择分析。这类结果通常整理成一个包含很多行的表，我们关心其中“显著”的部分，就会经常用到awk进行筛选。awk命令的具体用法我就不做介绍了，下面直接进行问题演示：

现在有一个文件，第一列是行号，第二列为p-value，内容如下：

$ cat test.txt
1    0.2
2    3.44346e-111
3    6.53124151732e-318

如果我们要筛选出“p-value < 0.05的行”，理论上应该输出第2行和第3行，但你会发现使用awk命令只输出了第2行：

$ awk '$2<0.05' test.txt
2    3.44346e-111

这是精度问题吗？

在科学运算中，精度是个非常常见的问题，很多人看到这个异常可能会马上想到它。计算机中很多浮点数值以非精确的二进制的形式保存，比如十进制的“0.1”对应的二进制是一个无限循环的数。我们不妨先做个测试，看看一个明显超出双精度浮点取值范围的数1e-3000在awk中会有什么样的表现：

$ awk 'BEGIN{if(1e-3000==0){print "True"}else{print "False"}}'
True # 1e-3000超出了awk的双精度浮点（64位）限制，近似成了0
$ awk 'BEGIN{if(1e-3000<0.05){print "True"}else{print "False"}}'
True # 很好理解，并不会影响这个近似成0的数和0.05之间进行大小比较的判断

我们直接用文件中的6.53124151732e-318，再来试试awk：

$ awk 'BEGIN{if(6.53124151732e-318==0){print "True"}else{print "False"}}'
False # 说明这个数没有超出双精度的范围（最小正数）
$ awk 'BEGIN{if(6.53124151732e-318<0.05){print "True"}else{print "False"}}'
True # 在直接输入这个数给awk时，也能正确判断这个数小于0.05

也就是说，在awk中，6.53124151732e-318这个数，并没有超出双精度浮点的取值范围，也不会影响到和0.05的大小判断。那为什么数据通过管道或者文件输入会导致awk判断大小出问题，而直接输入到awk中又不会呢？似乎这不是精度问题？下面我将一步一步解释。

双精度浮点到底是怎么储存的？

事实上，不管是awk、Python、R、还是现代计算机系统中其他程序，数值的储存方式都遵循“IEEE 754”标准，其中双精度浮点（binary64）在内存中储存形式如下图所示（图片来自Wiki词条：Double-precision floating-point format）：

• sign为符号位，占1 bit，代表正负；
• exponent为指数位，占11 bits，范围从-1024~1023（2¹¹ = 2048）；
• fraction为有效数字，占52 bits，使得精度约为15~17位有效数字（2⁵³ ≈ 1.11e−16）。

其中保证15~17位“完整有效数字精度”的取值范围见下表:

但为什么6.53124151732e-318这个数，没有超出双精度浮点的取值范围呢？因为在对有效数字精度“妥协”的情况下，最小非零正数可以达到4.94066e−324（2⁻¹⁰²²·2⁻⁵² = 2⁻¹⁰⁷⁴ ≈ 4.9e−324）。

我们在awk中可以做个测试：

$ awk 'BEGIN{print 5e-324}'
4.94066e-324 # 在这个最小非零正数附近的数都会近似为这个数
$ awk 'BEGIN{print 3e-324}'
4.94066e-324 # 在这个最小非零正数附近的数都会近似为这个数
$ awk 'BEGIN{print 2e-324}'
0 # 如果小太多，会近似为零
$ echo 2.226e-308 | awk '{if($1<0.05){print "True"}else{print "False"}}'
True # 没有超过“完整有效数字精度”的数，结果就是正确的
$ echo 2.225e-308 | awk '{if($1<0.05){print "True"}else{print "False"}}'
False # 超过“完整有效数字精度”的数，就会立即引发这个问题

因此，可以肯定，这的确是精度引发的问题。

从管道、文件中读入的数，和直接在awk中输入的有何不同？

我们先来看看，awk中的数据类型是怎么指定的：

$ echo Hello | awk '{print typeof($1)}'
string # 从管道输入的"Hello"很明显是一个字符串类型（typeof函数需要新版本的awk）
$ echo 1990 | awk '{print typeof($1)}'
strnum # 从管道输入的1990这类纯数字组成的，被认为是strnum类型
$ awk 'BEGIN{print typeof(1990)}'
number # 直接输入的1990，是数字类型

那么awk究竟以怎样的规则指定数据类型的呢：

• 数值常量或者数值运算结果指定为是数字类型number
• 字符串常量或字符串运算结果指定为字符串类型string
• 从管道、逐行读入、FILENAME、ARGV、ENVIRON以及由match()、split()、patsplit()产生的数组的元素。为数字的指定为strnum类型，否则为字符串类型string。未初始化的变量也为strnum类型
• 数据类型通过赋值传递给变量。但不管如何使用变量，它的数据类型不再变化

我们再来看看，数据超过“完整有效数字精度”范围后，awk怎么指定数据类型：

$ echo 2.226e-308 | awk '{print typeof($1)}'
strnum # 管道读入的、没有超出“完整精度”范围的数，指定为strnum类型（第3条规则）
$ echo 2.225e-308 | awk '{print typeof($1)}'
string # 管道读入的、超出了“完整精度”范围的数，就被指定为了string类型

那么问题就来了，超出了“完整精度”范围的string，和number类型的0.05之间如何比较大小？这点，awk的官方文档也给出了说明：

也就是说，string和number之间的大小比较，会被当做两个string之间的比较。一切真相大白！

总结

第一个例子中发生了什么？

• 6.53124151732e-318这个数从文件或者管道输入时，因为超过了“完整有效数字精度”的最小非零正数范围，所以被awk当作了string类型
• 成为了string类型的6.53124151732e-318，和number类型的0.05进行比较时，被当成了两个string之间的比较
• string类型的6.53124151732e-318比string类型的0.05大，所以此行不输出

可能有人会怀疑——这种超过“完整有效数字精度”的数，在日常生信分析中真的会出现吗？直接给出例子：

> phyper(940, 1133, 1000, 999, lower.tail=F)
[1] 5.780568e-322 # 常规的统计检验，比如这里的超几何分布，也可能得到这种数
> p.adjust(c(5.780568e-322, 0.02, 0.1), method="fdr")
[1] 1.73417e-321  3.00000e-02  1.00000e-01 # 通过多重检验校正，也可能存在于最后结果中

正因为这类数并没有超出最小非零正数，所以在做统计检验的时候能被正常输出，但是这个数又已经超过了“完整有效数字精度”，所以引发了awk的问题。

那么awk正确的使用姿势是？

$ awk 'strtonum($2)<0.05' test.txt # 通过strtonum()函数把p-value转换为number类型
2    3.44346e-111
3    6.53124151732e-318
$ awk '$2+0<0.05' test.txt # 通过数值运算“+0”，把p-value转换为number类型
2    3.44346e-111
3    6.53124151732e-318
$ awk '$2-0.05<0' test.txt # 类似，通过数值运算转换number类型
2    3.44346e-111
3    6.53124151732e-318

这3种方法都可以，不过strtonum()作为扩展函数，可能不是在所有版本的awk中都支持（没做过调查）。

写在最后

作为最常用的Linux命令之一，我们经常会用到awk做一些简单的格式化、筛选文本信息的工作。如果理解不透彻，极有可能“错误使用”。在做生物信息学乃至其他方向的科学运算时一定要小心避免踩坑，建议课题组和从事数据分析工作的公司仔细排查下现有分析流程中可能存在的问题。

参考资料

• https://www.gnu.org/software/gawk/manual/gawk.html
• https://en.wikipedia.org/wiki/Double-precision_floating-point_format
• https://www.gnu.org/software/gawk/manual/html_node/Variable-Typing.html#Variable-Typing