NC单细胞文章复现(二)：标准化
分享是一种态度
我们开始看看作者是怎么标准化数据的。考虑到各种混杂因素对单个细胞中定量表达水平的影响，混杂因素包括dropout events的频率、单个细胞中的低表达量mRNA、不同类型细胞捕获的高可变性或技术噪声，因此我们需要对QC后的数据进行标准化。接下来我们的分析结果会与文章有所不同，但是我们重点要学的是对数据的处理过程，所以我们继续学习下去。
二.标准化
标准化分为三个步骤进行，
- (1) 使用Census algorithm将TPM值转换为relative counts(用Monocle包的函数Relative2abs)
- (2) 对relative counts进行去卷积化（scran包）；
- (3) 用RUVSeq去除scrannormalized Census counts 额外的变异噪音。
#我们对过滤好的1326个细胞，总共13280个基因重新进行QC
#同样地，把cacultecell改成addPerCellQC
sceset <- sceset_all[keep.genes, ]
sceset<- addPerCellQC(sceset,exprs_values = "exprs",
detection_limit = min_thresh_log_tpm,flatten=FALSE,
subsets= list("regular" = which(colData(sceset)$pool_H12 == 0), "pool" = which(colData(sceset)$pool_H12 == 1)))
#创建对象
min_thresh_tpm <- 1
pd_2 <- new("AnnotatedDataFrame", data = as.data.frame(colData(sceset)))
featureNames(pd_2) <- rownames(pd_2)
fd_2 <- new("AnnotatedDataFrame", data = as.data.frame(rowData(sceset)))
featureNames(fd_2) <- rowData(sceset)$gene_short_name
HSMM <- newCellDataSet(assays(sceset)$tpm, phenoData = pd_2, featureData = fd_2,
expressionFamily = tobit(), lowerDetectionLimit = min_thresh_tpm)
HSMM <- detectGenes(HSMM, min_expr = min_thresh_tpm)
colData(sceset)$num_genes_expressed <- pData(HSMM)$num_genes_expressed
接下来这一步，使用relative2abs函数中的Census algorithm把TPM values 转换为 relative counts。也就是每个细胞中的TPM值除以每个细胞中估计的Total mRNA，重新计算TPM值。这样子，TPM值被转换成负二项分布的relative Census counts，这样相比于edgeR，DESeq2和SCDE 能更好地反映了单个细胞中mRNA的实际数量。我们看看monocle包中对relative2abs的描述。
对于我们来说，学到了另一种转录组测序数据处理方法。比如是否可以拿GEO/TCGA或者自己的数据TPM/FPKM转换成relative counts，然后进行下游分析？这个留给大家自己探索了....
# 将 TPM values转换为relative counts
rpc_matrix <- relative2abs(HSMM)
#创建对象
HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"),
phenoData = phenoData(HSMM),
featureData = featureData(HSMM),
lowerDetectionLimit = min_thresh_log_tpm,
expressionFamily = negbinomial.size())
#归一化处理
HSMM <- estimateSizeFactors(HSMM)
HSMM <- estimateDispersions(HSMM)
assays(sceset)$monocle <- as.matrix(exprs(HSMM))
assays(sceset)$log2_tpm <- log2(assays(sceset)$tpm + 1)
assays(sceset)$log2_monocle <- as.matrix(log2(exprs(HSMM) + 1))
assays(sceset)$log2_tpm_median <- t(scale(t(assays(sceset)$exprs), scale = FALSE))
接下来这一步会过将1326个samples过滤掉137个，剩下1189个samples。但是我们按照作者的代码跑的话，无法过滤掉这137个样本。我的理解是可能是在前面的calculateQCMetrics换成了addPerCellQC，QC的算法略有改变，对结果会有一点影响。
# 将每个细胞的表达量除以其因子大小，减少了干扰效应
exprs_filtered <- t(t(exprs(HSMM)/pData(HSMM)$Size_Factor))
#挑选基因值不是0的，此时没有1个samples的值为0。
nz_genes <- which(exprs_filtered != 0)
#再次log2转换
exprs_filtered[nz_genes] <- log2(exprs_filtered[nz_genes] + 1)
assays(sceset)$norm_log2_monocle <- as.matrix(exprs_filtered)
#此时，我们得到的samples仍然是1376个
length(colnames(as.matrix(exprs_filtered)))
> length(colnames(as.matrix(exprs_filtered)))
[1] 1326
接下来使用scan包来remove additional sources of unwanted variation。
RUVgb函数认为在样本中管家基因具有近似恒定表达量，并使用广义线性模型去回归管家基因表达的估计变异。
sceset_mon_norm <- sceset
assays(sceset_mon_norm)$counts <- assays(sceset)$monocle
clusters_mon_norm <- quickCluster(sceset_mon_norm, min.size = 100)
sceset_mon_norm <- computeSumFactors(sceset_mon_norm, cluster = clusters_mon_norm, sizes = 200, positive = TRUE)
if (length(which(sizeFactors(sceset_mon_norm) == 0)) > 0) {
sceset_mon_norm <- sceset_mon_norm[,-which(sizeFactors(sceset_mon_norm) == 0)]
assays(sceset_mon_norm)$log2_tpm_median <- t(scale(t(assays(sceset_mon_norm)$exprs), scale = FALSE))
}
sceset_mon_norm <- logNormCounts(sceset_mon_norm)
assays(sceset_mon_norm)$exprs <- assays(sceset_mon_norm)$logcounts
很遗憾，我们没能揪出那 137个samples，将1326个samples过滤成1189个samples，但是作者的代码大概意思懂了就好，毕竟单细胞发展这么快，相应的R包更新也快，有些函数里面的算法可能已经改变了，这样的话有时候就会没法完全重现作者的数据。作者也在给出了这样的说明：
# expressionFamily = gaussianff()) not supported from VGAM 1.0.6 onwards
# original code (working slightly differently now with the more recent versions of quickCluster and computeSumFactors from scran)
# due to updates/changes in quickCluster and computeSumFactors (at a minimum), the resulting values are slightly different
# than the ones obtained when carrying out the original analysis. to reproduce the results from the paper, we can simply
# read in the normalized data and the phenoData, as preprocessed in the original analysis
可喜的是，我们学习到了作者是如何对数据进行质控的，这一部分的内容在文献方法那里有详细的描述，有兴趣的小伙伴们可以去钻研。并且我们可以拿到作者上传到GSE118390上的已经normalized data来继续下游分析学习，下一节将学习细胞亚群的注释和一张复杂的热图绘制，我会结合文献的内容跟大家分享。
往期回顾
OSCA单细胞数据分析笔记8—Dimensionality reduction
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
