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/