R语言做生信Single Cell RNA-seq生信分析流程

实验记录3：用R包Seurat进行QC、PCA分析与t-SNE聚

2018-12-22  本文已影响568人  MC学公卫

梗概

  1. 将Cellranger中的filtered_gene_bc_matrices用于分析。
  2. 进行质量控制（QC），以删除异常细胞；
  3. 标准化与归一化，消除技术噪音与批次效应；
  4. 主成分分析（PCA）与挑选
  5. t-SNE聚类

参考网站：https://satijalab.org/seurat/pbmc3k_tutorial.html

上次结果：

经过Cellranger的数据整理之后，得到：

Seurat是一种R包，设计用于QC，分析和探索单细胞RNA-seq数据。 Seurat旨在使用户能够从单细胞转录组测量中识别和解释异质性来源，并整合不同类型的单细胞数据。

由于Seurat包在90.9服务器中无法安装，故用90.36进行分析

登录xiongjh@202.116.90.36
将数据传输过去

scp /data/zhengll/project/HCA/download/cellranger/HCATisStabAug177078016/outs/filtered_gene_bc_matrices zhengll@202.116.90.9 /data/zll/project/deepBase3/HCA
#检查数据是否传输成功
cd /data1/zll/project/deepBase3/HCA/
ls
filtered_gene_bc_matrices
cd filtered_gene_bc_matrices/
ls
GRCh38
cd GRCh38/
ls
barcodes.tsv  genes.tsv  matrix.mtx
#传输成功

运行R，并且加载这两个包

library(Seurat)
library(dplyr)

读取数据

spleen.data <- Read10X(data.dir = '/data1/zll/project/deepBase3/HCA/filtered_gene_bc_matrices/GRCh38/')

dim(spleen.data)
[1] 33694  1960

原始数据的基因数为33694，细胞数为1960.

比较普通与疏松矩阵的内存使用：

> dense.size <- object.size(x = as.matrix(x = spleen.data))
> dense.size
530488272 bytes

#转化为疏松矩阵，查看大小
> sparse.size <- object.size(x = spleen.data)
> sparse.size
45955656 bytes

> dense.size/sparse.size
11.5 bytes

初始化Seurat对象：
命令CreateSeuratObject
输入数据spleen.data
留下所有在>=3个细胞中表达的基因min.cells = 3；
留下所有检测到>=200个基因的细胞min.genes = 200。
(为了除去一些)

spleen <- CreateSeuratObject(raw.data = spleen.data, min.cells = 3, min.genes = 200, project = "10X_spleen")

spleen
An object of class seurat in project 10X_spleen 
15655 genes across 1959 samples.

剩下15655 基因和 1959 个细胞

质量控制

以下步骤包括Seurat中scRNA-seq数据的标准预处理工作流程。这些代表了Seurat对象的创建，基于QC指标的细胞选择和过滤，数据标准化和缩放，以及高度可变基因的检测。

mito.genes <- grep(pattern = "^MT-", x = rownames(x = spleen@data), value = TRUE)
percent.mito <- Matrix::colSums(spleen@raw.data[mito.genes, ])/Matrix::colSums(spleen@raw.data)
spleen <- AddMetaData(object = spleen, metadata = percent.mito, col.name = "percent.mito")
VlnPlot(object = spleen, features.plot = c("nGene", "nUMI", "percent.mito"), nCol = 3)
VlnPlot_of_spleen.png 
> par(mfrow = c(1, 2))
> GenePlot(object = spleen, gene1 = "nUMI", gene2 = "percent.mito")
> GenePlot(object = spleen, gene1 = "nUMI", gene2 = "nGene")
GenePlot_of_spleen.png

过滤细胞，根据上面的两幅图，去除异常值，这里选择基因数从300-5000，线粒体基因占比小于0.1的细胞。（主要看小提琴图1和图3）

spleen <- FilterCells(spleen, subset.names = c("nGene", "percent.mito"), low.thresholds = c(300, -Inf), high.thresholds = c(5000,0.10))

查看过滤掉剩下多少细胞：

spleen
An object of class seurat in project 10X_spleen 
 15655 genes across 1940 samples.

剩下15655个基因，1940个细胞。

数据标准化

加个log

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

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

spleen <- FindVariableGenes(object = spleen, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)

Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
TEXT_SHOW_BACKTRACE environmental variable.
> length(x=spleen@var.genes)
[1] 1829
高度变异基因.png

官网给出的筛选差异基因中的一些问题：

How do I select my set of variable genes?
FindVariableGenes() plots dispersion (a normalized measure of to cell-to-cell variation) as a function of average expression for each gene. Our goal is to identify a set of high-variance genes, and we recommend setting the cutoff parameters in this function by visually evaluating this plot to define outliers. However, particularly when using UMI data, we often obtain similar results including much larger gene sets, including the entire transcriptome.

缩放数据并删除不需要的变体来源

您的单细胞数据集可能包含“不感兴趣”的变异来源。这不仅包括技术噪音，还包括批次效应，甚至包括生物变异来源（细胞周期阶段）。正如(Buettner, et al NBT，2015)中所建议的那样，从分析中回归这些信号可以改善下游维数减少和聚类。为了减轻这些信号的影响，Seurat构建线性模型以基于用户定义的变量预测基因表达。这些模型的缩放得分残差存储在Scale.data槽中，用于降维和聚类。

我们可以消除由批次（如果适用）驱动的基因表达中的细胞 - 细胞变异，细胞比对率（由Drop-seq数据的Drop-seq工具提供），检测到的分子数量和线粒体基因表达。对于循环细胞，我们还可以学习“细胞周期”评分（参见此处的示例）并对其进行回归。在这个有丝分裂后血细胞的简单例子中，我们回归了每个细胞检测到的分子数量以及线粒体基因含量百分比。

spleen <-ScaleData(spleen, vars.to.regress = c("nUMI","percent.mito"))

Regressing out: nUMI, percent.mito
  |=========================================================================================| 100%
Time Elapsed:  18.0711550712585 secs
Scaling data matrix
  |=========================================================================================| 100%

PCA分析

主成分分析是什么？

主成分分析，是考察多个变量间相关性一种多元统计方法，研究如何通过少数几个主成分来揭示多个变量间的内部结构，即从原始变量中导出少数几个主成分，使它们尽可能多地保留原始变量的信息，且彼此间互不相关.通常数学上的处理就是将原来P个指标作线性组合，作为新的综合指标。

将数据集降维，利用低阶的变量去反应整体的结果。

spleen <- RunPCA(spleen, pc.genes = spleen@var.genes, do.print = TRUE, pcs.print = 1:5, genes.print = 5)

[1] "PC1"
[1] "CD69"  "CD79A" "TRAC"  "CD3D"  "MS4A1"
[1] ""
[1] "FCN1"          "LYZ"           "SERPINA1"      "CSTA"          "RP11-1143G9.4"
[1] ""
[1] ""
[1] "PC2"
[1] "CD79A"    "MS4A1"    "VPREB3"   "CD79B"    "HLA-DQB1"
[1] ""
[1] "NKG7" "CST7" "GZMA" "CD7"  "CCL5"
[1] ""
[1] ""
[1] "PC3"
[1] "TRDC"  "KLRF1" "MS4A1" "CD79B" "IRF8" 
[1] ""
[1] "IL7R" "TRAC" "CD3D" "CD2"  "CD3G"
[1] ""
[1] ""
[1] "PC4"
[1] "GIMAP7" "GZMB"   "FGFBP2" "SPON2"  "PRF1"  
[1] ""
[1] "BAG3"    "HSPD1"   "FKBP4"   "DNAJA1"  "ZFAND2A"
[1] ""
[1] ""
[1] "PC5"
[1] "UBE2C" "TYMS"  "MKI67" "TOP2A" "AURKB"
[1] ""
[1] "FCGR3A" "FGFBP2" "SPON2"  "GNLY"   "GZMB"  
[1] ""
[1] ""

PCElbowPlot(spleen)
碎石图.jpeg

选择了前10个PC成分

spleen <- FindClusters(spleen, reduction.type = "pca", dims.use = 1:10, resolution = 0.6, print.output = 0, save.SNN = TRUE)
PrintFindClustersParams(spleen)

Parameters used in latest FindClusters calculation run on: 2018-10-01 21:59:55
=============================================================================
Resolution: 0.6
-----------------------------------------------------------------------------
Modularity Function    Algorithm         n.start         n.iter
     1                   1                 100             10
-----------------------------------------------------------------------------
Reduction used          k.param          prune.SNN
     pca                 30                0.0667
-----------------------------------------------------------------------------
Dims used in calculation
=============================================================================
1 2 3 4 5 6 7 8 9 10

细胞聚类

spleen <- RunTSNE(spleen, dims.use = 1:10, do.fast= TRUE)
TSNEPlot(spleen)
TSNE.jpeg 
> saveRDS(spleen, file = "/Users/shinianyike/Desktop/zll/Seurat/spleen_results/spleen_1.rds")

保存了变量

一些补充：
1、过滤低质量细胞
在 scRNA-seq 分析中，有些细胞质量比较低，比如细胞处于凋亡状态,细胞中 RNA 发生降解等,这些细胞的存在会影响分析，因此我们第一步需要对细胞进行过滤。主要可分为三类:
①利用细胞检测到的基因数或者是 reads 比对率来判断技术噪音,但不管是基因检测数目还是比对率都跟实验方法有很大相关性。 如果比对率太低,表明 RNA 可能发生了降解,或者文库有污染或者细胞裂解不完全。
②如果实验中加入了 spike-ins,可以通过计算比对到内源性 RNA 和外源性 RNA(spike-ins)的 reads 比例来过滤低质量细胞,比值偏低表明细胞中的 RNA 数量较低,细胞可丢弃。但是也需要注意 其实当细胞状态不一样,比如处于不同细胞周期时,细胞的 RNA 数 量是具有很大差异的。不过我们依然认为在一大群细胞中,spike-ins 比例特别高的细胞在很大概率上应该被排除在外。这里软件 SinQC (Single-cell RNA-seq Quality Control)就可以根据比对率和检测到 的基因数来过滤细胞(图 1)（本实验没有）。
③根据整体的基因表达谱来定义技术噪音。比如对细胞进行聚类分析，PCA 分析等，将 outlier 细胞删除，或者细胞表达中位值低于某一设定阈值时将该细胞过滤掉。当然这种方法也存在误删具有真正生物学差异的细胞,因此在删除细胞时需要小心，可与上述另外两种方法连用。

上一篇下一篇

猜你喜欢

热点阅读