Codeplot
国家基因库单细胞训练营 -- 单细胞RNA数据基础分析

单细胞RNA-seq数据基础分析

作者:宋雨默

一、数据处理

1. 文库构建

  下图为DNBelab C系列高通量单细胞RNA文库构建测序流程。首先,制备高质量的单细胞悬液,保证活率在80%以上,细胞大小不超过40um;然后,将制备好的单细胞悬液、用于捕获单细胞mRNA的barcode标记的磁珠和液滴生成油上样到便携式微流控装置中。该装置能生成大量油包水液滴,每个液滴中理想状态下包含一个beads和一个cell,cell膜在液滴中破裂释放出mRNA,被磁珠捕获,同一个cell的mRNA经过反转录会被带上相同的barcode分子标签,实现了单个细胞的标记。对反转录及扩增之后的cDNA进行质控,浓度和片段分布合格的cDNA用于后续的打断建库,形成独有的环状文库,对该文库进行线性扩增,拷贝出具成千上万个重复单元的DNA 纳米球结构DNB,进而进行后续的测序。

2. 下机数据质控

数据结构

  下机数据包括双端cDNA以及Oligo fastq文件。

  DNA数据结构:Read1为10bp barcode + 10bp barcode + 10bp UMI, Read2为100bp cDNA;

  Oligo数据结构:Read1为10bp barcode + 10bp barcode,Read2为10bp UMI + barcode + 10bp barcode。

  过滤处理后的fastq文件如下(以cDNA fastq文件为例):

3. 参考基因组比对、基因组注释及UMI校正

步骤

  • STAR (version 2.7.1a)比对参考基因组;

  • 根据gtf文件注释reads比对基因组的区域;

  • UMI校正:一个细胞中每个基因的UMI的多样性被用来评估基因表达水平,但由于测序或PCR错误可能导致UMI的错误,引入偏差。基于Hamming距离来校正UMI。在默认情况下,来自一个细胞的同一基因的两组UMI,如果Hamming距离等于1,将被视为来自同一转录本。[PISA]

4. Beads calling及beads合并

  采用EmptyDrops方法根据统计的原始beads x genes矩阵进行有效beads选取。[EmptyDrops]

  根据beads的oligo种类及数目计算beads间余弦相似度,将每个bead与其他beads的相似值由高到低排序,每个bead根据和该bead相似度最高的空bead作为阈值来过滤,多beads合并算法阈值以下求连通分量,阈值以上求强连通分量。

5. 细胞基因表达

  采用PISA count模块对细胞、基因进行统计,输出用于下游分析的三个标准矩阵文件,该格式由于存储非零count值,能够有效减小文件的大小。三个文件格式解释如下:

  • 'barcodes.tsv.gz':细胞id文件;

  • 'features.tsv.gz':基因名文件;

  • 'matrix.mtx.gz':细胞、基因umi count文件,前两行为注释行,第三行分别为总基因数、总细胞数以及总非零count数,第四行至最后分别为基因索引信息,细胞索引信息以及对应索引的非零umi count值。


二、基于表达矩阵的基本分析

  Seurat 对象是储存单细胞数据集的数据(如计数矩阵)和分析(如 PCA 或聚类结果)的容器。

In []:

library(dplyr)
library(cowplot)
library(Seurat)
library(ggplot2)
library(psych)
library(qgraph)
library(igraph)
library(Matrix)
library(SeuratWrappers)
library(pryr)

In [2]:

# 设置工作路径
 input_dir <- "/opt/example/scRNA-Example3/Example3/"

# 读取表达矩阵
 data <- Read10X(paste(input_dir,"PBMC-1/FilterMatrix",sep=''),gene.column = 1)

# 读取.h5文件
# data <- Read10X_h5(file)

# 创建object对象
object <- CreateSeuratObject(data)

object

In [3]:

# 统计矩阵行列
dim(object)
, ,
  1. 32356
  2. 2000
,

In [4]:

# 统计细胞表达分布
hist(colSums(object),
     breaks = 150, main = "UMI count per cell",
     xlab = "UMI count per cell")

rm(data)
gc()

, , , , , , , , , ,
A matrix: 2 × 6 of type dbl
used(Mb)gc trigger(Mb)max used(Mb)
Ncells 3352711179.1 5645538301.6 5645538301.6
Vcells10508291 80.229552883225.523997029183.1
,

img

1. Cell QC

  细胞质控是下游分析的前提条件,它决定了分析的准确性。常用的质控指标包括细胞的基因数、UMI数及线粒体比例。

  • 单个细胞中检测到的基因的数目

    • 低质量的细胞以及空胞液滴中一般检测到很少的基因

    • 包含多个细胞的液滴会检测到异常多的基因

  • 单个细胞UMI数目

  • 单个细胞线粒体比例

    • 低质量或死细胞(细胞质mRNA泄漏)常存在线粒体污染

  建议t:

  • 过滤离群值及低质量细胞;

  • 不同样本的过滤参数不同。
    In [5]:

# 计算线粒体比例
object[["percent.mt"]] <- PercentageFeatureSet(object, pattern = "^MT-")

# 过滤前
# 通过小提琴图展示过滤指标
VlnPlot(object, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.2)
# 通过散点图展示数据情况
FeatureScatter(object, feature1 = "nCount_RNA", feature2 = "percent.mt")
FeatureScatter(object, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")



img
img
img
In [6]:

# 过滤后
pbmc <- subset(object, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 5)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.2)

FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

rm(object)
gc()

img
img

, , , , , , , , , ,
A matrix: 2 × 6 of type dbl
used(Mb)gc trigger(Mb)max used(Mb)
Ncells 3595052192.0 5645538301.6 5645538301.6
Vcells15144983115.635543459271.229551672225.5
,

img

2. Normalization

  去除测序深度的影响,是细胞间mRNA数目一致,对counts进行normalization,默认使用"LogNormalize" 的方法,即将每个细胞基因的counts除以该细胞细胞的总counts数,乘10,000,再进行对数转换(对数转换可使数据在不同数量级上可比)。
In [7]:

# normalization.method: "CLR": centered log ratio transformation
#                       "RC": equals to "LogNormalize" without log-transformation

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

# 数据保存在pbmc[["RNA"]]@data
pbmc[["RNA"]]@data[1:4, 1:4]
pbmc[["RNA"]]@counts[1:4,1:4]

# 标准化后细胞表达分布 趋近正态分布
hist(colSums(pbmc$RNA@data),
     breaks = 150,
     main = "UMI count per cell after normalization",
     xlab = "UMI count per cell") 

img

3. Feature selection

  选择细胞与细胞之间具有高度变异性的基因(例如在某些细胞高表达,而其他细胞不表达)进行后续分析,这些基因可以代表细胞与细胞间的主要生物学差异,并用于PCA等下游分析。建议nfeatures为2000-5000。
In [8]:

# 默认使用vst,根据均值、方差关系拟合曲线,对特征值标准化,降序排列取n个基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000,verbose=FALSE)
VariableFeaturePlot(pbmc)

Out [8]:
Warning message:
“Transformation introduced infinite values in continuous x-axis”
Warning message:
“Removed 3649 rows containing missing values (geom_point).”