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/