比较不同的对单细胞转录组数据normalization方法

使用CPM去除文库大小影响

之所以需要normalization,就是因为测序的各个细胞样品的总量不一样,所以测序数据量不一样,就是文库大小不同,这个因素是肯定需要去除。最简单的就是counts per million (CPM),所有样本的所有基因的表达量都乘以各自的文库reads总数再除以一百万即可。(一般miRNA-seq数据结果喜欢用这个) 代码如下:

calc_cpm <-
function (expr_mat, spikes = NULL) 
{
    norm_factor <- colSums(expr_mat[-spikes, ])
    return(t(t(expr_mat)/norm_factor)) * 10^6
}

但是CPM方法有一个很严重的缺陷,那些高表达并且在细胞群体表达差异很大的基因会严重影响那些低表达基因。

RPKM, FPKM and TPM去除基因或者转录本长度影响

最常见的是下面3个:

  • RPKM - Reads Per Kilobase Million (for single-end sequencing)
  • FPKM - Fragments Per Kilobase Million (same as RPKM but for paired-end sequencing, makes sure that paired ends mapped to the same fragment are not counted twice)
  • TPM - Transcripts Per Kilobase Million (same as RPKM, but the order of normalizations is reversed - length first and sequencing depth second)

这些normalization方法并不适合单细胞转录组测序数据,因为有一些scRNA-seq建库方法具有3端偏好性,一般是没办法测全长转录本的,所以转录本的长度跟表达量不是完全的成比例。

对于这样的数据,需要重新转换成 reads counts 才能做下游分析。

适用于bulk RNA-seq的normalization方法

比较流行的有:

  • DESeq的size factor (SF)
  • relative log expression(RLE)
  • upperquartile (UQ)
  • weighted trimmed mean of M-values(TMM)

这些适用于 bulk RNA-seq data 的normalization方法可能并不适合 single-cell RNA-seq data ,因为它们的基本假设是有问题的。

特意为single-cell RNA-seq data 开发的normalization方法

  • LSF (Lun Sum Factors)
  • scran package implements a variant on CPM specialized for single-cell data

而scater包把这些normalization方法都包装到了normaliseExprs函数里面,可以直接调用。

并且通过plotPCA函数来可视化这些normalization的好坏。

工作环境

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

## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R") 
biocLite("scater") 
biocLite("scran") 
install.packages("devtools")
library("devtools") 
install_github("hemberg-lab/scRNA.seq.funcs") 

加载代码如下:

library(scRNA.seq.funcs)
library(scater)
library(scran)
set.seed(1234567)

加载测试数据

这里选取的是芝加哥大学Yoav Gilad lab实验的Tung et al 2017的单细胞测序文章的数据

options(stringsAsFactors = FALSE)
set.seed(1234567)
## 这里直接读取过滤好的数据,是一个SCESet对象,适用于scater包的
## http://www.biotrainee.com/jmzeng/scRNA/tung_umi.rds
umi <- readRDS("tung_umi.rds")

## 如果没有这个rds对象,就自己把read counts的表达矩阵读进去,变成这个适用于scater包的SCESet对象,代码如下;
if(F){
      # 这个文件是表达矩阵,包括线粒体基因和 ERCC spike-ins 的表达量,可以用来做质控
    molecules <- read.table("tung/molecules.txt", sep = "\t")
    ## 这个文件是表达矩阵涉及到的所有样本的描述信息,包括样本来源于哪个细胞,以及哪个批次。
    anno <- read.table("tung/annotation.txt", sep = "\t", header = TRUE)
    pheno_data <- new("AnnotatedDataFrame", anno)
    rownames(pheno_data) <- pheno_data$sample_id
    dat <- scater::newSCESet(
      countData = molecules,
      phenoData = pheno_data
    )
  ## 这个代码过时了,请注意!!!
    set_exprs(dat, "log2_counts") <- log2(counts(dat) + 1)

}

umi.qc <- umi[fData(umi)$use, pData(umi)$use] 
## counts(umi) 和  exprs(umi) 这里是不一样的。
## 前面的过滤信息,这里直接用就好了。
endog_genes <- !fData(umi.qc)$is_feature_control
dim(exprs( umi.qc[endog_genes, ]))
## [1] 13997   654
## 可以看到是过滤后的654个单细胞的13997个基因的表达矩阵。
umi.qc
## SCESet (storageMode: lockedEnvironment)
## assayData: 14063 features, 654 samples 
##   element names: counts, exprs, log2_counts 
## protocolData: none
## phenoData
##   rowNames: NA19098.r1.A01 NA19098.r1.A02 ... NA19239.r3.H12 (654
##     total)
##   varLabels: individual replicate ... outlier (47 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: ENSG00000237683 ENSG00000187634 ... ERCC-00171
##     (14063 total)
##   fvarLabels: mean_exprs exprs_rank ... use (13 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:

实践

Raw

先看看原始的表达值的分布情况,这里本来应该是对每一个样本画boxplot的,但是这里的样本数量太多了,这样的可视化效果很差, 就用PCA的方式,看看这表达矩阵是否可以把样本区分开,只有那些区分度非常好的normalization方法才是最优的。

不过scater包提供了一个plotRLE函数,可以画出类似于样本boxplot的效果。

plotPCA(
    umi.qc[endog_genes, ],
    colour_by = "batch",
    size_by = "total_features",
    shape_by = "individual",
    exprs_values = "log2_counts"
)

CPM

scater默认对表达矩阵做了cpm转换,所以可以直接提取里面的信息

plotPCA(
    umi.qc[endog_genes, ],
    colour_by = "batch",
    size_by = "total_features",
    shape_by = "individual",
    exprs_values = "exprs"
)

还可以看看CPM和原始的log转换的表达矩阵的区别

plotRLE(
    umi.qc[endog_genes, ], 
    exprs_mats = list(Raw = "log2_counts", CPM = "exprs"),
    exprs_logged = c(TRUE, TRUE),
    colour_by = "batch"
)

TMM

需要用函数 normaliseExprs 来对SCESet对象里面的表达矩阵做TMM转换,

umi.qc <- normaliseExprs(
    umi.qc,
    method = "TMM",
    feature_set = endog_genes
)
plotPCA(
    umi.qc[endog_genes, ],
    colour_by = "batch",
    size_by = "total_features",
    shape_by = "individual",
    exprs_values = "norm_exprs"
)

这次的转换会以normexprs的属性存储在里面,同时增加了一个normcpm属性。

plotRLE(
    umi.qc[endog_genes, ], 
    exprs_mats = list(Raw = "log2_counts", TMM = "norm_exprs"),
    exprs_logged = c(TRUE, TRUE),
    colour_by = "batch"
)

观察变化规律

到这里为止,表达矩阵已经有了 counts, exprs, log2counts, normcpm, norm_exprs 这些形式。

## 最开始读入是 基于read counts的表达矩阵
counts(umi.qc)[1:10,1:3]
##                 NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683              0              0              0
## ENSG00000187634              0              0              0
## ENSG00000188976              3              6              1
## ENSG00000187961              0              0              0
## ENSG00000187608              1              7              0
## ENSG00000188157              2              3              2
## ENSG00000131591              0              0              0
## ENSG00000078808              3              2              0
## ENSG00000176022              0              0              1
## ENSG00000160087              3              3              0
## 默认的CPM转换后的矩阵
exprs(umi.qc)[1:10,1:3]
##                 NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683       0.000000       0.000000       0.000000
## ENSG00000187634       0.000000       0.000000       0.000000
## ENSG00000188976       2.167284       3.008380       1.400312
## ENSG00000187961       0.000000       0.000000       0.000000
## ENSG00000187608       1.113650       3.204929       0.000000
## ENSG00000188157       1.734589       2.177376       2.097332
## ENSG00000131591       0.000000       0.000000       0.000000
## ENSG00000078808       2.167284       1.743673       0.000000
## ENSG00000176022       0.000000       0.000000       1.400312
## ENSG00000160087       2.167284       2.177376       0.000000
## 通过set_exprs进行简单的对数转换后的表达矩阵。
log2(counts(umi.qc) + 1)[1:10,1:3]
##                 NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683       0.000000       0.000000       0.000000
## ENSG00000187634       0.000000       0.000000       0.000000
## ENSG00000188976       2.000000       2.807355       1.000000
## ENSG00000187961       0.000000       0.000000       0.000000
## ENSG00000187608       1.000000       3.000000       0.000000
## ENSG00000188157       1.584963       2.000000       1.584963
## ENSG00000131591       0.000000       0.000000       0.000000
## ENSG00000078808       2.000000       1.584963       0.000000
## ENSG00000176022       0.000000       0.000000       1.000000
## ENSG00000160087       2.000000       2.000000       0.000000
## 通过normaliseExprs函数指定 TMM 转换
norm_exprs(umi.qc)[1:10,1:3]
##                 NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683       0.000000       0.000000       0.000000
## ENSG00000187634       0.000000       0.000000       0.000000
## ENSG00000188976       2.167284       3.008380       1.400312
## ENSG00000187961       0.000000       0.000000       0.000000
## ENSG00000187608       1.113650       3.204929       0.000000
## ENSG00000188157       1.734589       2.177376       2.097332
## ENSG00000131591       0.000000       0.000000       0.000000
## ENSG00000078808       2.167284       1.743673       0.000000
## ENSG00000176022       0.000000       0.000000       1.400312
## ENSG00000160087       2.167284       2.177376       0.000000
## 对TMM转换后,再进行cpm转换的表达矩阵。
norm_cpm(umi.qc)[1:10,1:3]
##                 NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683        0.00000        0.00000        0.00000
## ENSG00000187634        0.00000        0.00000        0.00000
## ENSG00000188976       47.28989       95.43385       22.20531
## ENSG00000187961        0.00000        0.00000        0.00000
## ENSG00000187608       15.76330      111.33949        0.00000
## ENSG00000188157       31.52660       47.71693       44.41062
## ENSG00000131591        0.00000        0.00000        0.00000
## ENSG00000078808       47.28989       31.81128        0.00000
## ENSG00000176022        0.00000        0.00000       22.20531
## ENSG00000160087       47.28989       47.71693        0.00000
# PS: 记住,这个时候是没有norm_counts(umi.qc) 函数的。

scran

这个scran package implements a variant on CPM specialized for single-cell data,所以需要特殊的代码

qclust <- quickCluster(umi.qc, min.size = 30)
umi.qc <- computeSumFactors(umi.qc, sizes = 15, clusters = qclust)
umi.qc <- normalize(umi.qc)
plotPCA(
    umi.qc[endog_genes, ],
    colour_by = "batch",
    size_by = "total_features",
    shape_by = "individual",
    exprs_values = "exprs"
)

也可以比较它相当于最粗糙的对数转换,效果好在哪里。

plotRLE(
    umi.qc[endog_genes, ], 
    exprs_mats = list(Raw = "log2_counts", scran = "exprs"),
    exprs_logged = c(TRUE, TRUE),
    colour_by = "batch"
)

Size-factor (RLE)

这个normalization方法最初是DEseq包提出来的。

umi.qc <- normaliseExprs(
    umi.qc,
    method = "RLE", 
    feature_set = endog_genes
)
plotPCA(
    umi.qc[endog_genes, ],
    colour_by = "batch",
    size_by = "total_features",
    shape_by = "individual",
    exprs_values = "norm_exprs"
)

Upperquantile

umi.qc <- normaliseExprs(
    umi.qc,
    method = "upperquartile", 
    feature_set = endog_genes,
    p = 0.99
)
plotPCA(
    umi.qc[endog_genes, ],
    colour_by = "batch",
    size_by = "total_features",
    shape_by = "individual",
    exprs_values = "norm_exprs"
)

Downsampling

最后要介绍的这个去除文库大小差异的方法是从大的文库样本里面随机抽取部分reads使之文库大小缩减到跟其它文库一致。它的优点是抽样过程中会造成一些基因表达量为0,这样人为创造了dropout情况,弥补了系统误差。但是有个很重要的缺点,就是每次抽样都是随机的,这样结果无法重复,一般需要多次抽样保证结果的鲁棒性。

抽样函数如下:

Down_Sample_Matrix <-
function (expr_mat) 
{
    min_lib_size <- min(colSums(expr_mat))
    down_sample <- function(x) {
        prob <- min_lib_size/sum(x)
        return(unlist(lapply(x, function(y) {
            rbinom(1, y, prob)
        })))
    }
    down_sampled_mat <- apply(expr_mat, 2, down_sample)
    return(down_sampled_mat)
}

抽样后的counts矩阵赋值给SCESet对象的新的属性。

norm_counts(umi.qc) <- log2(Down_Sample_Matrix(counts(umi.qc)) + 1)
plotPCA(
    umi.qc[endog_genes, ],
    colour_by = "batch",
    size_by = "total_features",
    shape_by = "individual",
    exprs_values = "norm_counts"
)
umi.qc
## SCESet (storageMode: lockedEnvironment)
## assayData: 14063 features, 654 samples 
##   element names: counts, exprs, log2_counts, norm_counts, norm_cpm, norm_exprs 
## protocolData: none
## phenoData
##   rowNames: NA19098.r1.A01 NA19098.r1.A02 ... NA19239.r3.H12 (654
##     total)
##   varLabels: individual replicate ... size_factor (48 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: ENSG00000237683 ENSG00000187634 ... ERCC-00171
##     (14063 total)
##   fvarLabels: mean_exprs exprs_rank ... use (13 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:

同样的,也可视化一下表达矩阵,看看这个normalization的效果如何。

plotRLE(
    umi.qc[endog_genes, ], 
    exprs_mats = list(Raw = "log2_counts", DownSample = "norm_counts"),
    exprs_logged = c(TRUE, TRUE),
    colour_by = "batch"
)

文中提到的测试数据:

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-20

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏专知

【专知-PyTorch手把手深度学习教程08】NLP-PyTorch: 用字符级RNN生成名字

【导读】主题链路知识是我们专知的核心功能之一,为用户提供AI领域系统性的知识学习服务,一站式学习人工智能的知识,包含人工智能( 机器学习、自然语言处理、计算机视...

4535
来自专栏专知

【论文推荐】最新六篇情感分析相关论文—深度上下文、支持向量机、两级LSTM、多模态情感分析、软件工程、代码混合

【导读】专知内容组整理了最近六篇情感分析(Sentiment Analysis)相关文章,为大家进行介绍,欢迎查看! 1. Deep contextualize...

1.4K13

预测模型的计算时间

在周二我给精算师上的5小时机器学习速成课结束时,皮埃尔问了我一个有趣问题,是关于不同技术的计算时间的。我一直在介绍各种算法的思想,却忘了提及计算时间。我想在数据...

4547
来自专栏深度学习自然语言处理

【概率笔记】条件概率这样学才快啦

比如,一个上学期间整天鬼混的学沫,根本就不好好学习,对于他而言,选择题的四个选项ABCD被他选取的概率就为1/4。而对于大学霸来说,题题都会,那么他选取每一个选...

1033
来自专栏Echo is learning

machine learning 之 Anomaly detection

1321
来自专栏数据之美

相似文档查找算法之 simHash 简介及其 java 实现

传统的 hash 算法只负责将原始内容尽量均匀随机地映射为一个签名值,原理上相当于伪随机数产生算法。产生的两个签名,如果相等,说明原始内容在一定概 率 下是相...

85210
来自专栏小小挖掘机

DQN三大改进(二)-Prioritised replay

Prioritised replay原文:https://arxiv.org/pdf/1511.05952.pdf 代码地址:https://github.co...

1.5K4
来自专栏IT大咖说

死磕一周算法,我让服务性能提高 50%

内容来源:haolujun,https://www.cnblogs.com/haolujun/p/9527776.html

1915
来自专栏IT派

Python 实现简单的导弹自动追踪

自动追踪算法,在我们设计2D射击类游戏时经常会用到,这个听起来很高大上的东西,其实也并不是军事学的专利,在数学上解决的话需要去解微分方程,

2023
来自专栏ATYUN订阅号

使用Scikit-Learn进行命名实体识别和分类(NERC)

命名实体识别和分类(NERC)是识别名称等信息单元的过程(包括人员,组织和位置名称),以及包括非结构化文本中的时间,日期,钱和百分比表达式等数值表达式。目标是开...

1.6K6

扫码关注云+社区

领取腾讯云代金券