有一种生意双方都觉得亏
生意这样的商业活动，理论上可以激励创造，让参与交易的双方都受益，才有可能持续。
比如你不可能花费半年时间去系统性学习R语言和Linux操作，处理fastq的单细胞测序数据，做统计可视化图表，就为了一辈子的一个项目。所以理论上你的最优解决方案是委托专业的生信工程师，比如我们就在单细胞天地发布过：单细胞转录组下游分析这样报价合理吗？ ，试图促进大家合作，可惜的是虽然几万人参与讨论，但是最后却不欢而散。专业的工程师觉得为客户学习一个R包收费2000合情合理，但是委托者觉得一个项目全套分析收2000才合理。
正好，最近我的学徒也分享了他在分析他们课题组数据项目的一个感悟，值得分享给大家，大家可以思考一下，为什么会有这样的生意，交易双方都觉得血亏呢？
起初是求助一段perl代码
学徒在文献里面找到了下面这段perl代码，不太看得懂，根据大致猜测，翻译成了后面这段R代码，请曾老师帮忙看一下猜得对不对：
# calculate a discrete logarithmically-stepped integer index from RPKM values
INDEX_FUNCTION='function getIndex(rpkm){ \
ff=log('MAX_RPKM'/'MIN_PRKM') \
i=rpkm/'$MIN_RPKM'; \
i=i<1?0:i; \
i=i>xrpkm?xrpkm:i; \
i=i>0?log(i)/log(ff):-1; \
i=int(i)+1; \
return i; \
}'
getIndex(rpkm)=function(rpkm){
i=rpkm/MIN_RPKM
ff=log(MAX_RPKM/MIN_PRKM)
if(i<1){i=0}
if(i>xrpkm){i=xrpkm}
if(i>0) {log(i)/log(ff)=log(i)/log(ff)-1} #化简之后相当于i=rpkm/MIN_RPKM/ff
i=floor(i)+1
return(i)
}
其实，如果不是擅长这两种编程语言，这两个代码看起来都是天书。
我还没有来得及帮忙解决这个问题，第二天，学徒主动发给我了他的解决方案！
是隐马尔可夫模型划分基因组单元
昨天的perl代码看懂了，今天来把坑填上。
昨天的问题来自一个叫 segment.pl 的软件，用来将连续的基因组划分成离散的片段。划分的依据是用户给定的一些特征，比如新生RNA转录情况，或者CNV。 perl代码：
# calculate a discrete logarithmically-stepped integer index from RPKM values
MIN_RPKM=0.001 # the minimum RPKM value; bins with RPKM<$MIN_RPKM all have index=0
MAX_RPKM=100 # the maximum RPKM value; bins with RPKM>$MAX_RPKM all have the same maximum index
OBS_FOLD_FACTOR=2 # observation indices stepped logarithmically every $OBS_FOLD_FACTOR-fold from $MIN_RPKM to $MAX_RPKM
#以上这些是软件文档里的默认参数和解释
INDEX_FUNCTION='function getIndex(rpkm){ \
xrpkm='$MAX_RPKM'/'$MIN_RPKM'
ff=OBS_FOLD_FACTOR
i=rpkm/'$MIN_RPKM'; \
i=i<1?0:i; \ #若i<1,则i=0;否则i=i，即不变
i=i>xrpkm?xrpkm:i; \ #若i>xrpkm,则i=xrpkm;否则不变
i=i>0?log(i)/log(ff):-1; \ #若i>0，则i=log(i)/log(ff);否则i=-1
i=int(i)+1; \ #i向下取整后加1
return i; \
}'
?代表判断，若为真，则取:前的值；若为假，则取:后的值
翻译成R：
MIN_RPKM=0.001
MAX_RPKM=100
OBS_FOLD_FACTOR=2
getIndex=function(rpkm){
xrpkm=MAX_RPKM/MIN_RPKM
ff=OBS_FOLD_FACTOR
i=rpkm/MIN_RPKM
if(i<1){i=0}
if(i>xrpkm){i=xrpkm}
if(i>0) {i=log(i,base=ff)}
else {i=-1}
i=floor(i)+1
return(i)
}
代码作用
我花了九牛二虎之力终于明白这段代码的作用了： （在这之前已经把基因组划分成了1kb的bin，计算了每个bin的nascent RNA表达量RPKM） 将每个bin的RPKM换算成index，其换算方法为：若RPKM低于设定的下界，则index=0；若RPKM高于设定的上界，则令RPKM=设定的上界 然后对于RPKM不低于下界的bin，做以下计算：先将RPKM处以设定的下界，然后以ff变量为底取对数，再向下取整后加1 用人话翻译一下：对在一定范围内的表达量，取对数后进行等分，从而将表达量转为index（正整数）；表达量超过范围时，将index设定为index的最小或最大值。 原文是这么描述这段代码的作用的： Corrected bin data were then indexed into a series of bounded integer observation values from 0 through 17, logarithmically spaced across bin RPKM values from <0.0005 to >100. 以及这么描述的： Score the bins by rounding each into one of 17 logarithmically distributed RPKM input states. 想通了以后觉得作者这么写也没问题，但是对于不懂的人来说就是天书啊orz
故事背景——用隐马尔可夫模型对基因组进行分割
上面这段代码看起来已经挺复杂了，但它其实只是整篇文章的一个小插曲。如果有兴趣的话，可以看看下面的故事背景：
目的：切割基因组
需求：将连续的基因组分割成离散的转录单元。
输入：全基因组每个1kb的bin的RPKM 核心：隐马尔可夫模型（等下讲） 直接输出：每个bin的转录状态，即index（有限范围的正整数）。 间接输出：index高于阈值的bin认为是转录的，index低于阈值认为是不转录的 最终输出：相邻的转录的bin连接起来，构成转录单元（transcription unit）。这是nascent RNA测序中一个重要的概念。
整个workflow就是以上这么回事，难点在于隐马尔可夫模型（HMM），也就是如何根据我们测到的RPKM，来判断一个bin有没有表达。
（我内心：这不是直接划个阈值就完事了吗？？？好吧，存在即合理，简单地划个阈值肯定会有诸多问题的）
下面我尽可能通俗地讲一下HMM是怎么回事：
做法：隐马尔可夫模型
隐马尔可夫模型的四要素：真实状态（state），观察到的状态（observed value），发射概率（emission possibility），转换概率（transition possibility） 对于nascent RNA来说，测序测到的表达量是observed value，而我们想知道的转录单元是state。为什么要用HMM？是因为测序（observed value）不可能完完全全准确地反映真实状态（state），所以要用HMM来帮助我们判断真实状态。 发射概率—— 在某一种真实状态下，产生某一种观察状态的概率，它是一个条件概率：P ( value=j | state=i )。在我们这个问题中，发射概率是这样算的： 先把每个bin都注释到唯一的基因（如果没注释到，或者注释到多个基因，则舍去这个bin）。 同时把每个bin和每个gene的RPKM都转换成index（也就是开头的那段perl代码了） 然后计算index=j的bin落在index=i的基因里的频率。i和j都遍历各自的所有可能，然后就会得到一个条件概率矩阵，这个矩阵就是emission possibility file。
转换概率—— 从上一个单位到下一个单位时（在这里就是上一个bin到下一个bin），真实状态从a变成b（也就是转录状态发生改变）的概率。作者设定了转换概率为0.005，也就是说一个转录本有0.995的概率会继续转录下去，有0.005的概率会终止。反之，一个不转录的区域有0.005的概率会开始转录。 四大要素都准备好了之后，就可以用segment软件来计算了。（至于其中的算法，用的是Viterbi，应该只有理工科学生才会去关注了）
# input文件中是每个bin的信息，包含三列，分别是index，chr，起始位置
# EMISS_PROB_FILE中是发射概率矩阵
# -z 0.5 代表全基因组的bin中，假设不转录的bin占50%
# -p 0.995 代表转换概率为0.005
以上就是隐马尔可夫模型的故事。我相信它这么复杂，没有几个人会用它来研究转录问题的。
所以做nascent RNA和transcription unit研究的人也不是很多，国内应该也没有公司会做这个测序。
毕竟吃力不讨好，看懂这个模型花了我将近两周时间，写这个文档来分享都花了一个小时，谁乐意跳这个大坑呢？
你现在觉得专业的工程师为了你的课题帮忙看一个算法看一个R包或者perl代码，收费多少合适呢？
你以为这样就完了吗？实际上还有进阶资料，关于HMM计算转录单元，不过，实在是太冷门了，如果你确实不从事这方面就不需要再看了哈。
文末友情宣传
全国巡讲全球听（买一得五） ，你的生物信息学入门课 生信技能树的2019年终总结 ，你的生物信息学成长宝藏 2020学习主旋律，B站74小时免费教学视频为你领路
HMM计算转录单元进阶资料
可试读80%