那到底该如何因地制宜呢

Posted on by

前面我在教程：标准是需要因地制宜的，提到过根据inferCNV结果进行细胞的拷贝数打分，主要是读取一个文件 infercnv.observations.txt ， 然后根据里面的的数值进行归类，阈值是 0.3，0.7，1.3，1.5，2 分别代表拷贝数缺失或者扩展的程度界限。其中位于 0.7-1.5之间就是拷贝数正常。但是呢，这个标准在Smart-seq2和10X的是有差异的。在Smart-seq2来源的单细胞转录组数据集，我们可以：

cnv_score_table[cnv_score_mat > 0 & cnv_score_mat < 0.3] <- "A" #complete loss. 2pts
cnv_score_table[cnv_score_mat >= 0.3 & cnv_score_mat < 0.7] <- "B" #loss of one copy. 1pts
cnv_score_table[cnv_score_mat >= 0.7 & cnv_score_mat < 1.3] <- "C" #Neutral. 0pts
cnv_score_table[cnv_score_mat >= 1.3 & cnv_score_mat <= 1.5] <- "D" #addition of one copy. 1pts
cnv_score_table[cnv_score_mat > 1.5 & cnv_score_mat <= 2] <- "E" #addition of two copies. 2pts
cnv_score_table[cnv_score_mat > 2] <- "F" #addition of more than two copies. 2pts

但是10X来源的单细胞转录组，如果仍然是采用这个标准，就会出现找不到有拷贝变异的细胞，非常尴尬。

那么，我们就需要探索，这个 阈值是 0.3，0.7，1.3，1.5，2 分别代表拷贝数缺失或者扩展的程度界限，是如何得到的呢？

我们读取一个文件 infercnv.observations.txt 做判断，所以就可以读取infercnv.references.txt来获得阈值。这个时候我想起来了，前两天的教程：R语言的各种统计分布函数，你应该了解的都在这！，提到了sigma法则，对于正态分布的x，x取值在(mean-3sd,mean+3sd)几乎就是极端值啦，因为pnorm(3)-pnorm(-3)=0.9973002，这个概率外的事情基本上不可能发生！

所以我就尝试了一下，读取这个 infercnv.observations.txt 文件，对全部的的数值选取两个极端值，因为 infercnv.observations.txt 里面记录的是我们设定的正常2倍体单细胞的转录组表达量矩阵。

tmp=read.table("inferCNV_output/infercnv.references.txt", header=T)
down=mean(rowMeans(tmp)) - 3 * mean( apply(tmp, 1, sd))
up=mean(rowMeans(tmp)) + 3 * mean( apply(tmp, 1, sd))
oneCopy=up-down
oneCopy
a1= down- 2*oneCopy;a1
a2= down- 1*oneCopy;a2
down;up
a3= up + 1*oneCopy;a3
a4= up + 2*oneCopy ;a4

得到的结果如下：

> a2= down- 1*oneCopy;a2
[1] 0.06005809
> down;up
[1] 0.6905682
[1] 1.321078
> a3= up + 1*oneCopy;a3
[1] 1.951588
> a4= up + 2*oneCopy ;a4
[1] 2.582098

上面的代码如果你想运行，需要follow CNS图表复现09—上皮细胞可以区分为恶性与否，系列教程。

真的是蛮有意思的哦！

​

Comments are closed.