比较不同单细胞转录组数据寻找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")

这个M3DropM3DropCleanData函数会自动过滤那些表达基因数量很少的细胞,过滤低表达基因,然后把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

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-01-22

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏SDNLAB

5G革命的技术,一个都不能少

第五代移动网络简称5G是产业界即将实现的移动技术革命,是LTE-A网络的深层演进技术。5G网络中的关键技术包括MIMO、OFDM、SC-FDMA等。 超密集微型...

42212
来自专栏新智元

【20张图玩转机器学习】深度学习、神经网络和大数据信息梳理(下载)

【新智元导读】ChatbotLife 的创始人兼编辑 Stefan Kojouharov 收集并整理了一系列 AI 相关的信息图示,为了便于使用,还附带了注释和...

3965
来自专栏程序你好

程序员算法基础——贪心算法

1182
来自专栏机器之心

入门 | 玩转词向量:用fastText预训练向量做个智能小程序

5779
来自专栏大数据文摘

妙解谷歌压箱底面试题:如何正确的从楼上抛鸡蛋

1533
来自专栏AI科技大本营的专栏

计算机如何理解我们的语言?NLP is fun!

【导读】我们从日常每天都会用到的推荐系统到现在研究火热的开放性聊天、对话机器人,越来越多的产品与应用的背后都需要自然语言处理(NLP)和知识图谱的技术。也有越来...

913
来自专栏数据魔术师

干货 | 变邻域搜索算法(VNS)求解TSP(附C++详细代码及注释)

上次变邻域搜索的推文发出来以后,看过的小伙伴纷纷叫好。小编大受鼓舞,连夜赶工,总算是完成了手头上的一份关于变邻域搜索算法解TSP问题的代码。今天,就在此...

9782
来自专栏每日一篇技术文章

opengL ES _ 入门_05

ID是漫反射的强度,Ii是光的入射光的强度,和KD的漫反射,是对粗糙松散耦合对象材料。松散的意思是,在许多现实世界的材料,实际表面可能有点抛光,但半透明的,而层...

1303
来自专栏Python小屋

Python计算电场中两点间的电势差

根据组合数定义,需要计算3个数的阶乘,在很多编程语言中都很难直接使用整型变量表示大数的阶乘结果,虽然Python并不存在这个问题,但是计算大数的阶乘仍需要相当多...

751
来自专栏生信技能树

如何通过Google来使用ggplot2可视化

今天是大年初二,这篇文章我只想传达一点: 没有什么菜鸟级别的生物信息学数据处理是不能通过Google得到解决方案的,如果有,请换个关键词继续Google! 第一...

3088

扫码关注云+社区