生意这样的商业活动,理论上可以激励创造,让参与交易的双方都受益,才有可能持续。 比如你不可能花费半年时间去系统性学习R语言和Linux操作,处理fastq的单细胞测序数据,做统计可视化图表,就为了一辈子的一个项目。所以理论上你的最优解决方案是委托专业的生信工程师,比如我们就在单细胞天地发布过:单细胞转录组下游分析这样报价合理吗? ,试图促进大家合作,可惜的是虽然几万人参与讨论,但是最后却不欢而散。专业的工程师觉得为客户学习一个R包收费2000合情合理,但是委托者觉得一个项目全套分析收2000才合理。
正好,最近我的学徒也分享了他在分析他们课题组数据项目的一个感悟,值得分享给大家,大家可以思考一下,为什么会有这样的生意,交易双方都觉得血亏呢?
学徒在文献里面找到了下面这段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代码看懂了,今天来把坑填上。
# 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; \
}'
?代表判断,若为真,则取:前的值;若为假,则取:后的值
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)
}
上面这段代码看起来已经挺复杂了,但它其实只是整篇文章的一个小插曲。如果有兴趣的话,可以看看下面的故事背景:
整个workflow就是以上这么回事,难点在于隐马尔可夫模型(HMM),也就是如何根据我们测到的RPKM,来判断一个bin有没有表达。 (我内心:这不是直接划个阈值就完事了吗???好吧,存在即合理,简单地划个阈值肯定会有诸多问题的) 下面我尽可能通俗地讲一下HMM是怎么回事:
cat $input | segment -e $EMISS_PROB_FILE -z 0.5 -p 0.995 > segment_output.txt
# 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计算转录单元,不过,实在是太冷门了,如果你确实不从事这方面就不需要再看了哈。