比較不同的對單細胞轉錄組資料normalization方法

阿新 • • 發佈：2022-05-03

使用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

比較不同的對單細胞轉錄組資料normalization方法

使用CPM去除文庫大小影響

RPKM, FPKM and TPM去除基因或者轉錄本長度影響

適用於bulk RNA-seq的normalization方法

特意為single-cell RNA-seq data 開發的normalization方法

工作環境

載入測試資料

實踐

Raw

CPM

TMM

觀察變化規律

scran

Size-factor (RLE)

Upperquantile

Downsampling

抽樣後的counts矩陣賦值給SCESet物件的新的屬性。

比較不同的對單細胞轉錄組資料normalization方法

比較不同的對單細胞轉錄組資料尋找差異基因的方法

比較不同的對單細胞轉錄組資料聚類的方法

比較不同單細胞轉錄組資料尋找features方法

SC2disease：人類疾病的單細胞轉錄組的人工收集資料庫

轉錄組資料拼接之應用篇

單細胞轉錄組3大R包之monocle2

單細胞轉錄組3大R包之Seurat

單細胞轉錄組3大R包之scater

高階轉錄組分析和R語言資料視覺化第十二期（線上線下同時開課）

一個植物轉錄組專案的實戰

js實現樹形資料轉成扁平資料的方法示例

mysql實用技巧之比較兩個表是否有不同資料的方法分析

MySql查詢兩張相同表，合併成一組資料，並區分資料的不同

doccano標註後的序列標註任務資料轉錄為BIO形式

python 如何判斷一組資料是否符合正態分佈

efcore 跨表查詢，實現一個介面內查詢兩個不同資料庫裡各自的表資料

Django def clean()函式對錶單中的資料進行驗證操作

基於python實現計算兩組資料P值

Java對世界不同時區timezone之間時間轉換的處理方法

比較不同的對單細胞轉錄組資料normalization方法

使用CPM去除文庫大小影響

RPKM, FPKM and TPM去除基因或者轉錄本長度影響

適用於bulk RNA-seq的normalization方法

特意為single-cell RNA-seq data 開發的normalization方法

工作環境

載入測試資料

實踐

Raw

CPM

TMM

觀察變化規律

scran

Size-factor (RLE)

Upperquantile

Downsampling

抽樣後的counts矩陣賦值給SCESet物件的新的屬性。

相關推薦