比较不同单细胞转录组数据寻找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 条评论
登录 后参与评论

相关文章

来自专栏一心无二用,本人只专注于基础图像算法的实现与优化。

双指数边缘平滑滤波器用于磨皮算法的尝试。

  说起为什么会看到这个东西,那还真的绕一圈。首先在写《Single Image Haze Removal Using Dark Channel Prior》一...

2176
来自专栏计算机视觉战队

简单车辆检测源码

简单车牌检测 main: clc; clear; close all; I=imread('1.png'); figure(1),imshow(I);titl...

3236
来自专栏生信技能树

【直播】我的基因组49:Y染色体的SNV不能用常规流程来找?

在上一次直播中,我们说到了一个不符合我们的认知的问题。就是我的全基因组测序数据里找到的SNV的纯合杂合比例失衡,这着实让我非常纠结。在朋友圈大量求助中,肿瘤所的...

2689
来自专栏数据结构与算法

AGC015 C Nuske vs Phantom Thnook(前缀和)

752
来自专栏生信技能树

【直播】我的基因组80:为什么有些基因的内部测序深度差异如此大

这一讲里,我们依旧根据统计的基因测序的深度进行一下讨论,来看看为什么有些基因的内部测序深度差异如此大? 在前面我们的计算中,s列表示的是基因的每一个坐标的...

3277
来自专栏Y大宽

Cytoscape插件3:Enrichment Map(1)

早期的基因列表解释依赖于选择一系列高得分的基因,然后建立相当主观奇怪的关系。富集分析是一个自动的,基于严格的统计学的方法来分析和解释很大的基因列表,使用的是先验...

772
来自专栏人工智能头条

利用GPU和Caffe训练神经网络

1385
来自专栏新智元

从“London”出发,8步搞定自然语言处理(Python代码)

【新智元导读】自然语言处理是AI的一个子领域,从人们日常沟通所用的非结构化文本信息中提取结构化数据,以便计算机理解。本文用通俗易懂的语言深入浅出的介绍了自然语言...

532
来自专栏LET

glTF(二):PBR

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

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

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

783

扫描关注云+社区