scRNA-seq分析

强推这篇文章,非常详尽,各种图例,眼花缭乱。 Analysis of single cell RNA-seq data

单细胞测序就是获取单个细胞遗传信息的测序技术,能够检出混杂样品测序所无法得到的异质性信息。

安装R包

SC3是一个非常好的用来分析单细胞测序的R包。

source("https://bioconductor.org/biocLite.R")
biocLite("SC3")
biocLite("scater")
biocLite("SingleCellExperiment")

该R包有些参数是依赖perl的,所以还需要有perl的环境。

前数据处理

单细胞测序,如果原始数据只有一个fastq,要先把每个细胞的测序结果分开。 首先使用trim_galore来去接头。 使用trim_galore之前,需要安装cutadapt。

pip install --user upgrade cutadapt

然后接下来是根据barcode来分开每个细胞的测序结果。 可以使用barcode_splitter

# 安装
pip install barcode_splitter
# 使用
barcode_splitter --bcfile barcodes.txt read1.fq.gz read2.fq.gz --idxread 2

当然,如果能确认一开始每个文件就是一个细胞的测序数据,就不用做这一步了。

这里推荐一个公开的不错的练手数据。GSE36552 一共是124个细胞,每个细胞单独为一个fq文件。

然后就是比对还有统计counts数/RPKM。 这里推荐参考上次写的RNA-seq流程,查看Hisat2用法查看featuteCounts用法

得到的结果,就是一个counts矩阵了。

使用SC3进行分析

SC3 package manual

这里是使用RPKM值的分析流程。

library(SingleCellExperiment)
library(SC3)
library(scater)

这里使用的是yan的数据集,作为SC3的示例数据。 这个数据是人胚胎发育不同时期的单细胞测序数据。 主要是两个表格 一个是样本信息

head(ann)
##                 cell_type1
## Oocyte..1.RPKM.     zygote
## Oocyte..2.RPKM.     zygote
## Oocyte..3.RPKM.     zygote
## Zygote..1.RPKM.     zygote
## Zygote..2.RPKM.     zygote
## Zygote..3.RPKM.     zygote

另外一个是RPKM值

yan[1:3, 1:3]
##          Oocyte..1.RPKM. Oocyte..2.RPKM. Oocyte..3.RPKM.
## C9orf152             0.0             0.0             0.0
## RPS11             1219.9          1021.1           931.6
## ELMO2                7.0            12.2             9.3

实际操作时,根据需求选择究竟时用RPKM值还是用counts来做分析。 这里有文章说明了counts和FPKM值(pair ends数据为FPKM;single end数据为RPKM)的差别。 Counts vs. FPKMs in RNA-seq

# 创建SingleCellExperiment对象
sce <- SingleCellExperiment(
    assays = list(
        counts = as.matrix(yan),
        logcounts = log2(as.matrix(yan) + 1)
    ), 
    colData = ann
)

# 创建基因名
rowData(sce)$feature_symbol <- rownames(sce)
# 去重复
sce <- sce[!duplicated(rowData(sce)$feature_symbol), ]

# 想要特别关注的基因
isSpike(sce, "ERCC") <- grepl("ERCC", rowData(sce)$feature_symbol)

# PCA
plotPCA(sce, colour_by = "cell_type1")

PCA

可以看出基本上胚胎发育不同时期的细胞已经有明显差异了。

运行SC3

将数据在ks(簇的数量)范围内探索数据的聚类,这里采取的范围是2:4。

sce <- sc3(sce, ks = 2:4, biology = TRUE)

输出结果到excel,这一步需要perl环境。

sc3_export_results_xls(sce)

和谐矩阵,计算聚类的和谐度(?),当对角线为红色而其他是蓝色时,说明聚类的效果好。

sc3_plot_consensus(
    sce, k = 3, 
    show_pdata ="cell_type1"
)

consensus

基因表达矩阵

sc3_plot_expression(
    sce, k = 3, 
    show_pdata = "cell_type1"
)

expression

差异表达

sc3_plot_de_genes(
    sce, k = 3, 
    show_pdata = "cell_type1"
)

de

Marker基因

sc3_plot_markers(
    sce, k = 3, 
    show_pdata = "cell_type1"
)

markers

SC3这个包的特性很多,而且作者也在论文里讲了和另外一些包的分析结果对比。实际上,看作者前面的那篇终极教案,非常厉害。 争取一个月内学习完吧。

参考文献

Analysis of single cell RNA-seq data

SC3: consensus clustering of single-cell RNA-seq data

Counts vs. FPKMs in RNA-seq

Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells

单细胞测序扫盲:是什么?为什么?怎么做?