什么？你做的差异基因方法不合适？
單細胞系列教程
收藏 北大生信平臺” 單細胞分析、染色質分析” 視頻和PPT分享
Science:?小鼠腎臟單細胞轉錄組+突變分析揭示腎病潛在的細胞靶標
10X單細胞測序分析軟件:Cell?ranger，從拆庫到定量
Hemberg-lab單細胞轉錄組數據分析（一）- 引言
Hemberg-lab單細胞轉錄組數據分析（二）- 實驗平臺
Hemberg-lab單細胞轉錄組數據分析（三）- 原始數據質控
Hemberg-lab單細胞轉錄組數據分析（四）- 文庫拆分和細胞鑒定
Hemberg-lab單細胞轉錄組數據分析（五）- STAR, Kallisto定量
Hemberg-lab單細胞轉錄組數據分析（六）- 構建表達矩陣，UMI介紹
Hemberg-lab單細胞轉錄組數據分析（七）- 導入10X和SmartSeq2數據Tabula Muris
Hemberg-lab單細胞轉錄組數據分析（八）- Scater包輸入導入和存儲
Hemberg-lab單細胞轉錄組數據分析（九）- Scater包單細胞過濾
Hemberg-lab單細胞轉錄組數據分析（十）- Scater基因評估和過濾
Hemberg-lab單細胞轉錄組數據分析（十一）- Scater單細胞表達譜PCA可視化
Hemberg-lab單細胞轉錄組數據分析（十二）- Scater單細胞表達譜tSNE可視化
單細胞分群后，怎么找到Marker基因定義每一類群？
DESeq2差異基因分析和批次效應移除
如何火眼金睛鑒定那些單細胞轉錄組中的混雜因素
表達數據標準化理論
簡介
上一步，我們鑒定出了重要的干擾因素和解釋變量可能對表達定量帶來影響。scater允許在后續統計模型中引入這些變量來屏蔽技術操作帶來的影響，或者可以給函數normaliseExprs()提供一個設計矩陣design matrix來直接移除干擾因素的影響。在這一章先不涉及這些。
相反，我們探索下簡單的量化因子size-factor標準化如何在校正文庫大小的同時移除部分干擾因素引入的檢測偏差。
文庫大小
scRNA-seq數據的文庫大小變化很大是因為混樣測序，來源于每個細胞的總reads叔差別較大。一些定量方法，如Cufflinks, RSEM) 在估算基因表達時已經考慮了文庫大小的影響因此不需要這一步標準化。
然而，如果采用的是其它的定量方法就必須首先通過某種方法估算一起比較的每個樣品的文庫大小也稱為量化因子 (ormalization factor)，然后原始表達量乘以或除以量化因子矩陣獲得標準化后的表達結果。很多開發出來用于bulk RNA-seq的結果，也可以用于scRNA-seq，比如UQ, SF, CPM, RPKM, FPKM, TPM.
標準化方法
CPM
最簡單的標準化方式是原始reads除以樣品總的可用reads數乘以1,000,000獲得每百萬reads的count數 (counts per million(__CPM__)。因為考慮的是總的細胞內的reads，計算總的reads數時只考慮內源性基因而排除spike-ins部分的reads。CPM計算的R代碼是:
calc_cpm <- function (expr_mat, spikes = NULL){norm_factor <- colSums(expr_mat[-spikes,])return(t(t(expr_mat)/norm_factor)) * 10^6 }這種計算方式的缺點是容易受到極高表達且在不同樣品中存在差異表達的基因的影響；這些基因的打開或關閉會影響到細胞中總的分子數目，可能導致這些基因標準化之后就不存在表達差異了，而原本沒有差異的基因標準化 之后卻有差異了。
RPKM、FPKM和TPM是CPM按照基因或轉錄本長度歸一化后的表達，也會受到這一影響。
為了解決這一問題，研究者提出了其它的標準化方法。
Relative Log Expression RLE (SF)
量化因子 (size factor, SF)是由DESeq [@Anders2010-jr]提出的。其方法是首先計算每個基因在所有樣品中表達的幾何平均值。每個細胞的量化因子(size factor)是所有基因與其在所有樣品中的表達值的幾何平均值的比值的中位數。由于幾何平均值的使用，只有在所有樣品中表達都不為0的基因才能用來計算，所以不適合大批量低深度的scRNA-seq數據。edgeR & scater 把這種方法稱為”relative log expression” (RLE)。其在R中計算函數是:
calc_sf <- function (expr_mat, spikes=NULL){geomeans <- exp(rowMeans(log(expr_mat[-spikes,])))SF <- function(cnts){median((cnts/geomeans)[(is.finite(geomeans) & geomeans >0)])}norm_factor <- apply(expr_mat[-spikes,],2,SF)#return(t(t(expr_mat)/norm_factor))return(norm_factor) }上四分位數 (upperquartile, UQ)
上四分位數 (upperquartile, UQ)是樣品中所有基因的表達除以該樣品處于上四分位數的基因的表達值 [@Bullard2010-eb]。同時為了保證絕對表達水平的相對穩定，計算得到的上四分位數值要除以所有樣品中上四分位數值的中位數。對低深度scRNA-seq數據，這個方法的一個缺點是可能處于上四分位數的基因的表達值為0或接近0。這個限制可以通過采用更高的分位數如99%分位數 (scater的默認值)或排除表達值為0的基因后剩余基因的上四分位數。示例如下：
calc_uq <- function (expr_mat, spikes=NULL){UQ <- function(x) {quantile(x[x>0],0.75)}uq <- unlist(apply(expr_mat[-spikes,],2,UQ))norm_factor <- uq/median(uq)#return(t(t(expr_mat)/norm_factor))return(norm_factor) }TMM (M-值的加權截尾均值)
TMM是M-值的加權截尾均值 [@Robinson2010-hz]。選定一個樣品為參照，其它樣品中基因的表達相對于參照樣品中對應基因表達倍數的log2值定義為M-值。隨后去除M-值中最高和最低的30%，剩下的M值計算加權平均值。每一個非參照樣品的基 因表達值都乘以計算出的TMM。這個方法的兩個可能問題是，一是Trim后沒有足夠的非0基因，另外該方法假設大部分基因的表達是沒有差異的。
scran
scran采用為scRNA-seq設計的CPM方法的變種 [@LLun2016-pq]. 該方法通過把多組細胞合并到一起來屏蔽較多的0值問題，然后采用類似_CPM的方式計算標準化因子。因為一個細胞會出現在多個合并的集合里面 (pool)，細胞特異的因子可以采用線性代數從非特異性因子中去卷積計算得來。(Since each cell is found in many different pools, cell-specific factors can be deconvoluted from the collection of pool-specific factors using linear algebra.)
Downsampling
最后一個校正文庫大小的方式是對表達矩陣進行向下抽樣使得每個細胞檢測到的總分子數相同。這個方法的優勢是計算過程中會引入0值進而消除不同細胞檢測到的基因數不同引入的偏差。該方法最大的缺點是其非確定性，每次downsampling獲得的表達矩陣都會有些細微不同。通常需要重復多次保證結果的穩定性。其操作代碼是:
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) }標準化效果評估
我們將通過PCA方法和計算樣品范圍的 relative log expression (scater::plotRLE())來比較不同標準化方法的效果。含有更多reads的細胞，其大部分基因的表達比所有細胞的中值表達水平也更高，得到RLE值為正值；含有更少reads的細胞，其大部分基因的表達比所有細胞的中值表達水平更低，得到的RLE為負值。而標準化后的RLE值應該為0。計算_RLE_的函數如下：
calc_cell_RLE <- function (expr_mat, spikes = NULL) {RLE_gene <- function(x) {if (median(unlist(x)) > 0) {log((x + 1)/(median(unlist(x)) + 1))/log(2)}else {rep(NA, times = length(x))}}if (!is.null(spikes)) {RLE_matrix <- t(apply(expr_mat[-spikes, ], 1, RLE_gene))}else {RLE_matrix <- t(apply(expr_mat, 1, RLE_gene))}cell_RLE <- apply(RLE_matrix, 2, median, na.rm = T)return(cell_RLE) }注意 1： RLE, TMM, 和 UQ 量化因子方法是為bulk RNA-seq開發的，依賴于實驗條件，不一定適合scRNA-seq數據，尤其是潛在假設不成立時。
注意 2： scater提供了一個函數calcNormFactors實現了幾個文庫大小標準化方法方便后續使用。
注意 3： edgeR對一些標準化方法做了額外的調整使得它可能獲得與原始方法不同的結果。edgeR和scater的”RLE”方法基于DESeq的量化因子方法，但可能計算結果會不同于DESeq/DESeq2包的estimateSizeFactorsForMatrix計算的結果。另外，一些版本的edgeR只有在所有細胞的lib.size為設置為1時才計算標準化因子。
注意 4： CPM標準化使用的是scater包的calculateCPM()函數。scater包的normaliseExprs()函數用于 RLE, UQ 和 TMM 標準化計算。scran 標準化使用的是scran包計算量化因子 (基于SingleCellExperiment數據對象)和scater包的normalize()函數。所有標準化函數把計算結果存儲到SCE對象的logcounts通道 (slot)。downsampling 標準化使用的是前面展示的方法。
標準化實戰 (UMI)
繼續使用tung數據集進行后續分析：
library(scRNA.seq.funcs) library(scater) library(scran) options(stringsAsFactors = FALSE) set.seed(1234567) # umi <- readRDS("tung/umi.rds") # umi.qc <- umi[rowData(umi)$use, colData(umi)$use] # endog_genes <- !rowData(umi.qc)$is_feature_controlumi <- readRDS("tung/umi.rds") umi_endog_genes <- !rowData(umi)$is_feature_control umi_endog <- umi[umi_endog_genes,] umi.qc <- umi[rowData(umi)$use, colData(umi)$use] umi_qc_endog_genes <- !rowData(umi.qc)$is_feature_control umi.qc_endog <- umi.qc[umi_qc_endog_genes,]原始值
umi.qc_endog <- runPCA(umi.qc_endog, exprs_values = "logcounts_raw") scater::plotPCA(umi.qc_endog,by_exprs_values = "logcounts_raw",colour_by = "batch",size_by = "total_features_by_counts",shape_by = "individual" ) # plotPCA( # umi.qc[endog_genes, ], # exprs_values = "logcounts_raw", # colour_by = "batch", # size_by = "total_features", # shape_by = "individual" # )CPM
logcounts(umi.qc) <- log2(calculateCPM(umi.qc, use_size_factors = FALSE) + 1)scater::plotPCA(umi.qc[umi_qc_endog_genes,],by_exprs_values = "logcounts",colour_by = "batch",size_by = "total_features_by_counts",shape_by = "individual" )library(cowplot) # plotRLE( # umi.qc[endog_genes,], # exprs_mats = list(Raw = "logcounts_raw", CPM = "logcounts"), # exprs_logged = c(TRUE, TRUE), # colour_by = "batch" # )cpmlogcounts <- plotRLE(umi.qc[umi_qc_endog_genes,],exprs_values = "logcounts",exprs_logged = c(TRUE),colour_by = "batch" ) + ggtitle("Log CPM")rawlogcounts <- plotRLE(umi.qc[umi_qc_endog_genes,],exprs_values = "logcounts_raw",exprs_logged = c(TRUE),colour_by = "batch" ) + ggtitle("Raw log count")plot_grid(rawlogcounts, cpmlogcounts, nrow=2, ncol=1)Size-factor (RLE)
# umi.qc <- normaliseExprs( # umi.qc, # method = "RLE", # feature_set = umi_qc_endog_genes, # return_log = TRUE, # return_norm_as_exprs = TRUE # ) sizeFactors(umi.qc) <- calc_sf(counts(umi.qc),spikes=rowData(umi.qc)$is_feature_control_ERCC) umi.qc <- normalize(umi.qc)scater::plotPCA(umi.qc[umi_qc_endog_genes,],by_exprs_values = "logcounts",colour_by = "batch",size_by = "total_features_by_counts",shape_by = "individual" )rlelogcounts <- plotRLE(umi.qc[umi_qc_endog_genes,],exprs_values = "logcounts",exprs_logged = c(TRUE),colour_by = "batch" ) + ggtitle("RLE size factor")plot_grid(rawlogcounts, rlelogcounts, nrow=2, ncol=1)Upperquantile
# umi.qc <- normaliseExprs( # umi.qc, # method = "upperquartile", # feature_set = endog_genes, # p = 0.99, # return_log = TRUE, # return_norm_as_exprs = TRUE # ) sizeFactors(umi.qc) <- calc_uq(counts(umi.qc),spikes=rowData(umi.qc)$is_feature_control_ERCC) umi.qc <- normalize(umi.qc) plotPCA(umi.qc[umi_qc_endog_genes, ],colour_by = "batch",size_by = "total_features_by_counts",shape_by = "individual" )uqlogcounts <- plotRLE(umi.qc[umi_qc_endog_genes,],exprs_values = "logcounts",exprs_logged = c(TRUE),colour_by = "batch" ) + ggtitle("UQ")plot_grid(rawlogcounts, uqlogcounts, nrow=2, ncol=1)TMM
umi.qc <- normaliseExprs(umi.qc,method = "TMM",feature_set = umi_qc_endog_genes,return_log = TRUE,return_norm_as_exprs = TRUE ) plotPCA(umi.qc[umi_qc_endog_genes, ],colour_by = "batch",size_by = "total_features",shape_by = "individual" )plotRLE(umi.qc[umi_qc_endog_genes, ],exprs_mats = list(Raw = "logcounts_raw", TMM = "logcounts"),exprs_logged = c(TRUE, TRUE),colour_by = "batch" )scran
library(scran) qclust <- quickCluster(umi.qc, min.size = 30) umi.qc <- computeSumFactors(umi.qc, sizes = 15, clusters = qclust) umi.qc <- normalize(umi.qc) plotPCA(umi.qc[umi_qc_endog_genes, ],colour_by = "batch",size_by = "total_features_by_counts",shape_by = "individual" )scranlogcounts <- plotRLE(umi.qc[umi_qc_endog_genes,],exprs_values = "logcounts",exprs_logged = c(TRUE),colour_by = "batch" ) + ggtitle("Scran")plot_grid(rawlogcounts, scranlogcounts, nrow=2, ncol=1)scran有時會獲得負或零量化因子，這將會嚴重干擾標準化后的表達矩陣，需要采用下面的方法確認沒有問題：
summary(sizeFactors(umi.qc))## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.4671 0.7787 0.9483 1.0000 1.1525 3.1910對這個數據集來說，量化因子是合理。如果計算時發現scran給出的量化因子是非正值嘗試增加cluster和pool的大小，直到獲取正值。
Downsampling
logcounts(umi.qc) <- log2(Down_Sample_Matrix(counts(umi.qc)) + 1) plotPCA(umi.qc[umi_qc_endog_genes, ],colour_by = "batch",size_by = "total_features_by_counts",shape_by = "individual" )dslogcounts <- plotRLE(umi.qc[umi_qc_endog_genes,],exprs_values = "logcounts",exprs_logged = c(TRUE),colour_by = "batch" ) + ggtitle("Downsampling")plot_grid(rawlogcounts, dslogcounts, nrow=2, ncol=1)易生信系列培訓課程，掃碼獲取免費資料
更多閱讀
畫圖三字經?生信視頻?生信系列教程?
心得體會?TCGA數據庫?Linux?Python?
高通量分析?免費在線畫圖?測序歷史?超級增強子
生信學習視頻?PPT?EXCEL?文章寫作?ggplot2
海哥組學?可視化套路?基因組瀏覽器
色彩搭配?圖形排版?互作網絡
自學生信
后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集
總結
以上是生活随笔為你收集整理的什么？你做的差异基因方法不合适？的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 「R」ggplot2拼图包patchwo
- 下一篇： “嘿，我们又见面了！”