查看: 884|回复: 37

跟NBT学Scissor| bulk RNA + scRNA鉴定与目标表型相关的细胞亚群

  [复制链接]

版主

Rank: 10Rank: 10Rank: 10

主题
43
注册时间
2016.3.21
在线时间
185 小时

活跃会员突出贡献


发表于 2021.12.21 22:23:43 | 显示全部楼层 |阅读模式

本文首发于“生信补给站”,https://mp.weixin.qq.com/s/SKoaR1YQvbgOUx1Iv069MA



Nature Biotechnology上发表了题为“Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data”的研究文章。研究团队开发了Scissor算法,可利用大量单细胞数据和表型信息识别与疾病高度相关的细胞亚群。



本文会尝试介绍一下(1)Scissor算法的工作流程和部分文章内容 ,以及
(2)使用R实现方法。

一 Scissor算法的工作流程
Scissor的工作流程如图所示

图a:首先需要三个数据文件:单细胞表达矩阵、bulk表达矩阵 和 bulk样本的表型文件。
图b:其次计算每对细胞和bulk样本的Pearson相关系数构建相关系数矩阵(关键步骤) ,然后 Scissor会通过优化样本表型Y与相关矩阵S的回归模型。
图c:由上述优化模型求解的非零系数β用于选择与目标表型相关的细胞亚群,其中Scissor+ 表示所选择的细胞与感兴趣的表型呈正相关,Scissor-为负相关。表型信息可以是连续变量、二分向量或临床生存数据,会分别对应不同的回归方法。
图d,e :为控制假阳性后续会有一个可靠性显著性检验以及 差异表达基因分析、功能富集分析和motif分析等下游分析。
以上,Scissor就可以利用单细胞表达数据,bulk表达矩阵和表型信息,从单细胞数据中自动识别与给定表型相关度最高的细胞亚群。 二 文章TCGA-LUAD数据集结果

文章使用Scissor在模拟数据集,TCGA肺腺癌数据集,TP53 表型, 黑色素瘤scRNA数据集 , 阿尔茨海默症(AD)数据集的测试中均有很好的表现。

下面简单的介绍下TCGA-LUAD数据集的结果,
如图j 所示通过计算得到了361 Scissor+ cell 和534 Scissor - cell,对照图i 发现大部分为cancer cell ,只有零星几个是T cell 和 B cell 。
图k给了详细的数值,其中Scissor+ cell 中有98.1%为cancer cell 。

表明在有表型数据和bulk seq数据的情况下,Scissor能够很好地区分肿瘤细胞和正常细胞,从单细胞数据中准确识别与目标表型相关的大部分细胞。 注意Scissor+ 和 Scissor- 是相对的,取决于表型结局中0 和 1的设置。
更多数据集的结果请详细阅读原文(https://www.nature.com/articles/s41587-021-01091-3)。

三 R实现Scissor算法
3.1 安装,加载
下载Scissor包和示例数据

  1. #Scissor安装和加载
  2. devtools::install_github('sunduanchen/Scissor')
  3. library(Scissor)

  4. #加载LUAD 的scRNA-seq count数据
  5. location <- "https://xialab.s3-us-west-2.amazonaws.com/Duanchen/Scissor_data/"
  6. load(url(paste0(location, 'scRNA-seq.RData')))

  7. #加载LUAD 的bulk-seq 数据 和 生存数据
  8. load(url(paste0(location, 'TCGA_LUAD_exp1.RData')))
  9. load(url(paste0(location, 'TCGA_LUAD_survival.RData')))
复制代码
3.2 查看数据情况(1)查看三个数据集的基本情况
  1. dim(sc_dataset)
  2. #[1] 33694  4102

  3. dim(bulk_dataset)
  4. #[1] 56716   471

  5. bulk_dataset[1:4,1:4]
  6. #       TCGA-05-4249 TCGA-05-4250 TCGA-05-4382 TCGA-05-4384
  7. #TSPAN6     57.52398   66.4940573    30.905120     35.10434
  8. #TNMD        0.00000    0.1748819     0.000000      0.00000
  9. #DPM1       98.83813  135.2046911    84.362971     89.10553
  10. #SCYL3      16.23339    6.0050820     6.226792     16.93050

  11. head(bulk_survival)
  12. #  TCGA_patient_barcode OS_time Status
  13. #1         TCGA-05-4249    1158      0
  14. #2         TCGA-05-4250     121      1
  15. #3         TCGA-05-4382     607      0
  16. #4         TCGA-05-4384     426      0
  17. #5         TCGA-05-4389    1369      0
  18. #6         TCGA-05-4390    1126      0
复制代码
(2)准备一个seurat数据( 最好包含预处理数据和构造网络数据)
  1. sc_dataset <- Seurat_preprocessing(sc_dataset, verbose = F)

  2. names(sc_dataset)
  3. #[1] "RNA"     "RNA_nn"  "RNA_snn" "pca"     "tsne"    "umap"   

  4. #umap可视化
  5. DimPlot(sc_dataset, reduction = 'umap', label = T, label.size = 10)
复制代码


(3)检查bulk seq数据 和 表型数据表型文件的第一列为示例id,注意其顺序需要与bulk表达矩阵中的列名顺序相同。此处以生存数据为例,表型文件的第二列是一个二元变量,'1'表示事件(如癌症复发或死亡),'0'表示右截尾 (0,1对应最终的Scissor- Scissor+)。
  1. #检查对应顺序
  2. all(colnames(bulk_dataset) == bulk_survival$TCGA_patient_barcode)
  3. #TRUE

  4. phenotype <- bulk_survival[,2:3]
  5. colnames(phenotype) <- c("time", "status")
  6. head(phenotype)
  7. #  time status
  8. #1 1158      0
  9. #2  121      1
  10. #3  607      0
  11. #4  426      0
  12. #5 1369      0
  13. #6 1126      0
复制代码

3.3 运行 Scissor
此处为生存数据,所以family = "cox" ,如果是logistic回归则使用family = "binomial"。
  1. infos1 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = 0.05,
  2.                   family = "cox", Save_file = 'Scissor_LUAD_survival.RData')
复制代码


可以看到Scissor首先输出相关系数的五分位数的summary,相关系数均为正且值都不接近于0。如果自己使用的数据集的中值相关性过低(< 0.01),则scissors会给出警告,表型-细胞关联的结果可能不太可靠。

3.4 查看Scissor结果Scissor 鉴定出预后较差的 201 Scissor+ cells 和预后较好的 4个 Scissor- cells ,可以通过infos1数据中查看这205个cell 的具体信息

  1. names(infos1)
  2. #[1] "para"        "Coefs"       "Scissor_pos" "Scissor_neg"

  3. length(infos1$Scissor_pos)
  4. #[1] 201
  5. infos1$Scissor_pos[1:4]
  6. #[1] "AAAGTAGAGGAGCGAG_19" "AACCATGCATCTCCCA_19" "AACCGCGAGCTGCGAA_20" "AACTCAGTCCGCGGTA_19"

  7. length(infos1$Scissor_neg)
  8. #[1] 4
  9. infos1$Scissor_neg
  10. #[1] "ACGCCAGTCCTCCTAG_20" "ACGGGCTAGTGGCACA_20" "CCGGTAGGTACCCAAT_15" "GACGCGTAGTGGTCCC_20"
复制代码
3.5 Scissor 结果 umap可视化
  1. #Scissor+ Scissor- 定义
  2. Scissor_select <- rep("Background cells", ncol(sc_dataset))
  3. names(Scissor_select) <- colnames(sc_dataset)
  4. Scissor_select[infos1$Scissor_pos] <- "Scissor+ cell"
  5. Scissor_select[infos1$Scissor_neg] <- "Scissor- cell"

  6. #metadata 中添加 Scissor信息
  7. sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "scissor")
  8. DimPlot(sc_dataset, reduction = 'umap', group.by = 'scissor', cols = c('grey','royalblue','indianred1'),pt.size = 2)
复制代码


可以发现clusters 1 和 3  是与预后差表型 最相关的细胞亚群,还可以做一下差异分析以及富集分析等后续分析。

3.6 优化模型参数(1)alpha 设置为NULL
  1. infos2 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = NULL, cutoff = 0.03,
  2.                  family = "cox", Load_file = 'Scissor_LUAD_survival.RData')
复制代码
(2)alpha 设置梯度
  1. infos3 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = seq(1,10,2)/1000, cutoff = 0.2,
  2.                  family = "cox", Load_file = 'Scissor_LUAD_survival.RData')
复制代码
更多参数,模型详见 https://sunduanchen.github.io/Sc ... issor_Tutorial.html


参考资料:
[1] https://sunduanchen.github.io/Sc ... issor_Tutorial.html[2] Sun D, Guan X, Moran AE, et al. Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat Biotechnol. 2021;10.1038/s41587-021-01091-3


PS:有个交流的讨论组,想沟通交流的,后台回复”入群“。
◆ ◆ ◆  ◆ ◆
精心整理(含图PLUS版)|R语言生信分析,可视化(R统计,ggplot2绘图,生信图形可视化汇总)

               







本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2021.12.24 08:09:34 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
76 小时

发表于 2021.12.24 10:34:28 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
31
注册时间
2016.1.8
在线时间
556 小时

发表于 2021.12.24 13:40:59 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
31
注册时间
2016.1.8
在线时间
556 小时

发表于 2021.12.24 13:42:10 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2021.12.25 12:11:11 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2021.12.26 09:31:46 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
76 小时

发表于 2021.12.26 16:02:28 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.9.3
在线时间
1 小时

发表于 2021.12.26 22:01:15 | 显示全部楼层
太喜欢了
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2021.12.27 08:57:36 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2021.12.28 10:28:34 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2021.12.29 09:31:54 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2021.12.30 08:10:39 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2021.12.31 17:50:56 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
76 小时

发表于 2022.1.1 22:05:11 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2018.11.22
在线时间
14 小时

发表于 2022.1.3 17:57:33 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
76 小时

发表于 2022.1.3 18:36:29 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2022.1.4 08:17:46 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
76 小时

发表于 2022.1.4 14:38:38 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
79 小时

发表于 2022.1.5 09:17:57 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表