单细胞数据复现-肺癌文章代码复现1

原创

小胡子刺猬的生信学习123

发布于 2022-05-06 10:14:33

2.4K2

发布于 2022-05-06 10:14:33

由于老板给定的研究方向是做单细胞组学的，跟以往的组学都是不太一样的，导致学的相当吃力，同时网上的一部分的代码是分段的，有的时候学习不到相关的分析的思路，因此看到这篇文章，作者把全部的代码都整理到一个网站上了，然后提供给相关的研究人员进行后面的复现，真是太优秀了，必须打call，复现好香，自己不用调试代码，自己最近调试到要疯，希望自己变成哪吒，三个脑子。

doi：10.1038/s41388-021-02054-3

代码及数据下载网站：https://codeocean.com/capsule/8321305/tree/v1

因为作者的R代码分为7个大块，最后一个部分是选用的TCGA上的数据进行的再次验证，因为我是搞植物的，TCGA我也用不到，所以最后一个部分我会对代码进行拆解，看一下作者的分析思路，就放弃了进行相关的复现。

首先是对这篇文章的内容进行大致的拆解。

Introduction

肺癌在诊断后治愈性较低，属于大部分疾病中死亡率较高的疾病。目前的组学分析不能在肿瘤微环境的层面解析细胞的多样性，因此，选用单细胞组学对其解析成为评估单个细胞的临床相关性的一个重要依据。

因此，本研究中作者通过单细胞rna的测序技术，解析了肺癌肿瘤微环境中的肿瘤上皮细胞和相关的恶性肿瘤细胞，丰富了以前的单细胞组学研究。通过研究结果获得患者间的肿瘤微环境的异质性细胞组成遵循与癌细胞分化分级相关的特定模式，可以在单细胞的层面上揭示微环境模式的潜在预后相关性。因此，作者得出结论。对肺癌微环境的全面解析可能有助于揭示基于癌细胞和微环境特征的新的临床相关的肿瘤亚型。

代码解析

基本上作者所处理的第一部分的内容主要是针对于结果的第一部分，是基本的单细胞的RNA分析的流程。

首先是将所需要的R包进行加载，对后面要绘图的元素进行颜色的选取及加载，保证全文的配色统一，作者选用的配色相对而言，是比较好看的，我们可以考虑在后面进行摘取，选用。

同时发现作者在前面就对过滤条件进行了配置，其实我们在自己调试参数的时候，可以考虑在后面过滤的时候进行设置。

环境需要的包进行加载

// ### load libraries
library(Seurat)
library(dplyr)
library(reticulate)
library(sctransform)
library(cowplot)
library(ggplot2)
library(viridis)
library(tidyr)
library(magrittr)
library(reshape2)
library(readxl)
library(progeny)
library(readr)
library(stringr)

nFeature_lower <- 500
nFeature_upper <- 10000
nCount_lower <- 1000
nCount_upper <- 100000
pMT_lower <- 0
pMT_upper <- 30
pHB_lower <- 0
pHB_upper <- 5

theme_set(theme_cowplot())

#color scheme
use_colors <- c(
  Tumor = "brown2",
  Normal = "deepskyblue2",
  G1 = "#46ACC8",
  G2M = "#E58601",
  S = "#B40F20",
  Epithelial = "seagreen",
  Immune = "darkgoldenrod2",
  Stromal = "steelblue",
  p018 = "#E2D200",
  p019 = "#46ACC8",
  p023 = "#E58601",
  p024 = "#B40F20",
  p027 = "#0B775E",
  p028 = "#E1BD6D",
  p029 = "#35274A",
  p030 = "#F2300F",
  p031 = "#7294D4",
  p032 = "#5B1A18",
  p033 = "#9C964A",
  p034 = "#FD6467")

data加载

下面就是对数据进行加载。看到作者是选用的cellranger后的三个gz文件进行的分析，其实通过项目经验来看，我们在cellranger那里可选择的东西很少，但是还是有些小细节需要注意的。

第一个就是需要保证我们的参考基因组和注释文件的注释高度完整，这样可以保证后面得到的东西很全。

第二个就是在进行umi定量的时候，也是需要不断进行调整的，来得到自己最适的细胞数量。

这篇文献的作者是将所有的sample整合到一个excel表里面，选用的for循环，进行文件的读取；如果后面我们要用的话，需要保证的是路径正确及sample的表格尽量与作者对应一致。

同时在医学的单细胞分析的时候有线粒体reads、核糖体reads和血细胞reads，在前面过滤的时候都需要考虑进行，对后续的分析是不是有影响。

//# Data loading and QC

### sample list
samples <- read_excel("../data/metadata/patients_metadata.xlsx", range = cell_cols("A:A")) %>% .$sample_id

### import cellranger files from different data sets
for (i in seq_along(samples)){
  assign(paste0("scs_data", i), Read10X(data.dir = paste0("../data/cellranger/", samples[i], "/filtered_feature_bc_matrix")))
}

### create seurat objects from cellranger files
for (i in seq_along(samples)){
  assign(paste0("seu_obj", i), CreateSeuratObject(counts = eval(parse(text = paste0("scs_data", i))), project = samples[i], min.cells = 3))
}

### merge data sets
seu_obj <- merge(seu_obj1, y = c(seu_obj2, seu_obj3, seu_obj4, seu_obj5, seu_obj6, seu_obj7, seu_obj8, seu_obj9, seu_obj10, seu_obj11, seu_obj12, seu_obj13, seu_obj14, seu_obj15, seu_obj16, seu_obj17, seu_obj18, seu_obj19, seu_obj20), add.cell.ids = samples, project = "lung")

### calculate mitochondrial, hemoglobin and ribosomal gene counts
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^MT-", col.name = "pMT")
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^HBA|^HBB", col.name = "pHB")
seu_obj <- PercentageFeatureSet(seu_obj, pattern = "^RPS|^RPL", col.name = "pRP")

qcparams <- c("nFeature_RNA", "nCount_RNA", "pMT", "pHB", "pRP")
for (i in seq_along(qcparams)){
  print(VlnPlot(object = seu_obj, features = qcparams[i], group.by = "orig.ident", pt.size = 0))
}
for (i in seq_along(qcparams)){
  print(RidgePlot(object = seu_obj, features = qcparams[i], group.by = "orig.ident"))
}

VlnPlot(seu_obj, features = c("nFeature_RNA", "nCount_RNA", "pMT"), pt.size = 0, group.by = "orig.ident", ncol = 1, log = T)
ggsave2("SuppFig1B.pdf", path = "../results", width = 30, height = 20, units = "cm")

去除多余的环境变量

//可以将自己不需要的环境变量移除
remove(seu_rna1)   ##文章里面是有17个病人，移除的较长，因此我只是拿其中的一个示例

Data Filtering

作者为了看刚刚在前面提到了三类reads对整体过滤的影响，先用散点图查看自己一开始设的阈值是不是符合要求。

qc_std_plot_helper <- function(x) x + 
  scale_color_viridis() +
  geom_point(size = 0.01, alpha = 0.3)

qc_std_plot <- function(seu_obj) {
  qc_data <- as_tibble(FetchData(seu_obj, c("nCount_RNA", "nFeature_RNA", "pMT", "pHB", "pRP")))
  plot_grid(
    
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), log2(nFeature_RNA), color = pMT))) + 
      geom_hline(yintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_hline(yintercept = log2(nFeature_upper), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), log2(nFeature_RNA), color = pHB))) + 
      geom_hline(yintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_hline(yintercept = log2(nFeature_upper), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), log2(nFeature_RNA), color = pRP))) + 
      geom_hline(yintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_hline(yintercept = log2(nFeature_upper), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), pMT, color = nFeature_RNA))) + 
      geom_hline(yintercept = pMT_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pMT_upper, color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), pHB, color = nFeature_RNA))) + 
      geom_hline(yintercept = pHB_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pHB_upper, color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nCount_RNA), pRP, color = nFeature_RNA))) + 
      geom_vline(xintercept = log2(nCount_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nCount_upper), color = "red", linetype = 2),
    
    
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nFeature_RNA), pMT, color = nCount_RNA))) + 
      geom_hline(yintercept = pMT_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pMT_upper, color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nFeature_RNA), pHB, color = nCount_RNA))) + 
      geom_hline(yintercept = pHB_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pHB_upper, color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_upper), color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(log2(nFeature_RNA), pRP, color = nCount_RNA))) + 
      geom_vline(xintercept = log2(nFeature_lower), color = "red", linetype = 2) +
      geom_vline(xintercept = log2(nFeature_upper), color = "red", linetype = 2),
    
    qc_std_plot_helper(ggplot(qc_data, aes(pRP, pMT, color = nCount_RNA))) + 
      geom_hline(yintercept = pMT_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pMT_upper, color = "red", linetype = 2),
    qc_std_plot_helper(ggplot(qc_data, aes(pRP, pMT, color = nFeature_RNA))) + 
      geom_hline(yintercept = pMT_lower, color = "red", linetype = 2) +
      geom_hline(yintercept = pMT_upper, color = "red", linetype = 2),
    
    
    ggplot(gather(qc_data, key, value), aes(key, value)) +
      geom_violin() +
      facet_wrap(~key, scales = "free", ncol = 5),
    
    ncol = 3, align = "hv"
  )
}

下面是根据一开始设置的阈值开始进行过滤，可以发现在这篇文章的代码中，作者经常会对RDSdata进行保存，方便后面进行不断的参数调试，我最近也是在各种调试参数，才明白这个是有多重要。

// ## Before filtering

seu_obj_unfiltered <- seu_obj

#saveRDS(seu_obj_unfiltered, file = "seurat_objects/all_unfiltered.RDS")

qc_std_plot(seu_obj_unfiltered)

ggsave2("SuppFig1A.png", path = "../results", width = 30, height = 30, units = "cm")


## After filtering

#seu_obj_unfiltered <- readRDS("all_unfiltered.RDS")

seu_obj <- subset(seu_obj_unfiltered, subset = nFeature_RNA > nFeature_lower & nFeature_RNA < nFeature_upper & nCount_RNA > nCount_lower & nCount_RNA < nCount_upper & pMT < pMT_upper & pHB < pHB_upper)

qc_std_plot(seu_obj)

seu_obj_unfiltered
seu_obj

红线区域表示作者准备在后面保留的区域。

可以考虑在过滤后再次用小提琴图或者散点图看一下过滤后的质量。一般是如果后面想要分析一些比较稀有的细胞亚群的时候，建议前面过滤的不要那么狠，后面看一下分群的效果，前面过滤的时候在进行优化。

VlnPlot(seu_obj, features = c("nFeature_RNA", "nCount_RNA", "pMT"), pt.size = 0, group.by = "orig.ident", ncol = 1, log = T)
ggsave("afterfilter.pdf", path = "../results", width = 30, height = 20, units = "cm")
FeatureScatter(seu_obj, features = c("nFeature_RNA", "nCount_RNA", "pMT"))
ggsave("pearplot-after-qc.pdf.pdf", width = 15, height = 6)