前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >比较不同单细胞转录组数据寻找features方法

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

作者头像
生信技能树
发布2018-03-09 10:45:43
3K0
发布2018-03-09 10:45:43
举报
文章被收录于专栏:生信技能树

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

背景介绍

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

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

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

包的安装

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

代码语言:javascript
复制
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") 

加载代码如下:

代码语言:javascript
复制
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 个细胞里面的表达量

代码语言:javascript
复制
# http://www.biotrainee.com/jmzeng/scRNA/usoskin1.rds
usoskin1 <- readRDS("../usoskin/usoskin1.rds")
dim(usoskin1)
代码语言:javascript
复制
## [1] 25334   622
代码语言:javascript
复制
table(colnames(usoskin1))
代码语言:javascript
复制
## 
##  NF  NP PEP  TH 
## 139 169  81 233

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

代码语言:javascript
复制
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)
代码语言:javascript
复制
## [1] 15708   532
代码语言:javascript
复制
celltype_labs <- uso_list$labels # filtered cell-type labels
cell_colors <- brewer.pal(max(3,length(unique(celltype_labs))), "Set3")

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

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

寻找highly variable genes (HVG)

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

代码语言:javascript
复制
plot(rowMeans(expr_matrix),rowVars(expr_matrix),log='xy')

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

代码语言:javascript
复制
Brennecke_HVG <- M3Drop::BrenneckeGetVariableGenes(
    expr_matrix,
    fdr = 0.01,
    minBiolDisp = 0.5
)

探究 High Dropout Genes

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

代码语言:javascript
复制
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包版本,下面就不运行了。

代码语言:javascript
复制
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包版本,下面就不运行了。

代码语言:javascript
复制
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])

基因表达相关性

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

代码语言:javascript
复制
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个基因

代码语言:javascript
复制
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
代码语言:javascript
复制
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])

检查挑选的基因集的效果

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

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

代码语言:javascript
复制
M3Drop::M3DropExpressionHeatmap(
    PCA_genes , ## 或者 M3Drop_genes 等其它方法挑到的基因
    uso_list$data,
    cell_labels = uso_list$labels
)

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

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

代码语言:javascript
复制
# 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

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2018-01-22,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 背景介绍
  • 包的安装
  • 加载测试数据
    • 用M3Drop对表达矩阵进行一站式过滤
      • 寻找highly variable genes (HVG)
        • 探究 High Dropout Genes
          • Depth-Adjusted Negative Binomial (DANB)
            • 基因表达相关性
              • PCA挑选
              • 检查挑选的基因集的效果
              • 挑选的基因集跟DEseq得到的差异基因列表看交集
              领券
              问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档