了解scRNAseq包内置的单细胞转录组数据

首先需要下载并且安装这个包：

source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("scRNAseq")
#biocLite("scone")


这个包内置的是Pollen et al. 2014 数据集，而且只挑选 65个人类细胞，分成3类，分别是 pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ，还有 “GW16” and “GW21” ，“GW21+3” 这种孕期细胞。

大小是50.6 MB，下载需要一点点数据，先安装加载它们。

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm)
# Set assay to RSEM estimated counts
assay(fluidigm) = assays(fluidigm)$rsem_counts
# List all qc fields (accessible via colData())
metadata(fluidigm)$which_qc
##  [1] "NREADS"                       "NALIGNED"                    
##  [3] "RALIGN"                       "TOTAL_DUP"                  
##  [5] "PRIMER"                       "INSERT_SZ"                  
##  [7] "INSERT_SZ_STD"                "COMPLEXITY"                  
##  [9] "NDUPR"                        "PCT_RIBOSOMAL_BASES"        
## [11] "PCT_CODING_BASES"             "PCT_UTR_BASES"              
## [13] "PCT_INTRONIC_BASES"           "PCT_INTERGENIC_BASES"        
## [15] "PCT_MRNA_BASES"               "MEDIAN_CV_COVERAGE"          
## [17] "MEDIAN_5PRIME_BIAS"           "MEDIAN_3PRIME_BIAS"          
## [19] "MEDIAN_5PRIME_TO_3PRIME_BIAS"
## 可以看到，默认可以做QC的项目非常多。
## 值得注意的是 RALIGN 代表比对率， NREADS 代表文库大小。
# Joint distribution of "biological condition"" and "coverage type""
table(colData(fluidigm)$Coverage_Type,
     colData(fluidigm)$Biological_Condition)
##      
##        GW16 GW21 GW21+3 NPC
##   High   26    8     16  15
##   Low    26    8     16  15


这里面的表达矩阵是由 RSEM (Li and Dewey 2011) 软件根据 hg38 RefSeq transcriptome 得到的，总是130个文库，每个细胞测了两次，测序深度不一样。

可以用scran包的cyclone函数进行单细胞转录组的细胞周期状态推断

代码如下，需要注意是针对人类的测序数据，然后还有一些ID的转换。

library(scRNAseq)
library(RColorBrewer)
library(scran)
all.counts=assays(fluidigm)$rsem_counts
dim(all.counts)
head(rownames(all.counts))
sce <- SingleCellExperiment(list(counts=all.counts))
library(org.Hs.eg.db)
ensembl <- mapIds(org.Hs.eg.db, keys=rownames(sce), keytype="SYMBOL", column="ENSEMBL")
mm.pairs <- readRDS(system.file("exdata", "human_cycle_markers.rds", package="scran"))
assigned <- cyclone(sce, pairs=mm.pairs, gene.names=ensembl)
head(assigned$scores)
plot(assigned$score$G1, assigned$score$G2M,
    xlab="G1 score", ylab="G2/M score", pch=16)

table(assigned$phases)


结果是：

G1 G2M   S
12  38  80

因为是真实数据，所以细胞周期的推断比较明显：

细胞周期比较分明：

当然，这里还没有用5个阶段的细胞周期推断：

使用scran包的cyclone函数进行单细胞转录组的细胞周期状态推断

后期我会出补充教程，欢迎大家继续围观。

点击可以加入单细胞数据处理学习交流小组