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/

秒客网