比较不同单细胞转录组数据寻找features方法

挑选到的跟feature相关的基因集，有点类似于在某些组间差异表达的基因集，都需要后续功能注释。

背景介绍

单细胞转录组测序的确可以一次性对所有细胞都检测到上千个基因的表达，但是，大多数情况下，只有其中的少部分基因是有生物学意义的，比如可以区分不同的细胞类型，或者分化发育相关的基因，或者细胞应对外界刺激的。

而且大多数基因之所以在不同的细胞里面表达有差异，其实是技术限制，背景噪音。这些技术限制，包括批次效应，都会阻碍我们发现那些真正的有生物学意义的基因。

所以做 feature selection 分析来去除那些技术噪音相关基因，可以显著的提高信噪比，降低后续分析的复杂度。

包的安装

所以需要安装并且加载一些包,安装代码如下；

install.packages('ROCR')
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("M3Drop") 

install.packages("devtools")
library("devtools")
install_github("BPSC","nghiavtr")
install_github("hemberg-lab/scRNA.seq.funcs") 

加载代码如下：

library(scRNA.seq.funcs)
library(matrixStats)
library(M3Drop)
library(RColorBrewer)
set.seed(1)

加载测试数据

这里选取的是Usoskin et al 文章单细胞测序文章的数据，包含4种细胞类型：

• NP = non-peptidergic nociceptors

• PEP = peptidergic nociceptors

• NF = neurofilament containing

• TH = tyrosine hydroxylase containing neurons.

对应着25334个基因在 622 个细胞里面的表达量

# http://www.biotrainee.com/jmzeng/scRNA/usoskin1.rds
usoskin1 <- readRDS("../usoskin/usoskin1.rds")
dim(usoskin1)

## [1] 25334   622

table(colnames(usoskin1))

## 
##  NF  NP PEP  TH 
## 139 169  81 233

用M3Drop对表达矩阵进行一站式过滤

uso_list <- M3Drop::M3DropCleanData(
    usoskin1,
    labels = colnames(usoskin1),
    min_detected_genes = 2000,
    is.counts = TRUE
)
expr_matrix <- uso_list$data # Normalized & filtered expression matrix
dim(expr_matrix)

## [1] 15708   532

celltype_labs <- uso_list$labels # filtered cell-type labels
cell_colors <- brewer.pal(max(3,length(unique(celltype_labs))), "Set3")

这个M3Drop的M3DropCleanData函数会自动过滤那些表达基因数量很少的细胞，过滤低表达基因，然后把reads counts的表达矩阵转换为 counts per million (CPM) 。

过滤之后，只剩下 15708个基因在532个细胞 的表达了。

寻找highly variable genes (HVG)

那些在样本群体里面表达量变异比较大的基因可能是真正的生物学现象，也有可能是技术误差，而且变异程度总是跟基因的表达量成正相关。如下图所示：

plot(rowMeans(expr_matrix),rowVars(expr_matrix),log='xy')

Brennecke et al.提出了算法来矫正这一影响，这个方法也被包装成了Brennecke_getVariableGenes(counts, spikes) 函数，但是这个数据并没有ERCC spike-in，所以直接对整个表达矩阵处理即可。

Brennecke_HVG <- M3Drop::BrenneckeGetVariableGenes(
    expr_matrix,
    fdr = 0.01,
    minBiolDisp = 0.5
)

探究 High Dropout Genes

另外一个寻找HVGs是查看它们是否有非常多的0表达量情况，这种多0表达的情况叫做dropout rate，通常单细胞转录组表达矩阵里面过半数的基因都是0表达的。因为单细胞里面的mRNA很多无法被反转录，这种情况可以用Michaelis-Menten等式来模拟，如下图所示：

K = 49
S_sim = 10^seq(from=-3, to=4, by=0.05) # range of expression values
MM = 1-S_sim/(K+S_sim)
plot(S_sim, MM, type="l", lwd=3, xlab="Expression", ylab="Dropout Rate", xlim=c(1,1000))
S1 = 10; P1 = 1-S1/(K+S1) # Expression & dropouts for cells in condition 1
S2 = 750; P2 = 1-S2/(K+S2) # Expression & dropouts for cells in condition 2
points(c(S1,S2),c(P1,P2), pch=16, col="grey85", cex=3)
mix = 0.5; # proportion of cells in condition 1
points(S1*mix+S2*(1-mix), P1*mix+P2*(1-mix), pch=16, col="grey35", cex=3)

用来M3Drop包的M3DropFeatureSelection函数来挑选那些显著偏离了Michaelis-Menten曲线的基因，这里的阈值取1% FDR.

但是这个函数M3DropFeatureSelection依赖于正确的M3Drop包版本，下面就不运行了。

M3Drop_genes <- M3Drop::M3DropFeatureSelection(
    expr_matrix,
    mt_method = "fdr",
    mt_threshold = 0.01
)
title(main = "Usoskin")
M3Drop_genes <- M3Drop_genes$Gene

Depth-Adjusted Negative Binomial (DANB)

下面这个 NBumiConvertToInteger 也依赖于正确的M3Drop包版本，下面就不运行了。

usoskin_int <- NBumiConvertToInteger(usoskin1)
DANB_fit <- NBumiFitModel(usoskin_int) # DANB is fit to the raw count matrix
# Perform DANB feature selection
DropFS <- NBumiFeatureSelectionCombinedDrop(DANB_fit)
DANB_genes <- names(DropFS[1:1500])

基因表达相关性

这个表达矩阵很大，所以计算所有基因之间的相关性耗时很长，为节约时间，也不运行了。

cor_mat <- cor(t(expr_matrix), method="spearman") #Gene-gene correlations
diag(cor_mat) <- rep(0, times=nrow(expr_matrix))
score <- apply(cor_mat, 1, function(x) {max(abs(x))}) #Correlation of highest magnitude
names(score) <- rownames(expr_matrix);
score <- score[order(-score)]
Cor_genes = names(score[1:1500])

PCA挑选

PCA 速度还行，挑选第1，2主成分的前1500个基因

pca <- prcomp(log(expr_matrix+1)/log(2)); 
# PCA is typically performed on log-transformed expression data

plot(pca$rotation[,1], pca$rotation[,2], pch=16, col=cell_colors[as.factor(celltype_labs)]) # plot projection

score <- rowSums(abs(pca$x[,c(1,2)])) 
# calculate loadings for components 1 and 2
names(score) <- rownames(expr_matrix)
score <- score[order(-score)]
PCA_genes = names(score[1:1500])

检查挑选的基因集的效果

热图+聚类可以看看基因是否在各个细胞类型差异表达，并且把细胞类型比较好的分开。

这个热图非常耗时，如无必要，请不要运行这个代码

M3Drop::M3DropExpressionHeatmap(
    PCA_genes , ## 或者 M3Drop_genes 等其它方法挑到的基因
    uso_list$data,
    cell_labels = uso_list$labels
)

挑选的基因集跟DEseq得到的差异基因列表看交集

载入用DEseq得到的差异基因列，跟前面得到的M3Drop_genes比较一下。

# Load DE genes
# http://www.biotrainee.com/jmzeng/scRNA/DESeq_table.rds
DESeq_table <- readRDS("../usoskin/DESeq_table.rds")
DE_genes = unique(DESeq_table$Gene)

# Calculate precision
# sum(M3Drop_genes %in% DE_genes)/length(M3Drop_genes)
sum(PCA_genes %in% DE_genes)/length(PCA_genes)
## [1] 0.382

文中提到的测试数据：

http://www.biotrainee.com/jmzeng/scRNA/DESeq_table.rds

http://www.biotrainee.com/jmzeng/scRNA/pollen.rds

http://www.biotrainee.com/jmzeng/scRNA/tung_umi.rds

http://www.biotrainee.com/jmzeng/scRNA/usoskin1.rds

前面都忘记说了，后台回复单细胞，是有惊喜的！