有一種生意雙方都覺得虧

  • 2020 年 3 月 11 日
  • 筆記

生意這樣的商業活動，理論上可以激勵創造，讓參與交易的雙方都受益，才有可能持續。 比如你不可能花費半年時間去系統性學習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，應該只有理工科學生才會去關注了）
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計算轉錄單元，不過，實在是太冷門了，如果你確實不從事這方面就不需要再看了哈。

VirMach 便宜 VPS

QNews

QNews