copykat为什么没有infercnv直观呢
2023-04-18 14:28:08 时间
其实 copykat 仅仅是算法判别的时候不如人意，但是可视化的时候仍然是肉眼可以明显区分二倍体正常细胞和非整倍体的癌症细胞，所以我们想看看具体做什么改进，可以绕过这个bug，首选项我们把全部的上皮细胞按照病人进行了拆分，得到如下所示 的每个病人独立的文件夹以及每个文件夹下面的expFile.txt ！
81M 12 2 10:09 S11/expFile.txt
123M 12 2 10:22 S21/expFile.txt
84M 12 2 10:41 S47/expFile.txt
75M 12 2 10:54 S50/expFile.txt
78M 12 2 11:05 S53/expFile.txt
95M 12 2 11:20 S56/expFile.txt
95M 12 2 11:34 S57/expFile.txt
78M 12 2 11:49 S58/expFile.txt
107M 12 2 12:01 S63/expFile.txt
77M 12 2 12:17 S65/expFile.txt
87M 12 2 12:29 S66/expFile.txt
87M 12 2 12:42 S69/expFile.txt
79M 12 2 12:56 S71/expFile.txt
75M 12 2 13:08 S72/expFile.txt
76M 12 2 13:19 S74/expFile.txt
78M 12 2 13:34 S78/expFile.txt
77M 12 2 13:45 S79/expFile.txt
77M 12 2 13:57 S81/expFile.txt
84M 12 2 14:09 S82/expFile.txt
这样我只需要写一个简单的脚本 step4-run-copykat.R ，就可以并行全部的样品的方式进行曲线救国，也算是并行计算。
这个 简单的脚本 step4-run-copykat.R 下面的shell脚本 ：
ls -d S*|while read id;do
echo $id;
cd $id;
nohup Rscript.exe ../scripts/step4-run-copykat.R &
cd ../
done
很容易看到了每个病人的单细胞水平细胞表达量矩阵跑这个copykat流程耗费时间：
S11/nohup.out:Time difference of 54.33337 mins
S21/nohup.out:Time difference of 1.20906 hours
S47/nohup.out:Time difference of 54.88378 mins
S50/nohup.out:Time difference of 43.44914 mins
S53/nohup.out:Time difference of 51.72247 mins
S56/nohup.out:Time difference of 1.020941 hours
S57/nohup.out:Time difference of 1.001784 hours
S58/nohup.out:Time difference of 52.24233 mins
S63/nohup.out:Time difference of 1.096883 hours
S65/nohup.out:Time difference of 45.49075 mins
S66/nohup.out:Time difference of 56.40917 mins
S69/nohup.out:Time difference of 56.40654 mins
S71/nohup.out:Time difference of 52.58897 mins
S72/nohup.out:Time difference of 42.06179 mins
S74/nohup.out:Time difference of 47.09917 mins
S78/nohup.out:Time difference of 49.3231 mins
S79/nohup.out:Time difference of 50.86545 mins
S81/nohup.out:Time difference of 44.29561 mins
S82/nohup.out:Time difference of 45.59241 mins
而且我简单看了看这些病人的二倍体正常细胞和非整倍体细胞的判别差异：
tail -n 6 S*/nohup.out
==> S11/nohup.out <==
aneuploid diploid
epi 177 0
ref-Bcells 7 389
ref-Tcells 0 478
spike-Bcells 6 235
spike-Tcells 0 289
==> S21/nohup.out <==
aneuploid diploid
epi 935 23
ref-Bcells 9 393
ref-Tcells 0 480
spike-Bcells 4 241
spike-Tcells 0 289
==> S47/nohup.out <==
aneuploid diploid
epi 195 13
ref-Bcells 9 385
ref-Tcells 0 477
spike-Bcells 6 232
spike-Tcells 0 289
==> S50/nohup.out <==
aneuploid diploid
epi 27 8
ref-Bcells 316 77
ref-Tcells 20 457
spike-Bcells 184 55
spike-Tcells 24 265
==> S53/nohup.out <==
aneuploid diploid
epi 74 27
ref-Bcells 0 393
ref-Tcells 0 477
spike-Bcells 1 238
spike-Tcells 0 289
==> S56/nohup.out <==
aneuploid diploid
epi 348 114
ref-Bcells 0 399
ref-Tcells 0 477
spike-Bcells 0 241
spike-Tcells 0 289
==> S57/nohup.out <==
aneuploid diploid
epi 104 356
ref-Bcells 0 395
ref-Tcells 0 477
spike-Bcells 0 239
spike-Tcells 0 289
==> S58/nohup.out <==
aneuploid diploid
epi 96 2
ref-Bcells 311 82
ref-Tcells 16 461
spike-Bcells 180 60
spike-Tcells 18 271
==> S63/nohup.out <==
aneuploid diploid
epi 605 80
ref-Bcells 0 400
ref-Tcells 1 477
spike-Bcells 0 243
spike-Tcells 0 289
==> S65/nohup.out <==
aneuploid diploid
epi 56 28
ref-Bcells 307 84
ref-Tcells 32 445
spike-Bcells 190 45
spike-Tcells 29 257
==> S66/nohup.out <==
aneuploid diploid
epi 266 19
ref-Bcells 16 379
ref-Tcells 0 478
spike-Bcells 9 234
spike-Tcells 2 287
==> S69/nohup.out <==
aneuploid diploid
epi 290 3
ref-Bcells 7 391
ref-Tcells 0 479
spike-Bcells 5 238
spike-Tcells 0 289
==> S71/nohup.out <==
aneuploid diploid
epi 71 52
ref-Bcells 0 393
ref-Tcells 0 477
spike-Bcells 0 240
spike-Tcells 0 289
==> S72/nohup.out <==
aneuploid diploid
epi 28 2
ref-Bcells 317 75
ref-Tcells 29 448
spike-Bcells 188 47
spike-Tcells 32 254
==> S74/nohup.out <==
aneuploid diploid
epi 48 25
ref-Bcells 0 393
ref-Tcells 0 477
spike-Bcells 1 238
spike-Tcells 0 289
==> S78/nohup.out <==
aneuploid diploid
epi 94 8
ref-Bcells 317 75
ref-Tcells 35 442
spike-Bcells 192 44
spike-Tcells 31 255
==> S79/nohup.out <==
aneuploid diploid
epi 46 53
ref-Bcells 313 80
ref-Tcells 28 449
spike-Bcells 183 57
spike-Tcells 31 258
==> S81/nohup.out <==
aneuploid diploid
epi 75 3
ref-Bcells 309 83
ref-Tcells 21 456
spike-Bcells 187 49
spike-Tcells 20 266
==> S82/nohup.out <==
aneuploid diploid
epi 70 15
ref-Bcells 83 308
ref-Tcells 455 22
spike-Bcells 49 186
spike-Tcells 259 27
这样就很有意思，可以看到是有一半的病人里面的单细胞出现非常诡异的判别，比如这个S82/nohup.out ，它里面的大量的ref和spike都被判定为恶性细胞。让我们看看它的图表：
test_copykat_heatmap
确实有点奇特哦，这个时候肉眼跟软件其实是一样的！而且我去看了它的inferCNV结果，如下所示：
infercnv
可以看到，copykat 仅仅是没有infercnv直观，但是在这样的恶性细胞比例不高的病人数据里面，确实效果上没有太多区别，跟肉眼判断细胞恶性与否的结论也比较吻合！copykat 虽然把大量的ref和spike错误的判定为恶性细胞，但是很明显我们看图就会反过来把前面的恶性上皮细胞定义为正常细胞。这个时候根据有一些唯心主义的嫌疑了。
虽然 copykat 仅仅是没有infercnv直观，但是copykat至少给出来了 aneuploid 和 diploid的判断，inferCNV给出来的结果文件，仍然是需要自己读取，自己计算cnv打分，自己去给一个阈值来区分细胞的恶性与否。这个步骤也很考验大家的编程功底。至少需要熟读下面的5个教程：
- CNS图表复现13—使用inferCNV来区分肿瘤细胞的恶性与否
- CNS图表复现14—检查文献的inferCNV流程
- CNS图表复现15—inferCNV流程输入数据差异大揭秘
- CNS图表复现16—inferCNV结果解读及利用
- CNS图表复现17—inferCNV结果解读及利用之进阶
如果你也有自己的肿瘤相关单细胞数据的拷贝数变异分析需求， 而且自己多方尝试总是报错，欢迎留言委托我们哈。肿瘤相关单细胞数据的拷贝数变异单项分析2400，如果需要前面的降维聚类分群需要再多加800，提供区分样品的独立copykat和infercnv流程，以及合并的copykat和infercnv流程的全部结果和代码文件夹。
需要注意的是，我们并不会保留中间的txt等文件，因为它实在是太耗费磁盘空间了，比如我们演示的这个项目就20多个病人的不到7000个上皮细胞，走这个copykat和infercnv流程流程，得到了接近100G的文件：
11G ./4-3-only-Epi
478M ./Rdata-celltype
81G ./epi-one-by-one
相关文章
- scrapy框架爬虫
- 云原生场景下实现编译加速
- 读Java8函数式编程笔记07_设计和架构的原则
- 《RPC实战与核心原理》学习笔记Day14
- React中使用Redux (一) - 在React中直接使用Redux
- vue3移动端适配的解决方案
- 函数柯里化详解
- 读Java8函数式编程笔记06_Lambda表达式编写并发程序
- 洛谷oj题单【入门1】顺序结构-入门难度（Java）
- pycharm远程调试、开发（详细操作）
- 线程池的极简用法——内置线程池multiprocessing
- 若依前后端分离框架去掉首页 登录后跳转至动态路由的第一个路由
- SSL 证书基本概念扫盲
- springboot~openfeign开启熔断之后MDC为null的解决
- Zookeeper集群
- 运行项目报错 proxy error: could not proxy request...
- 读Java8函数式编程笔记05_数据并行化
- Google谷歌浏览器Post请求预见strict-origin-when-cross-origin跨域问题的 解决办法
- springboot~openfeign开启熔断之后MDC为null的理解
- elementUI使用el-upload上传文件写法总结及避坑，上传图片/视频到本地/服务器以及回显+删除