normalize
比较不同normalization方法
normalization目的是消除混杂因素的影响，这里主要是library size的影响。 bulk RNA-seq我们常用CPM（RPKM,FPKM,TPM）、SF、UQ、TMM来进行normmalization，但这些方法并不一定适用于single-cell RNA-seq数据，因为它们基于一个假设，即这些细胞中的RNA的量是一样多的。
library(scater)
## Warning: package 'scater' was built under R version 3.5.2
library(scran)
umi.qc <- readRDS("D:/paopaoR/normalize/tung/umi.qc.rds")
umi.qc
## class: SingleCellExperiment
## dim: 13904 639
## metadata(0):
## assays(3): counts logcounts_raw logcounts
## rownames(13904): ENSG00000237683 ENSG00000187634 ... ERCC-00170
## ERCC-00171
## rowData names(10): is_feature_control is_feature_control_MT ...
## log10_total_counts use
## colnames(639): NA19098.r1.A01 NA19098.r1.A02 ... NA19239.r3.H11
## NA19239.r3.H12
## colData names(51): individual replicate ...
## pct_counts_in_top_500_features_ERCC use
## reducedDimNames(0):
## spikeNames(2): ERCC MT
RAW
assay(umi.qc,"logcounts") <- log2(counts(umi.qc)+1)
assays(umi.qc)$logcounts[1:5,1:5]
## NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683 0 0.000000 0
## ENSG00000187634 0 0.000000 0
## ENSG00000188976 2 2.807355 1
## ENSG00000187961 0 0.000000 0
## ENSG00000187608 1 3.000000 0
## NA19098.r1.A04 NA19098.r1.A05
## ENSG00000237683 1 0.000000
## ENSG00000187634 0 0.000000
## ENSG00000188976 2 2.321928
## ENSG00000187961 0 0.000000
## ENSG00000187608 2 1.000000
plotPCA(
umi.qc,
colour_by = "batch",
size_by = "total_features_by_counts",
shape_by = "individual"
)
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.
image.png
plotRLE(umi.qc,exprs_values = "logcounts",exprs_logged = TRUE,colour_by = "batch")
image.png
CPM
最简单直接的方法，计算每百万数（CPM），方法是将每一列除以其总数，然后*10e6。 潜在缺陷：如果样本中包含的基因在细胞中高度表达或者差异表达会严重影响那些低表达的基因。 RPKM，FPKM，TPM是CPM的变体，对基因/转录本的长度进一步调节计数。
assay(umi.qc,"logcounts") <- log2(calculateCPM(umi.qc,use_size_factors = FALSE)+1)
assays(umi.qc)$logcounts[1:5,1:5]
## NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683 0.000000 0.000000 0.000000
## ENSG00000187634 0.000000 0.000000 0.000000
## ENSG00000188976 5.597831 6.567817 4.581699
## ENSG00000187961 0.000000 0.000000 0.000000
## ENSG00000187608 4.071250 6.788035 0.000000
## NA19098.r1.A04 NA19098.r1.A05
## ENSG00000237683 4.29030 0.000000
## ENSG00000187634 0.00000 0.000000
## ENSG00000188976 5.82525 5.845139
## ENSG00000187961 0.00000 0.000000
## ENSG00000187608 5.82525 3.918529
plotPCA(
umi.qc,
colour_by = "batch",
size_by = "total_features_by_counts",
shape_by = "individual"
)
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.
image.png
plotRLE(umi.qc,exprs_values = "logcounts",exprs_logged = TRUE,colour_by = "batch")
image.png
RLE (SF)
size factor(SF):首先计算所有细胞中每个基因的几何平均值，每个细胞的size factor是基因表达与基因几何平均值比值的中位数。 缺点：使用的是几何平均值，所以只能在所有细胞中使用非零表达值的基因进行计算，因此不适合进行大规模低深度的scRNAseq。在scater里面包装成normalizeEXprs()函数，但是我找了半天都报错说不存在这个函数，所以直接用function吧……
calc_sf <- function(exprmat){
geommeans <- exp(rowMeans(log(exprmat)))
SF <- function(cnts){
median((cnts/geommeans)[(is.finite(geommeans)&geommeans >0)])#非零
}
norm_factor <- apply(exprmat, 2,SF)
return(t(t(exprmat)/norm_factor))
}
assay(umi.qc,"logcounts") <- calc_sf(counts(umi.qc))
assays(umi.qc)$logcounts[1:5,1:5]
## NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683 0.0000000 0.000000 0.000000
## ENSG00000187634 0.0000000 0.000000 0.000000
## ENSG00000188976 2.9390913 5.641494 1.425718
## ENSG00000187961 0.0000000 0.000000 0.000000
## ENSG00000187608 0.9796971 6.581743 0.000000
## NA19098.r1.A04 NA19098.r1.A05
## ENSG00000237683 1.132252 0.0000000
## ENSG00000187634 0.000000 0.0000000
## ENSG00000188976 3.396756 3.4009587
## ENSG00000187961 0.000000 0.0000000
## ENSG00000187608 3.396756 0.8502397
plotPCA(
umi.qc,
colour_by = "batch",
size_by = "total_features_by_counts",
shape_by = "individual"
)
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.
image.png
plotRLE(umi.qc,exprs_values = "logcounts",exprs_logged = TRUE,colour_by = "batch")
image.png
UQ
upperquartile（UQ上四分位数），每个列除以每个library计数的75%分位数。通常，计算出来的分位数由细胞间的中位数进行缩放，以保持绝对表达水平相对一致。 缺点：在低深度scRNAseq实验中，大量的未检测到的基因可能导致75%分位数为0（或者接近于0）。
calc_sf <- function(exprmat){
UQ <- function(x){
quantile(x[x>0],0.75)
}
uq <- apply(exprmat,2,UQ)
norm_factor <- uq/median(uq)
return(t(t(exprmat)/norm_factor))
}
assay(umi.qc,"logcounts") <- calc_sf(counts(umi.qc))
assays(umi.qc)$logcounts[1:5,1:5]
## NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683 0.000000 0.000000 0.0
## ENSG00000187634 0.000000 0.000000 0.0
## ENSG00000188976 3.500000 7.000000 1.4
## ENSG00000187961 0.000000 0.000000 0.0
## ENSG00000187608 1.166667 8.166667 0.0
## NA19098.r1.A04 NA19098.r1.A05
## ENSG00000237683 1.4 0
## ENSG00000187634 0.0 0
## ENSG00000188976 4.2 4
## ENSG00000187961 0.0 0
## ENSG00000187608 4.2 1
plotPCA(
umi.qc,
colour_by = "batch",
size_by = "total_features_by_counts",
shape_by = "individual"
)
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.
image.png
plotRLE(umi.qc,exprs_values = "logcounts",exprs_logged = TRUE,colour_by = "batch")
image.png
TMM
m值是细胞间log2fc的变化。使用一个细胞作为参考，然后将每个细胞的m值与该参考进行比较。然后，通过删除顶部和底部~30%来调整这些值，然后通过加权计算剩余值的平均值，以考虑对数尺度对方差的影响。每个非参考单元格乘以计算因子。 潜在问题：修剪后留下的非零基因不足，以及大多数基因没有差异表达的假设。
downsampling
对表达矩阵进行下采样，使每个细胞的molecular数大致相同。 好处：下采样将引入零值，从而消除由于检测到的基因数量不同而产生的任何偏差。 缺点：这个过程不是确定性的，因此每次运行下采样时得到的表达矩阵是略有不同的。因此通常必须在多个下采样上运行分析，以确保结果是健壮的。
scater
sizeFactors(umi.qc) <- librarySizeFactors(umi.qc)
umi.qc <- computeSpikeFactors(umi.qc)
umi.qc <- normalize(umi.qc)
assays(umi.qc)$logcounts[1:5,1:5]
## NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683 0.0000000 0.000000 0.000000
## ENSG00000187634 0.0000000 0.000000 0.000000
## ENSG00000188976 1.7106114 2.702905 1.065341
## ENSG00000187961 0.0000000 0.000000 0.000000
## ENSG00000187608 0.8136601 2.893291 0.000000
## NA19098.r1.A04 NA19098.r1.A05
## ENSG00000237683 0.9627605 0.0000000
## ENSG00000187634 0.0000000 0.0000000
## ENSG00000188976 1.9437739 2.0912404
## ENSG00000187961 0.0000000 0.0000000
## ENSG00000187608 1.9437739 0.8601966
plotPCA(
umi.qc,
colour_by = "batch",
size_by = "total_features_by_counts",
shape_by = "individual"
)
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.
image.png
plotRLE(umi.qc,exprs_values = "logcounts",exprs_logged = TRUE,colour_by = "batch")
image.png
scran
scran包实现了专门用于单细胞数据的CPM变体。简单说来，这种方法通过将cells集中在一起计算一个标准化因子（类似于CPM）来处理每个cell中的大量零值问题。
qclust <- quickCluster(umi.qc, min.size = 30)
umi.qc <- computeSpikeFactors(umi.qc)
umi.qc <- computeSumFactors(umi.qc, sizes = 15, clusters = qclust)
umi.qc <- normalize(umi.qc)
assays(umi.qc)$logcounts[1:5,1:5]
## 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.078042 2.923782 1.406737
## ENSG00000187961 0.000000 0.000000 0.000000
## ENSG00000187608 1.052494 3.118755 0.000000
## NA19098.r1.A04 NA19098.r1.A05
## ENSG00000237683 1.203788 0.000000
## ENSG00000187634 0.000000 0.000000
## ENSG00000188976 2.295814 2.329960
## ENSG00000187961 0.000000 0.000000
## ENSG00000187608 2.295814 1.005025
plotPCA(
umi.qc,
colour_by = "batch",
size_by = "total_features_by_counts",
shape_by = "individual"
)
## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.
image.png
plotRLE(umi.qc,exprs_values = "logcounts",exprs_logged = TRUE,colour_by = "batch")
image.png
