扫码联系客服
“吾生也有涯，而知也无涯。”
大家好，我是生信人新来的小编青壬。坐标魔都东北，博士生，从事生物高通量数据分析工作。
我自己也关注了大量的生信方面的公众号。作为一个新兴领域，微信上各大公众号提供了诸多关于生信领域分析的见解与指导。所提供的信息也是多样的，小到一个热图怎么画，大到整套数据分析的流程。不过真正想入门生信乃至掌握生信，不经过长期的身体力行是不可能的。其中比较重要的是编程关，因为长期GUI界面使得大家习惯于以处理多种任务的可能性去换取操作的便捷性（两者就是一个跷跷板类的trade-off关系），甚至即使迁移到命令行也希望能拥有一个万金油脚本应对数据，每次跟着跑就可以了。这些观念往往让人在编程学习的进程上停滞不前。我将首先不定期更新关于编程的一些心得系列，希望对于广大初学者有所启示。
最近有闲的时候刷Rosalind上的生物编程题，其中有一道是关于斐波拉契数列的（有兴趣的同学可以到最下看一看想一想原题）。其实本来并不是太想做这题，毕竟这是一个经典的数学模型了，高中数学用一些技巧已经可以解出其通项公式。维基百科介绍如下：
至于递推公式的数学证明，只需要对Fn状态下的兔子略作分类讨论即可，在此不再赘述。
然而本着除恶务尽的负责态度，还是不情愿打开了这一题：
仔细一看还挺有意思，这里的兔子并不像之前经典模型中生存期无限长，而是会在3个月后死亡。【因为2个月就没法玩了...】
受到之前数学思维的影响，我还是倾向于解通项公式。然而在这种稍加改动的斐波拉契数列中，这样的分类讨论变得非常困难。我在简书上也看到某同学的做法（如果认头像没错的话这个同学应该是生信之光的小编，这个公众号的文章也不错），但很不幸他的做法是错的，因为其总结出来的递推公式甚至超过第三项就对不上实际数字了。
我觉得可以从后往前寻找递推公式，于是灵机一动就写出来了：
这个式子的意思是，你能活到现在的兔子，难道不是前三个月每一次新生兔子的和吗？
不过正当我兴冲冲去验证的时候，发现同样并不能符合实际模型【至于为什么这是一个坑，大家稍想一下应该可以明白】。这样问题又回到了起点。
于是我搬来了救兵，稍久他给出了证明...
这么朴素的证明我实在是很难读懂...于是发生了这样的Q&A(这只是一部分)。
当然事实证明这个方法是对的，利用了轮换对称的思路，最后所得的递推公式与实际情况完全相符。可是最后我转念一想，我们面对的是计算机，并不需要给出数学上的通项公式呀！于是一个简单的小脚本诞生了：
#!/usr/bin/env python3# -*- coding: utf-8 -*-
def count(n):
for i in range(n-1):
a = res[0]
b = res[1]
c = res[2]
res[0] = b + c
res[1] = a
res[2] = b
return(sum(res))
with open('3M-survival.txt','w') as f:
for i in range(1,31):
res = [1,0,0]
f.write("\t".join((str(i),str(count(i)),'\n')))在编程过程中我们只是描述了一下每个月发生“新陈代谢”的过程而已。因为兔子只能活三个月，那么像回转寿司一样只能占着对应的坑“流动”，这也就是a,b,c三分类的由来。最后利用for循环从头计算即可，毕竟你不能一步步算的，计算机完全可以呀！比如六个月输出的结果就是这样的：
和题目中的结果完全一致。 接下来我兴致来了，想进一步探究四个月，五个月生存期的模型种群数量的指数关系，这就很好办了，一样的思路，输出文件画个图拟合一下即可。R都不需要，在Excel做好了：
def count(n):
for i in range(n-1):
a = res[0]
b = res[1]
c = res[2]
d = res[3]
res[0] = b + c + d
res[1] = a
res[2] = b
res[3] = c
return(sum(res))
with open('4M-survival.txt','w') as f:
for i in range(1,31):
res = [1,0,0,0]
f.write("\t".join((str(i),str(count(i)),'\n')))
def count(n):
for i in range(n-1):
a = res[0]
b = res[1]
c = res[2]
d = res[3]
e = res[4]
res[0] = b + c + d + e
res[1] = a
res[2] = b
res[3] = c
res[4] = d
return(sum(res))
with open('5M-survival.txt','w') as f:
for i in range(1,31):
res = [1,0,0,0,0]
f.write("\t".join((str(i),str(count(i)),'\n'))) 最终可以得到如下的曲线图：
当然我还拿了标准的斐波拉契曲线做一下对比：
总结一下
通过这个问题应该明白，在我们编程解决实际问题的过程中，并不应该被传统的数学思维所限制。计算机在很多情况下是愿意做苦行僧的，利用这样特质可以有效减轻我们自己的思维负担。
期待下一次再见。
参考链接：
Rosalind原题（http://rosalind.info/problems/list-view/）
生信之光小编的解答（http://www.jianshu.com/p/01513bd10a41）
点击以下「关键词」，查看往期内容：
TCGA | 小工具 | 数据库 |组装| 注释 | 基因家族 | Pvalue
基因预测 |bestorf | sci | NAR | 在线工具 | 生存分析 | 热图
生信不死 | 初学者 | circRNA | 一箭画心| 十二生肖 | circos
舞台|基因组 | 黄金测序 | 套路 | 杂谈组装 | 进化 | 测序简史
