Featured

Linux shell/ python/ R: 生物信息的三把斧头

生物信息是一门典型的交叉学科,主要涉及分子生物学(基因组学)/ 计算机 / 统计。生物信息主要处理大规模测序数据,本质上是纯文本的字符串。人类基因组的参考序列就是一个由ATCG四种碱基排列的巨长字符串(3G bases)。很多时候, 生物信息分析要解决的问题是:
1、给定一条序列,就是一个短一点的字符串,比如50bp,150bp或800bp,需要知道该子字符串在总字符串中的位置,用生物的语言描述就是其在基因组中的染色体坐标。
2、给定很多条短序列,需要拼接成比较长的序列。因为参考基因组往往极长,不大可能通过一次实验从头测定到尾,所以只能通过构建大片段的重叠群,再把大片段继续打断成更小的片段进行测序,通过多重覆盖从而拼接成完整的基因组。这也是一种典型的Divide and Conquer策略。

相对于计算机科学可以处理的类型丰富多样的数据类型,比如图像,视频,声音,网页等等,生物信息的处理对象相对要纯粹和简单的多。但是,生物信息和其他计算机科学相比有个很大不同。就是生物信息的软件开发和数据分析绝大部分时候都是在Linux操作系统下进行。不同于主流的Windows系统电脑操作主要靠在图形界面下点击鼠标,Linux操作系统用的是纯命令行(在终端中敲命令)。生物信息的这个特点使得很多初学者入门比较困难。这就好比不会滑雪的人,即便给了他滑雪工具,他也会操作的很笨拙。但是一旦真正掌握滑雪工具,就会在雪地环境如履平地,挥洒自如。Linux操作系统就是这样方便和强大的系统,能提供异常丰富的文字处理命令,还可以通过管道(|)将前一个命令的输出传给下一个命令做输入,通过这样的命令组合可以实现非常强大的功能。此外,Linux系统还有一个重要特性是系统非常干净和简洁,不像Windows那么臃肿笨拙,那么消耗内存。

所以说,生物信息的从业人员的必备技能是要熟悉Linux操作系统,熟练使用常用的几十个命令,深刻理解Linux用户权限和文件管理。因为往往需要几百G甚至数T的数据,生物信息的从业人员还学要了解并行计算,熟悉常用的job提交系统比如Oracle Grid Engine (从前的SGE)。Linux命令的组合除了可以直接在终端里面一边敲一遍看结果外,也可以写成一个脚本(shell 脚本)以.sh为后缀名命名,然后进行执行。这个shell脚本在现有工具(软件)齐备的情况下可以完成大部分的原始数据预处理工作。数据分析的需求往往是个性化,因为实验设计不大可能千篇一律,否则科研的创新无从谈起。个性化的数据分析就意味着需要写程序,这个时候python就上场了。

当然,编程的语言五花八门,不过生物信息领域常用的主要有:python/ perl/ java/ C/ C++,这些语言都各有自己的用户和支持者。一般来讲,底层的算法比如大规模的序列比对或全基因组组装,这些任务计算量极大,往往需要用C/ C++来写,因为C是编译后执行的而且可以直接操控硬件,使得同样的算法用C/ C++实现时速度最快。但是,底层算法往往很早就被生物信息领域的顶尖专家(比如bwa的作者Heng Li)写好了,一个普通的生物信息从业人员需要写底层算法的机会并不常有。更多的时候,他需要写一些程序去收集,提取和总结shell脚本预处理好的数据和结果,在这个时候python可以非常好的解决此类问题。python是一种跨平台的解释型脚本语言,所有的Linux系统都自带python。而且,安装了python语言的Window电脑,也可以在CMD中运行python脚本。python的可读性强,学习时更容易上手。现在好多课堂都在尝试将python作为主要的教学语言,这个教学语言过去曾经是pascal/ basic/ VB/ C。此外python在科学计算,机器学习和深度学习领域都有非常强的社区支持。众多基于python开发的模块和程序包,都可以通过pipconda等工具轻松安装和卸载。写程序本质上是像搭积木一样的活,只有站在前人工作的基础上灵活使用别人的包才能提高程序的开发效率。

最后一把必备的工具是R。R其实是脱胎于S-plus的统计编程语言,有强大的统计分析和数据可视化(画图)的功能。但python也有一些可以做统计分析和画图的模块(如matplotlib),但是统计分析和可视化是R最擅长的领域。R的语法相对比较奇怪和不够直观,因此能用进入R console里面敲几个命令画画图并不难,但是要掌握着这么语言写大段的代码就没有那么容易了。不过,R的强大之处在于很多统计方法都封装好了,因此往往几行代码就可以实现非常酷的功能。R有非常强大的社区支持,在CRAN和Bioconductor上可以下载安装很多功能各异的算法包。

生物信息的这三把斧头各有分工,一个生物信息的数据分析项目最终所有的软件(python or R)和命令都可通过一个shell脚本一以贯之,从头跑到尾,成为一个流程 (pipeline)。

再论武汉市到底还有多少新型冠状病毒感染者?

首先这个病毒的传染力惊人从公开的新闻媒体报道有3个独立的案例可以说明这一点。

案例一:国家卫建委专家组成员北大第一医院呼吸和重症学科主任王广发主任被感染。虽然医护人员是高危人群,但是我们有理由相信王主任的专业素质和防范措施是到位的。即便这样,王主任还是因公感染了。

案例二:合肥6名年轻人因元月21日聚会,其中一人系19日从武汉回来,21日和同学聚会,导致之后6人全部被诊断为新型冠状病毒肺炎病人。

案例三:河南安阳鲁某,有武汉居住史,元月10日回安阳,一直没有症状却感染了3名亲属,并最终导致累计有5人感染。

案例二和案例三都充分说明此次新型冠状病毒可以导致无症状患者对外感染。这种超强的传播力和隐蔽性,为该病毒的防控提出了巨大考验。

武汉市到底有多少感染者?

先看一个新闻。元月29日也就是昨天,日本终于绷不住派了专机到武汉撤侨。第一架撤侨专机载206名日本人,结果5人有发烧咳嗽症状,其中3人被确诊为新型冠状病毒肺炎。我们根据日本这个小样本估算目前的感染比例是万150左右,这是非常惊人的数字。当然,小样本可能是非常不准确的。

几个方面的证据表明武汉市确诊病人被低估(武汉市应该还有大量尚未确诊的病人)。

证据一:从上海的小样本数据窥探湖北地区受感染的武汉人与武汉以外湖北人的构成。根据上海市卫计委1月30的播报,截止昨天(29日)上海市累计确诊101例。其中按居住地分,湖北武汉的有34人,湖北其它地方的有14人,显示二者比为2.4。须知这个是从湖北扩散到上海的一个小样本,对于我们窥探湖北的情况非常有价值。然而,根据腾讯新闻疫情实时播报,到2020年1月30日下午4点截止,全国确诊7783例,疑似12,167例,治愈129,死亡170。武汉确诊病人与武汉以外的湖北确诊病人的比为2261: 2325,居然小于1。这是不合理的,应该是武汉市的确诊病人被低估。

证据二:从武汉1月23日封城之前2周左右武汉进出人口数量估算武汉感染病人所占全国确诊病人的比例。因为武汉是该病毒的爆发之唯一源头,武汉以外的感染者不论是在湖北还是在其它省份都是因为武汉接触史而受感染的。之前新闻披露武汉封城之前共有500万人离开武汉,姑且不论这其中大部分人是取道武汉去往湖北其它各个城市。我们姑且认为武汉市还有1000万常住人口留守。那么武汉的确诊病人应该是武汉以外确诊病人的2倍。但是,根据今天的数字这个比例是2261:5522,居然小于0.5。表明武汉确诊病人严重不足。

证据三:根据湖北省内某市估算发病率从而对比武汉市的发病率。我们选择距离武汉200多公里的荆门市,该市有290万人口(湖北总人口5724万),有长荆铁路和武汉相通,是荆门人往来武汉最主要的出行方式。荆门和武汉之间每天有6对左右火车往返,一趟火车满载2000余人。这样计算每天有1万人左右往返荆门武汉,计算武汉封城前2周,大概有14万人。同时,武汉有100万大学生,如果其中70%为湖北籍学生,根据人口比例荆门籍的大学生应为100万*290/5724*70%=3万左右。我们估算和武汉有接触的荆门人口为15万左右,为荆门人口的5%。荆门目前确诊191例,则感染率为万12.7。而武汉地区的感染率为2261/1060万,仅为万2.1。这两个数字有近6倍的差异。

根据以上三个方面的观察,我们认为武汉市的感染者完全不足以被目前的确诊人数所充分代表。我们预见在未来的日子里,武汉市还会有更多病人被确诊,直到所有统计数据能够在逻辑上相互吻合为止。