跟着一区文章学通单细胞（十一）：（拓展篇）如何区分肿瘤组织中的癌上皮和正常上皮？（上）

跟着一区文章学通单细胞（十一）：（拓展篇）如何区分肿瘤组织中的癌上皮和正常上皮？（上） by 生信作曲家

前言

生信作曲家一直致力于构建国内和谐、开放、可持续的生物信息交流平台。生信作曲家的目标是为每个科研人扫清科研路上的一切障碍，让生物信息学人人可做。目前，生信作曲家已推出"跟着一区文章学通机器学习"，"跟着一区文章学通单细胞"，CIBERSORTx, NTP, TICPE，metaVIPER算法等多种强力教程，同时也兼具GEO数据库，ICGC数据库数据整理等核心入门技术教学，开办居家防疫7天发文A计划，生信中秋分享会、小年夜分享会、暑期集训营和冬季集训营等多次。

众所周知，单细胞系列教程在全网已经很多。相比于其他教程，生信作曲家希望为大家带来基于科研问题的，带着科研思考的单细胞转录组全套教程。这是本教程最具特色的之处，也是为什么你需要学习的原因。

本套教程来源于这篇今年2021年12月发表在Oncogene上的文章。doi: 10.1038/s41388-021-02054-3。之所以选择oncogene而不是CNS级别的刊物，是因为一向来的实践告诉我们，往往这些口碑较好的一区文章的思路、技术、格局，恰恰是我们普通科研人能最大化汲取的。

下辑的内容包括亚群细分，通路分析

拓展篇开始，我们将学习单细胞一些高级分析，不过，我们更希望你带着思考学习这些分析方法，而不是只为了“高级”而分析。

同时，值得一提的是，很多高级分析实际上并不一定非常合理！如果你带着思考去做，往往会发现一些很值得探讨的探讨的问题

本节内容带来的是一个命题：如何识别肿瘤样本里的上皮细胞哪些是真癌细胞，哪些是正常上皮，长时间不分析后，你需要对你的数据进行一个比较详细的回忆！因此本节恰好我们共同回忆上下辑的内容。

这是一个来自3个肺癌病人的单细胞数据集（3正常VS. 3肿瘤)，我们已经把他们拆分成了免疫、上皮、间质三部分，本节重点看上皮（epi）的组分。如果你对数据集陌生，请回顾第一节到第十节的内容（上下辑）

数据 阅读原文，密码pi8t获取。

代码实操

# 加载包
library(seurat)
library(ggplot2)
library(Seurat)
library(dplyr)
library(reticulate)
library(sctransform)
library(cowplot)
library(viridis)
library(tidyr)
library(magrittr)
library(reshape2)
library(readxl)
library(stringr)
library(cowplot)
library(scales)
library(tibble)
library(gplots)
library(RColorBrewer)

theme_set(theme_cowplot())

#配色方案
use_colors <- c(
  Tumor = "brown2",
  Normal = "deepskyblue2",
  G1 = "#46ACC8",
  G2M = "#E58601",
  S = "#B40F20",
  Epithelial = "seagreen",
  Immune = "darkgoldenrod2",
  Stromal = "steelblue",
  p032 = "#5B1A18",
  p033 = "#9C964A",
  p034 = "#FD6467",
  CNN = "chartreuse4",
  CNA = "orange")

epi <- readRDS("epi.RDS")
epi <- RunPCA(epi)
ElbowPlot(epi,  ndims = 50)


epi <- RunUMAP(epi, dims = 1:20)
epi <- FindNeighbors(epi, dims = 1:20)
for (i in c(0.2, 0.3, 0.4, 0.5, 1, 2)) {
  epi <- FindClusters(epi, resolution = i)
  print(DimPlot(epi, reduction = "umap", label = T) + labs(title = paste0("resolution: ", i)))
}

### 我们依然选择比较大的resolution，手动注释没有关系
### 控制在20个cluster左右即可
Idents(epi) <- epi@meta.data$SCT_snn_res.1

### 细胞类型注释对于上皮而言我们并不喜欢，因为很多上皮类型我们几乎是没听过的
### 我们更关注哪些是肿瘤细胞,我们要记住我们曾经认为，normal里比例高的细胞就是normal
### tumor里比例高的就是tumor
### 大家可以运行下面的代码回忆一下
  epi_clusters <- FetchData(epi, vars = c("SCT_snn_res.1", "tissue_type"))

  count_tumor <- epi_clusters %>% filter(tissue_type == "Tumor") %>% count() %>% as.numeric()
  count_normal <- epi_clusters %>% filter(tissue_type == "Normal") %>% count() %>% as.numeric()
  ### 每个cluster计数总数
  epi_counts <- epi_clusters %>% group_by(tissue_type) %>% count(SCT_snn_res.1)
  ### 除以总count数，是朴素的比例计算
  proportion_tumor <- epi_counts %>% filter(tissue_type == "Tumor") %>% mutate(proportion = n/count_tumor)
  proportion_normal <- epi_counts %>% filter(tissue_type == "Normal") %>% mutate(proportion = n/count_normal)

  ### 组合，最后一行很有意思，认为哪群比例高就属于哪群
  proportion_epi <- full_join(proportion_normal, proportion_tumor, by = "SCT_snn_res.1") %>% 
    mutate(proportion.x = ifelse(is.na(proportion.x), 0,  proportion.x)) %>%  
    mutate(proportion.y = ifelse(is.na(proportion.y), 0,  proportion.y)) %>%
    mutate(tissue_type.x = "Normal") %>%
    mutate(tissue_type.y = "Tumor") %>%
    mutate(cluster_type = ifelse(proportion.x > proportion.y, "Normal", "Tumor"))

  ### 每个细胞的详细归属cluster
  cluster_type_data <- left_join(x = epi_clusters, y = proportion_epi, by = "SCT_snn_res.1")
  rownames(cluster_type_data) <- rownames(epi_clusters)

  ### 加入metadata中
  epi <- AddMetaData(epi, select(cluster_type_data, cluster_type))


  ### 柱状图

  n1 <- select(proportion_epi, c(tissue_type.x, SCT_snn_res.1, proportion.x)) %>%
    mutate(tissue_type = tissue_type.x) %>% 
    mutate(proportion = proportion.x) %>%
    mutate(tissue_type.x = NULL) %>%
    mutate(proportion.x = NULL)
  t1 <- select(proportion_epi, c(tissue_type.y, SCT_snn_res.1, proportion.y)) %>%
    mutate(tissue_type = tissue_type.y) %>% 
    mutate(proportion = proportion.y) %>%
    mutate(tissue_type.y = NULL) %>%
    mutate(proportion.y = NULL)

  proportion_epi2 <- rbind(n1, t1)

  ggplot(proportion_epi2, aes(fill = tissue_type, y = proportion, x = SCT_snn_res.1)) +
    geom_bar(stat = "identity", position = "dodge", width = 0.7) +
    scale_fill_manual(values = use_colors)

  ggsave2("正常-肿瘤上皮比较.png", path = "result_3v3/", width = 40, height = 20, units = "cm")

  ### 根据下面的图，我们可以明确地识别出肿瘤上皮中哪些类群是新增加的
  ### 也就是异型增生的上皮细胞——癌细胞
  DimPlot(epi, group.by = "cluster_type",split.by = 'tissue_type' ,cols = use_colors, pt.size = 0.1)
  ggsave2("新增的肿瘤细胞.png", path = "result_3v3/", width = 8, height = 4)

### 但是方法的好处是快捷
### 但是坏处也很突出：在于理论上肿瘤细胞正常里应该没有，但是这边却还是有几个，这与分群的细节有关
### 因此上面这种比例法供大家备选

其他还有什么方法剥离出肿瘤上皮（癌变上皮）呢，拷贝数推断inferCNV可能是值得学习的方法，但是我们建议你运行inferCNV之前先领会本节的方法

后记

本节希望你能关注到区分肿瘤和正常上皮这个需求，并且后续想办法解决他！

ixxmu / mp_duty

跟着一区文章学通单细胞（十一）：（拓展篇）如何区分肿瘤组织中的癌上皮和正常上皮？（上） #2232

跟着一区文章学通单细胞（十一）：（拓展篇）如何区分肿瘤组织中的癌上皮和正常上皮？（上） by 生信作曲家

前言

代码实操

后记