copykat为什么没有infercnv直观呢

2021-12-04 10:28:54 浏览数 (2)

其实 copykat 仅仅是算法判别的时候不如人意,但是可视化的时候仍然是肉眼可以明显区分二倍体正常细胞和非整倍体的癌症细胞,所以我们想看看具体做什么改进,可以绕过这个bug,首选项我们把全部的上皮细胞按照病人进行了拆分,得到如下所示 的每个病人独立的文件夹以及每个文件夹下面的expFile.txt !

代码语言:javascript复制
   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脚本 :

代码语言:javascript复制
ls -d S*|while read id;do
echo $id;
cd $id;
nohup Rscript.exe ../scripts/step4-run-copykat.R &
cd ../
done

很容易看到了每个病人的单细胞水平细胞表达量矩阵跑这个copykat流程耗费时间:

代码语言:javascript复制
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

而且我简单看了看这些病人的二倍体正常细胞和非整倍体细胞的判别差异:

代码语言:javascript复制
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的文件:

代码语言:javascript复制
 11G ./4-3-only-Epi
478M ./Rdata-celltype
 81G ./epi-one-by-one

0 人点赞