ONCOCNV軟件思路分析之tumor處理

前期處理

• perl腳本統計RC(RC(read counts))
• 讀入control baseline 和 sigma（最後baseline 預測的mad值）
• 將gc < 0.28或gc > 0.68，sigma乘上1.5,後來又乘以6，對於小於0.01或者大於0.99分位數，sigma取0.01和0.99分位點的sigma
• 將sigma轉化爲權重，SigmaForWeights = 1/sigma^2/max(1/sigmaforWeithts^2)
• 根據mu值設置一些outlier的amplicon，threshold爲-2和2
• 若是1/3 amplicon都是NA，該樣品被拋棄

文庫大小校訂，（中位數校訂），GC含量校訂，長度校訂，ICA校訂

• 文庫大小校訂，（中位數校訂），GC含量校訂，長度校訂，多了一步中位數校訂，將NRC除以總的NRC值的中位數,做用是一開始對於0拷貝數進行校訂
• ICA標準化：取control樣品中的ICA1，ICA2，ICA3使用rlm創建線性模型求殘差做爲標準化 後的值，對於值爲0的logNRC，取最小logNRC

segmentByCBS

• 輸入標準化後的各個amplicon的染色體，起始位置，logNRC和權重（SigmaForWeights）,使用PSCBS中segmentByCBS（內部調用DNAcopy包）函數進行片斷劃分。劃分時，sigma越大權重越小，來避免斷點斷定在sigma較大的amplicon。輸出是每一個segment對應的拷貝數

輸出以下，若是該染色體中存在NA的數據，會保留一行NA的值

sampleName chromosome     start       end nbrOfLoci    mean
1        <NA>          1   2488068 244006378      1363  0.0934
2        <NA>         NA        NA        NA        NA      NA
3        <NA>          2   5833035 234681120      1008  0.0732
4        <NA>         NA        NA        NA        NA      NA
5        <NA>          3   3192502 195622096       939  0.0718
6        <NA>         NA        NA        NA        NA      NA
7        <NA>          4   1800963 153332857       423  0.0313
8        <NA>         NA        NA        NA        NA      NA
9        <NA>          5    223515 180058652       574 -0.3460
10       <NA>         NA        NA        NA        NA      NA
11       <NA>          6    393195 167275553      1378  0.0588
12       <NA>         NA        NA        NA        NA      NA
13       <NA>          7   2946229 152373106       945  0.0507
14       <NA>         NA        NA        NA        NA      NA
15       <NA>          8  30915961 145742416       852  0.0803
16       <NA>         NA        NA        NA        NA      NA
17       <NA>          9   5021975   5072501        21 -0.1446
18       <NA>          9   5072501 134100903       474 -0.3635
19       <NA>          9 134100903 139438447       156  0.0474

predictCluster,，肯定zero(copy number 0對應的logNRC)和copy number。注：這裏求出的copy number並非生物意義的copynumber 而是是否出現拷貝數異常，0爲未出現拷貝數異常

• 將每一個segment內的amplicon乘以權重求中位值獲得權重後獲得的weighted.median segmean，若是沒法求出weighted.median，取wighted.mean做爲segmean
• 原理：對應着相同copy數的segmean應該聚類在一塊兒
• 由於直接取segmean進行聚類可能會聚類錯誤，所以對數據添加高斯噪音：
   第一步，求出須要添加的sdError：取出在-2和2之間的amplicon（而且排除x和y染色體），將amplicon的logNRC - segmean獲得errors，對每一段segment的errors求方差，獲得該段sdNoise,再將sdNoisd比上該段amplicon長度開根獲得sdError，公式爲：sdError = sd(logNRC - segmean)/sqrt(n)，能夠理解爲sdError = ((logNRC - segmean) - mean(logNRC - segmean))^2/sqrt( n(segment)*n(in each segment) )
   第二步，添加噪音，噪音的均值爲0,方差爲sdError，每一個ampliocn 的值a = segmean + rnorm(n(in each segment),0,sdError)
• 使用mcluser聚類
  若是聚類結果大於1類，將該全部類別的方差取出，將異常簇給值NA，異常cluster判斷方法：每一個cluster都有一個對應的標準差（sigma），標準差大於中位cluster sigma加上7*mad（sigma）
• 計算每一個cluster的density
  計算方法：計算該cluster（高斯模型）的density，等於每一個cluster（所有高斯模型）的均值在這個cluster下的density的總和再乘以這個cluster的density比重。而每一個cluster的density = dnorm(x,mean of cluster,var of cluster) * proportion of cluster
• 將最大density的cluster做爲zero copy number
• 將最近的一個cluster（若是mean值差<0.04）兩個cluster乘比重後的均值賦給該cluster(zero)，該cluster比重變爲這兩個cluster相加
  
• 檢測預測時處於邊緣的值
  若是同一個amplicon a 對應的copy number mad值爲0.5,那麼a所對應的amplicon copy number 爲1，若是同一個amplicon a 對應的copy number mad值爲-0.5,那麼a所對應的amplicon copy number 爲-1
  若是常染色體只有1個cluster，那麼分析x,y染色體，若是segmean < log(0.75)，該amplicon copy number 給-1，若是segmean > log(1.25)，amplicon copy number給1
• 求出maxloss 和 maingain
  maxloss:
  copy number是-1的amplicon的segmean，
  copy number 是0 的amcplion segmean -(copy number是0的amplicon segmean - mean(copy number 是 -1的 amplicon segmean ))/2
  在上述中取最大值
  mingain:
  copy number是1的amplicon的segmean，
  copy number 是0 的amcplion segmean +(mean(copy number 是 1的 amplicon segmean )-copy number是0的amplicon segmean)/2
  在上述中取最小值
• 將x和y染色體segmean > mingain ，copy number賦予1，將x和y染色體segmean < maxloss，copy number 賦予-1
• 對predictCluster進行八次，取最小的zero值的一次，而後取出zero和預測的copy number

計算生物意義的copy number

• ratio = logNRC - zero，減去zero值(copy number 0 的cluster的均值)
• 若是存在正常細胞污染的比例，那麼減去這部分的影響: ratio = log((exp(ratio)-normalContamination)/(1-normalContamination))
• 求出copy number 0的amplicon求mad值，做爲所有的標準差(sample sigma)，全部amplicon的control sigma * sample sigma做爲校訂後的方差
• 將每一個segment內的amplicon ratio乘以權重求中位值獲得權重後獲得的weighted.median segmean，若是沒法求出weighted.median，取wighted.mean做爲segmean
• copy number計算方法
  

按照segment進行fixed variance test和t test

• fixed variance test
  ratio 應服從N(0,sigma)，對segment內的amplicon ratio 進行轉換，(logNRC-0)/sigma應該服從N(0,1)，根據中心極限定理，mean((logNRC-0)/sigma)服從N(0,1/sqrt(n))，計算該segment的mean((logNRC-0)/sigma)在N(0,1/sqrt(n))雙尾的p值
• t test
  (logNRC-0)/sigma應該服從N(0,1)，對每一個amplicon進行t test
  若是fixed variance test > 0.01 或者t test > 0.01，copy number賦予2
  CNV：copy >= 3 : Gain，copy = 2.5：PotentialGain，copy < 1：Loss，copy = 1.5：PotentialLoss

檢測segment的outlier

(ratio- segmean)/sigma在N(0,1)雙尾的p值
ratio/sigma在N(0,1)雙尾的p值
取二者之中較大的p值做爲該amplicon的outlier p值
全部的amplicon outlier值採用fdr方法計算q值
對於outlier q值<0.01的點加上註釋，outlier p值，outlier copy number採用round(exp(values[tt])*4)/2

按照基因進行t test

• 對每一個基因使用t檢驗,理論上 ratio/sample sigma/control sigma，服從均值爲0的分佈，若是p值<0.01，對每一個基因copynumber進行估計(perGeneEvaluation)
• 若是基因計算的結果不等於segment的結果，而且斷點不在基因內，ratio - segmean / (control sigma * sample sigma)應服從1/sqrt(n)正態分佈，做t檢驗，若是t檢驗顯著，predLarge Corrected給perGeneEvaluation
• 若是斷點在基因內，以第一個segmean(fragratio)做爲參考,若是出現segmean不同，將下一個amplicon segmean值改成參考segmean
• 對每一個基因內的ratio - segmean(fragratio)/(control sigma * sample sigma),做t檢驗，取出最大的p值做爲基因內全部amplicon的pvalRatioCorrected
• 但若是基因內存在copy number = 2 的amplicon，判斷是否存在round(exp(fragratio)*4)/2 = 2，若是有從這些amplicon中取出p值最大的p值做爲基因內全部amplicon的pvalRatioCorrected
• 若是該p值>0.05，說明所有點都被該fragratio解釋了
• 若是p值<0.05，說明不能所有解釋，對每一個fragratios相同的amplicon取出，將它們ratio/(control sigma * sample sigma)做t檢驗，若是最小值>0.01，那麼pergeneEvaluation=2

參考資料

ONCOCNV文獻：https://academic.oup.com/bioinformatics/article/30/24/3443/2422154
ONCOCNV代碼：http://boevalab.com/ONCOCNV/