14.1 前言

本篇是細(xì)胞注釋章節(jié)的更詳細(xì)的延續(xù)，該章節(jié)已經(jīng)介紹了差異基因表達(dá)（DGE）作為用細(xì)胞類型注釋簇的工具续誉。在這里，我們重點(diǎn)關(guān)注更復(fù)雜的實(shí)驗(yàn)設(shè)計(jì)中差異基因表達(dá)測試的更高級用例徙缴，這些實(shí)驗(yàn)設(shè)計(jì)涉及一種或多種條件，例如疾病词顾、基因敲除或藥物。在這種情況下碱妆，我們通常對感興趣的條件和參考條件之間基因表達(dá)模式差異的大小和重要性感興趣。該參考可以是一切昔驱，但通常是健康樣本疹尾。這種統(tǒng)計(jì)測試可以應(yīng)用于任意組，但在單細(xì)胞RNA-Seq的情況下通常應(yīng)用于細(xì)胞類型水平骤肛。

差異基因表達(dá)分析試圖推斷任何比較組之間（通常在每種細(xì)胞類型的健康組和狀況組之間）統(tǒng)計(jì)上顯著過度表達(dá)或表達(dá)不足的基因纳本。

這種分析的結(jié)果可能是影響并可能解釋任何觀察到的表型的基因組。然后可以更仔細(xì)地檢查這些基因組腋颠，例如受影響的途徑或誘導(dǎo)的細(xì)胞間通訊變化繁成。

差異基因表達(dá)測試通常會返回每個(gè)比較條件下每個(gè)比較基因的log2倍數(shù)變化和調(diào)整后的p值。然后可以按p值對該列表進(jìn)行排序并進(jìn)行更詳細(xì)的研究淑玫。

流行的學(xué)生t檢驗(yàn)是進(jìn)行此類檢驗(yàn)的一種方法巾腕。然而，它沒有考慮到一些單細(xì)胞RNA-seq的特殊性絮蒿，例如來自dropout的過多零或需要復(fù)雜的實(shí)驗(yàn)設(shè)計(jì)尊搬。更具體地說，在不匯集跨基因信息的情況下土涝，很少有人有足夠的樣本數(shù)量來準(zhǔn)確估計(jì)方差佛寿。此外，原始計(jì)數(shù)從來都不是給定樣本中特定基因表達(dá)的絕對測量值但壮。每個(gè)基因的實(shí)際讀取數(shù)取決于文庫制備的效率冀泻、非編碼轉(zhuǎn)錄本的污染量和測序深度。因此蜡饵，它確實(shí)缺乏單細(xì)胞RNA-seq的敏感性和特異性弹渔，更不用說實(shí)驗(yàn)設(shè)計(jì)靈活性了。

因此验残，差異基因表達(dá)測試是一個(gè)經(jīng)典的生物信息學(xué)問題捞附，已經(jīng)被許多工具解決。一般來說您没，目前從兩個(gè)角度來解決這個(gè)問題鸟召，即樣本級視圖，其中表達(dá)被聚合以創(chuàng)建“偽批量”氨鹏，然后使用最初為批量表達(dá)樣本設(shè)計(jì)的方法進(jìn)行分析欧募，例如edgeR或DEseq2以及細(xì)胞水平視圖，其中使用廣義混合效應(yīng)模型（例如 MAST[Finak et al., 2015]或glmmTMB[Brooks et al., 2017]）對細(xì)胞進(jìn)行單獨(dú)建模仆抵。DGE工具的跨數(shù)據(jù)集的共識和穩(wěn)健性較低跟继。如前所述种冬，盡管單細(xì)胞數(shù)據(jù)包含技術(shù)噪聲偽影，例如丟失舔糖、零膨脹和高細(xì)胞間變異性娱两， 與專門為scRNA-seq數(shù)據(jù)設(shè)計(jì)的方法相比，為批量RNA-seq數(shù)據(jù)設(shè)計(jì)的方法表現(xiàn)良好金吗。發(fā)現(xiàn)單細(xì)胞特異性方法特別容易將高表達(dá)基因錯(cuò)誤地標(biāo)記為差異表達(dá)十兢。

在本篇中，我們演示如何使用兩種工具進(jìn)行差異表達(dá)分析：具有擬似然檢驗(yàn)的edgeR和具有隨機(jī)效應(yīng)的MAST摇庙。由于edgeR和MAST都是在R中實(shí)現(xiàn)的旱物，因此我們使用anndata2ri包能夠同時(shí)使用Python中的AnnData對象和R中的SingeCellExperiment對象。

對于這兩種方法中的每一種卫袒，我們展示了兩個(gè)用例：如何對完整數(shù)據(jù)進(jìn)行分析以及如何對多種細(xì)胞類型和一種特定細(xì)胞類型進(jìn)行測試宵呛。

14.2 環(huán)境配置

import warnings

warnings.filterwarnings("ignore")

import matplotlib.pyplot as plt
import seaborn as sns
import scanpy as sc
import pandas as pd
import numpy as np
import random
import sc_toolbox
import pertpy 

import rpy2.rinterface_lib.callbacks
import anndata2ri
import logging

from rpy2.robjects import pandas2ri
from rpy2.robjects import r

sc.settings.verbosity = 0
rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)

pandas2ri.activate()
anndata2ri.activate()

%load_ext rpy2.ipython

%%R
library(edgeR)
library(MAST)

14.3 數(shù)據(jù)集準(zhǔn)備

我們將使用Kang數(shù)據(jù)集，它是來自8名狼瘡患者INF-β治療前后6小時(shí)（總共16個(gè)樣本）的10倍基于scRNA-seq的外周血單核細(xì)胞 (PBMC) [Kang et al., 2018]夕凝。 干擾素β以天然成纖維細(xì)胞或重組制劑（干擾素β-1a和干擾素β-1b）的形式使用宝穗，并發(fā)揮與干擾素α類似的抗病毒和抗增殖特性。干擾素β已被批準(zhǔn)用于治療復(fù)發(fā)緩解型多發(fā)性硬化癥和繼發(fā)進(jìn)展型多發(fā)性硬化癥码秉。

首先讽营，我們加載完整的數(shù)據(jù)集。

adata = pertpy.data.kang_2018()
adata

AnnData object with n_obs × n_vars = 24673 × 15706
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters'
    var: 'name'
    obsm: 'X_pca', 'X_umap'

我們將需要.obs的label（其中包含條件標(biāo)簽）泡徙、replicate和cell_type列橱鹏。

adata.obs[:5]

    nCount_RNA  nFeature_RNA    tsne1   tsne2   label   cluster cell_type   replicate   nCount_SCT  nFeature_SCT    integrated_snn_res.0.4  seurat_clusters
index                                               
AAACATACATTTCC-1    3017.0  877 -27.640373  14.966629   ctrl    9   CD14+ Monocytes patient_1016    1704.0  711 1   1
AAACATACCAGAAA-1    2481.0  713 -27.493646  28.924885   ctrl    9   CD14+ Monocytes patient_1256    1614.0  662 1   1
AAACATACCATGCA-1    703.0   337 -10.468194  -5.984389   ctrl    3   CD4 T cells patient_1488    908.0   337 6   6
AAACATACCTCGCT-1    3420.0  850 -24.367997  20.429285   ctrl    9   CD14+ Monocytes patient_1256    1738.0  653 1   1
AAACATACCTGGTA-1    3158.0  1111    27.952170   24.159738   ctrl    4   Dendritic cells patient_1039    1857.0  928 12  12

我們需要使用原始計(jì)數(shù)，因此我們檢查.X確實(shí)包含原始計(jì)數(shù)并將它們放入AnnData對象的counts層中堪藐。

np.max(adata.X)

3828.0

adata.layers["counts"] = adata.X.copy()

我們有8名對照患者和8名疾病患者莉兰。

print(len(adata[adata.obs["label"] == "ctrl"].obs["replicate"].cat.categories))
print(len(adata[adata.obs["label"] == "stim"].obs["replicate"].cat.categories))

8
8

我們過濾含有少于200個(gè)基因的細(xì)胞以及在少于3個(gè)細(xì)胞中發(fā)現(xiàn)的基因，以進(jìn)行基本的質(zhì)量控制礁竞。

sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata

AnnData object with n_obs × n_vars = 24562 × 15701
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'n_genes'
    var: 'name', 'n_cells'
    obsm: 'X_pca', 'X_umap'
    layers: 'counts'

14.4 偽批量Pseudobulk

edgeR是在R中實(shí)現(xiàn)的差異基因表達(dá)測試工具糖荒，最初是為批量基因表達(dá)數(shù)據(jù)而設(shè)計(jì)的。它實(shí)現(xiàn)了基于負(fù)二項(xiàng)分布模捂、經(jīng)驗(yàn)貝葉斯估計(jì)捶朵、精確檢驗(yàn)、廣義線性模型 (GLM) 和擬似然檢驗(yàn)的廣泛統(tǒng)計(jì)方法狂男。

在這里综看，我們將使用擬似然檢驗(yàn)，因?yàn)樗紤]了離散度估計(jì)的不確定性岖食。相反红碑，精確檢驗(yàn)假設(shè)估計(jì)的離差是真實(shí)值，這可能會導(dǎo)致一些不準(zhǔn)確。此外析珊，擬似然GLM在實(shí)驗(yàn)設(shè)計(jì)方面更加靈活羡鸥。

由于edgeR是作為批量數(shù)據(jù)的DE分析方法引入的，因此我們首先需要從單細(xì)胞數(shù)據(jù)集中創(chuàng)建偽批量樣本忠寻。對于每位患者惧浴，我們通過聚集每個(gè)亞群的細(xì)胞并獲取該亞群內(nèi)的平均基因表達(dá)來為每種細(xì)胞類型創(chuàng)建1個(gè)偽批量樣本。

如果您正在使用的數(shù)據(jù)沒有重復(fù)奕剃，則為每個(gè)患者創(chuàng)建多個(gè)（例如2-3個(gè)）pseudobulks以考慮患者變異性可能會有所幫助赶舆。在這里，我們選擇為每個(gè)患者創(chuàng)建1個(gè)偽批量祭饭，因?yàn)閷τ诿總€(gè)患者，我們有2個(gè)重復(fù)叙量，一個(gè)對照倡蝙，一個(gè)刺激。

無論我們是只想對少數(shù)細(xì)胞亞群進(jìn)行分析并分別為每個(gè)細(xì)胞亞群擬合一個(gè)模型绞佩，還是為所有細(xì)胞亞群擬合一個(gè)模型寺鸥，我們首先需要準(zhǔn)備數(shù)據(jù)，定義一個(gè)函數(shù)來創(chuàng)建pseudobulks并運(yùn)行edgeR管道品山。首先胆建，我們準(zhǔn)備數(shù)據(jù)。

由于我們需要為每個(gè)患者狀況組合創(chuàng)建偽批量肘交，因此我們首先需要通過連接replicate和label來創(chuàng)建這樣的列笆载。

adata.obs["sample"] = [
    f"{rep}_{l}" for rep, l in zip(adata.obs["replicate"], adata.obs["label"])
]

我們需要清理細(xì)胞類型名稱，即用下劃線替換空格并刪除+符號涯呻，以避免Python到R的轉(zhuǎn)換問題凉驻。

adata.obs["cell_type"] = [ct.replace(" ", "_") for ct in adata.obs["cell_type"]]
adata.obs["cell_type"] = [ct.replace("+", "") for ct in adata.obs["cell_type"]]

我們需要將分類元數(shù)據(jù)設(shè)置為真正分類的才能創(chuàng)建偽批量。

adata.obs["replicate"] = adata.obs["replicate"].astype("category")
adata.obs["label"] = adata.obs["label"].astype("category")
adata.obs["sample"] = adata.obs["sample"].astype("category")
adata.obs["cell_type"] = adata.obs["cell_type"].astype("category")

現(xiàn)在复罐，讓我們定義將單個(gè)細(xì)胞聚合成偽重復(fù)所需的函數(shù)：

• aggregate_and_filter是一個(gè)函數(shù)涝登，它從原始單細(xì)胞AnnData對象中創(chuàng)建一個(gè)AnnData對象，該對象為指定子群提供一個(gè)偽復(fù)制效诅。在這里胀滚，我們還篩選出特定群體的細(xì)胞數(shù)量少于30個(gè)的供體。
• 通過更改replicates_per_patent參數(shù)乱投，可以為每個(gè)樣本創(chuàng)建幾個(gè)（n）偽重復(fù)咽笼；然后細(xì)胞被分成n個(gè)大小大致相等的子集。

NUM_OF_CELL_PER_DONOR = 30


def aggregate_and_filter(
    adata,
    cell_identity,
    donor_key="sample",
    condition_key="label",
    cell_identity_key="cell_type",
    obs_to_keep=[],  # which additional metadata to keep, e.g. gender, age, etc.
    replicates_per_patient=1,
):
    # subset adata to the given cell identity
    adata_cell_pop = adata[adata.obs[cell_identity_key] == cell_identity].copy()
    # check which donors to keep according to the number of cells specified with NUM_OF_CELL_PER_DONOR
    size_by_donor = adata_cell_pop.obs.groupby([donor_key]).size()
    donors_to_drop = [
        donor
        for donor in size_by_donor.index
        if size_by_donor[donor] <= NUM_OF_CELL_PER_DONOR
    ]
    if len(donors_to_drop) > 0:
        print("Dropping the following samples:")
        print(donors_to_drop)
    df = pd.DataFrame(columns=[*adata_cell_pop.var_names, *obs_to_keep])

    adata_cell_pop.obs[donor_key] = adata_cell_pop.obs[donor_key].astype("category")
    for i, donor in enumerate(donors := adata_cell_pop.obs[donor_key].cat.categories):
        print(f"\tProcessing donor {i+1} out of {len(donors)}...", end="\r")
        if donor not in donors_to_drop:
            adata_donor = adata_cell_pop[adata_cell_pop.obs[donor_key] == donor]
            # create replicates for each donor
            indices = list(adata_donor.obs_names)
            random.shuffle(indices)
            indices = np.array_split(np.array(indices), replicates_per_patient)
            for i, rep_idx in enumerate(indices):
                adata_replicate = adata_donor[rep_idx]
                # specify how to aggregate: sum gene expression for each gene for each donor and also keep the condition information
                agg_dict = {gene: "sum" for gene in adata_replicate.var_names}
                for obs in obs_to_keep:
                    agg_dict[obs] = "first"
                # create a df with all genes, donor and condition info
                df_donor = pd.DataFrame(adata_replicate.X.A)
                df_donor.index = adata_replicate.obs_names
                df_donor.columns = adata_replicate.var_names
                df_donor = df_donor.join(adata_replicate.obs[obs_to_keep])
                # aggregate
                df_donor = df_donor.groupby(donor_key).agg(agg_dict)
                df_donor[donor_key] = donor
                df.loc[f"donor_{donor}_{i}"] = df_donor.loc[donor]
    print("\n")
    # create AnnData object from the df
    adata_cell_pop = sc.AnnData(
        df[adata_cell_pop.var_names], obs=df.drop(columns=adata_cell_pop.var_names)
    )
    return adata_cell_pop

我們還需要定義一個(gè)單獨(dú)的函數(shù)來擬合edgeR GLM：

• fit_model將SingleCellExperiment對象作為輸入戚炫，創(chuàng)建設(shè)計(jì)矩陣并輸出擬合的GLM褐荷。我們還輸出DGEList類的edgeR對象來進(jìn)行一些探索性數(shù)據(jù)分析(EDA)。

%%R
fit_model <- function(adata_){
    # create an edgeR object with counts and grouping factor
    y <- DGEList(assay(adata_, "X"), group = colData(adata_)$label)
    # filter out genes with low counts
    print("Dimensions before subsetting:")
    print(dim(y))
    print("")
    keep <- filterByExpr(y)
    y <- y[keep, , keep.lib.sizes=FALSE]
    print("Dimensions after subsetting:")
    print(dim(y))
    print("")
    # normalize
    y <- calcNormFactors(y)
    # create a vector that is concatentation of condition and cell type that we will later use with contrasts
    group <- paste0(colData(adata_)$label, ".", colData(adata_)$cell_type)
    replicate <- colData(adata_)$replicate
    # create a design matrix: here we have multiple donors so also consider that in the design matrix
    design <- model.matrix(~ 0 + group + replicate)
    # estimate dispersion
    y <- estimateDisp(y, design = design)
    # fit the model
    fit <- glmQLFit(y, design)
    return(list("fit"=fit, "design"=design, "y"=y))
}

現(xiàn)在我們定義了所需的所有函數(shù)嘹悼，因此我們可以繼續(xù)創(chuàng)建偽批量叛甫。我們可能希望稍后查看可用的元數(shù)據(jù)层宫，因此將其保留在AnnData對象中。

obs_to_keep = ["label", "cell_type", "replicate", "sample"]

我們需要將原始計(jì)數(shù)傳遞給edgeR其监。因此萌腿，我們將.X設(shè)置為counts層，以確保為原始計(jì)數(shù)創(chuàng)建偽重復(fù)抖苦。

adata.X = adata.layers["counts"].copy()

接下來毁菱，我們使用pseudobulks創(chuàng)建AnnData對象。

# process first cell type separately...
cell_type = adata.obs["cell_type"].cat.categories[0]
print(
    f'Processing {cell_type} (1 out of {len(adata.obs["cell_type"].cat.categories)})...'
)
adata_pb = aggregate_and_filter(adata, cell_type, obs_to_keep=obs_to_keep)
for i, cell_type in enumerate(adata.obs["cell_type"].cat.categories[1:]):
    print(
        f'Processing {cell_type} ({i+2} out of {len(adata.obs["cell_type"].cat.categories)})...'
    )
    adata_cell_type = aggregate_and_filter(adata, cell_type, obs_to_keep=obs_to_keep)
    adata_pb = adata_pb.concatenate(adata_cell_type)

Processing B_cells (1 out of 8)...
Dropping the following samples:
['patient_1039_ctrl']
    Processing donor 16 out of 16...

Processing CD14_Monocytes (2 out of 8)...
    Processing donor 16 out of 16...

Processing CD4_T_cells (3 out of 8)...
    Processing donor 16 out of 16...

Processing CD8_T_cells (4 out of 8)...
Dropping the following samples:
['patient_101_ctrl', 'patient_1039_ctrl', 'patient_1039_stim', 'patient_107_ctrl', 'patient_107_stim', 'patient_1244_ctrl', 'patient_1244_stim']
    Processing donor 16 out of 16...

Processing Dendritic_cells (5 out of 8)...
Dropping the following samples:
['patient_1016_ctrl', 'patient_1016_stim', 'patient_101_ctrl', 'patient_1039_ctrl', 'patient_1039_stim', 'patient_107_ctrl', 'patient_107_stim']
    Processing donor 16 out of 16...

Processing FCGR3A_Monocytes (6 out of 8)...
Dropping the following samples:
['patient_1039_ctrl', 'patient_1039_stim', 'patient_107_ctrl', 'patient_107_stim', 'patient_1244_stim']
    Processing donor 16 out of 16...

Processing Megakaryocytes (7 out of 8)...
Dropping the following samples:
['patient_1015_ctrl', 'patient_1015_stim', 'patient_1016_ctrl', 'patient_101_ctrl', 'patient_1039_stim', 'patient_107_stim', 'patient_1244_ctrl', 'patient_1244_stim', 'patient_1256_ctrl', 'patient_1256_stim', 'patient_1488_ctrl']
    Processing donor 11 out of 11...

Processing NK_cells (8 out of 8)...
Dropping the following samples:
['patient_1039_ctrl', 'patient_1039_stim']
    Processing donor 16 out of 16...

差異基因表達(dá)結(jié)果的有效性很大程度上取決于統(tǒng)計(jì)模型中變異主軸的捕獲锌历。中間數(shù)據(jù)探索步驟税灌，例如對偽大量樣本的主成分分析（PCA）或多維標(biāo)度（MDS），可以識別變異源隐轩，從而可以指導(dǎo)構(gòu)建數(shù)據(jù)建模的相應(yīng)設(shè)計(jì)和對比矩陣图云。

如果未能考慮包括生物重復(fù)在內(nèi)的實(shí)驗(yàn)的生物變異性的多種來源，將會導(dǎo)致FDR膨脹卤材。雖然增加每個(gè)個(gè)體的細(xì)胞數(shù)量可以提高精度遮斥，但它對檢測個(gè)體差異的能力影響有限。因此扇丛，提高統(tǒng)計(jì)功效的最佳方法是增加獨(dú)立實(shí)驗(yàn)樣本的數(shù)量术吗。

由于我們的數(shù)據(jù)已經(jīng)生成，我們無法進(jìn)一步增加獨(dú)立實(shí)驗(yàn)樣本的數(shù)量帆精。盡管如此较屿，我們現(xiàn)在將探索我們的數(shù)據(jù)以確定變化的主軸，以正確生成我們的設(shè)計(jì)矩陣卓练。

我們對創(chuàng)建的偽重復(fù)執(zhí)行非沉吡停基本的EDA，以檢查某些患者/偽批量是否是我們需要排除的異常值昆庇，以免DE結(jié)果出現(xiàn)偏差末贾。我們將原始計(jì)數(shù)保存在counts層中，然后對計(jì)數(shù)進(jìn)行歸一化并計(jì)算歸一化偽批量計(jì)數(shù)的PCA坐標(biāo)整吆。

adata_pb.layers['counts'] = adata_pb.X.copy()

sc.pp.normalize_total(adata_pb, target_sum=1e6)
sc.pp.log1p(adata_pb)
sc.pp.pca(adata_pb)

接下來拱撵，我們通過所有可用的元數(shù)據(jù)查看在PCA圖上創(chuàng)建的偽重復(fù)和顏色，以查看是否存在我們可能希望包含在設(shè)計(jì)矩陣中的任何混雜因素表蝙。我們還添加了lib_size和log_lib_size列來檢查庫大小和PC組件之間是否存在相關(guān)性拴测。

adata_pb.obs["lib_size"] = np.sum(adata_pb.layers["counts"], axis=1)
adata_pb.obs["log_lib_size"] = np.log(adata_pb.obs["lib_size"])

sc.pl.pca(adata_pb, color=adata_pb.obs, ncols=1, size=300)

我們觀察PCA圖上細(xì)胞類型的分離以及刺激細(xì)胞和未刺激細(xì)胞的分離。其他協(xié)變量（批次）似乎與PCA分量沒有明顯相關(guān)府蛇，因此我們不將它們中的任何一個(gè)包含到我們的設(shè)計(jì)矩陣中集索。

如上所述，edgeR將原始計(jì)數(shù)作為輸入，因此我們在繼續(xù)之前將計(jì)數(shù)放回.X字段务荆。

adata_pb.X = adata_pb.layers['counts'].copy()

14.4.1 單組分析

首先妆距，我們展示如何準(zhǔn)備數(shù)據(jù)、構(gòu)建設(shè)計(jì)矩陣以及對一種特定細(xì)胞類型進(jìn)行DE測試函匕。

我們在CD14+單核細(xì)胞數(shù)據(jù)子集上運(yùn)行管道娱据，正如論文中所示，在該子群中未確定最多數(shù)量的DE基因盅惜。

adata_mono = adata_pb[adata_pb.obs["cell_type"] == "CD14_Monocytes"]
adata_mono

View of AnnData object with n_obs × n_vars = 16 × 15701
    obs: 'label', 'cell_type', 'replicate', 'sample', 'batch', 'lib_size', 'log_lib_size'
    uns: 'log1p', 'pca', 'label_colors', 'cell_type_colors', 'replicate_colors', 'sample_colors', 'batch_colors'
    obsm: 'X_pca'
    varm: 'PCs'
    layers: 'counts'

清理樣本名稱以使繪圖不那么擁擠中剩。

adata_mono.obs_names = [
    name.split("_")[2] + "_" + name.split("_")[3] for name in adata_mono.obs_names
]

%%time
%%R -i adata_mono
outs <-fit_model(adata_mono)

[1] "Dimensions before subsetting:"
[1] 15701    16
[1] ""
[1] "Dimensions after subsetting:"
[1] 3709   16
[1] ""
CPU times: user 2.95 s, sys: 224 ms, total: 3.18 s
Wall time: 4.27 s

%%R
fit <- outs$fit
y <- outs$y

由于我們在進(jìn)入分析時(shí)并未事先假設(shè)特定基因?qū)⒈簧险{(diào)或下調(diào)，因此我們需要可視化來理解DGE結(jié)果抒寂。MDS圖可以提供高層次的概述结啼。通常，我們期望不同條件下的樣品之間存在分離屈芜，如下圖所示我們的結(jié)果郊愧。

%%R
plotMDS(y, col=ifelse(y$samples$group == "stim", "red", "blue"))

生物變異系數(shù) (BCV) 圖顯示了生物群體中每個(gè)基因的平均變異性，作為平均表達(dá)的函數(shù)沸伏。例如，bcv為0.3表示各組間基因表達(dá)平均存在30%的變異性动分。

低豐度基因通常表現(xiàn)出較大的BCV毅糟，因?yàn)榈拓S度基因的讀數(shù)計(jì)數(shù)測量更不確定。另一方面澜公，平均表達(dá)量高的基因的定量更可靠姆另，因此它們通常具有較低的變異性，因此 BCV 也較低坟乾。使用此圖可以檢測異臣７基因或查明可能需要在設(shè)計(jì)矩陣中反映的其他實(shí)驗(yàn)因素。例如甚侣，出現(xiàn)在圖右上角的一組具有高平均表達(dá)和高BCV的基因可以標(biāo)記實(shí)驗(yàn)壓力明吩、污染等，特別是如果它們屬于相似的基因家族殷费。

在下面的 BCV 圖中印荔，我們看到一些具有高bcv的低豐度基因，但沒有看到具有高bcv的高豐度基因详羡，這表明不應(yīng)在設(shè)計(jì)矩陣中建模任何進(jìn)一步的實(shí)驗(yàn)考慮因素仍律。

%%R
plotBCV(y)

接下來，我們進(jìn)行擬似然檢驗(yàn)以找到對照條件和刺激條件之間的DE基因实柠。讓我們檢查一下我們的設(shè)計(jì)矩陣必須指定哪些列來指定要測試的正確列水泉。

%%R
colnames(y$design)

[1] "groupctrl.CD14_Monocytes" "groupstim.CD14_Monocytes"
[3] "replicatepatient_107"     "replicatepatient_1015"   
[5] "replicatepatient_1016"    "replicatepatient_1039"   
[7] "replicatepatient_1244"    "replicatepatient_1256"   
[9] "replicatepatient_1488"   

%%R -o tt
myContrast <- makeContrasts('groupstim.CD14_Monocytes-groupctrl.CD14_Monocytes', levels = y$design)
qlf <- glmQLFTest(fit, contrast=myContrast)
# get all of the DE genes and calculate Benjamini-Hochberg adjusted FDR
tt <- topTags(qlf, n = Inf)
tt <- tt$table

讓我們檢查一下這個(gè)表：對于每個(gè)沒有被edgeR過濾掉的基因（在我們的例子中是 3709），該表包含了DE測試的結(jié)果。該表可以保存為.csv文件以便稍后使用并用于可視化草则。在對完整數(shù)據(jù)進(jìn)行分析后钢拧，我們將在本節(jié)末尾展示如何使用火山圖。

tt.shape

(3709, 5)

tt[:5]

logFC   logCPM  F   PValue  FDR
HESX1   8.345536    6.773420    1281.013295 1.837373e-15    2.766927e-12
CD38    7.126846    7.420668    1243.793133 2.266164e-15    2.766927e-12
NT5C3A  5.657050    8.327003    1218.102628 2.628780e-15    2.766927e-12
SOCS1   4.388247    6.943768    1191.289806 3.079524e-15    2.766927e-12
GMPR    6.943484    7.031832    1159.601183 3.730018e-15    2.766927e-12

順便說一句畔师，人們還可以使用glmTreat測試相對于指定倍數(shù)變化所提供的系數(shù)或?qū)Ρ榷炔町惐磉_(dá)的基因娶靡。

%%R
tr <- glmTreat(fit, contrast=myContrast, lfc=1.5)
print(head(topTags(tr)))

Coefficient:  -1*groupctrl.CD14_Monocytes 1*groupstim.CD14_Monocytes 
          logFC unshrunk.logFC    logCPM       PValue          FDR
HESX1  8.345536       8.702486  6.773420 2.108737e-14 3.535749e-11
CD38   7.126846       7.219184  7.420668 2.292459e-14 3.535749e-11
NT5C3A 5.657050       5.674640  8.327003 3.149544e-14 3.535749e-11
IL1RN  6.588583       6.596787 10.358946 4.249159e-14 3.535749e-11
GMPR   6.943484       7.047504  7.031832 4.766446e-14 3.535749e-11
DEFB1  6.654201       6.696759  8.067159 7.670147e-14 4.741429e-11

最后，我們可以看到有多少基因的FDR校正值小于0.01看锉。

涂片圖（也稱為MA姿锭、M值與A值圖）顯示基因的對數(shù)倍數(shù)變化與其平均豐度的函數(shù)關(guān)系。我們通常在低豐度范圍內(nèi)觀察到較高的logFC伯铣，因?yàn)樽x數(shù)計(jì)數(shù)在低豐度范圍內(nèi)變化更大呻此，導(dǎo)致logFC估計(jì)值較大。如果我們將loess曲線擬合到logFC和平均logCPM值腔寡，則趨勢應(yīng)以零為中心焚鲜。任何偏離此值的情況都可能表明數(shù)據(jù)尚未正確標(biāo)準(zhǔn)化。平均表達(dá)量大且絕對值 logFC 大的基因可以標(biāo)記生物學(xué)上感興趣的基因以供研究和隨訪放前。

%%R
plotSmear(qlf, de.tags = rownames(tt)[which(tt$FDR<0.01)])

14.4.2 多組分析

接下來忿磅，我們展示如何準(zhǔn)備數(shù)據(jù)、構(gòu)建設(shè)計(jì)矩陣并使用所有可用細(xì)胞類型的完整數(shù)據(jù)集的對比來執(zhí)行DE測試凭语。

%%time
%%R -i adata_pb
outs <-fit_model(adata_pb)

[1] "Dimensions before subsetting:"
[1] 15701    90
[1] ""
[1] "Dimensions after subsetting:"
[1] 2358   90
[1] ""
CPU times: user 11.9 s, sys: 39.7 ms, total: 12 s
Wall time: 12 s

%%R
fit <- outs$fit
y <- outs$y

現(xiàn)在我們使用對比對每種細(xì)胞類型進(jìn)行擬似然檢驗(yàn)葱她。由于沒有直接的方法將表從R循環(huán)內(nèi)的R移動(dòng)到Python，因此我們手動(dòng)獲取每種細(xì)胞類型的結(jié)果似扔。

%%R -i adata_pb -o de_per_cell_type
de_per_cell_type <- list()
for (cell_type in unique(colData(adata_pb)$cell_type)) {
    print(cell_type)
    # create contrast for this cell type
    myContrast <- makeContrasts(paste0("groupstim.", cell_type, "-groupctrl.", cell_type), levels = y$design)
    # perform QLF test
    qlf <- glmQLFTest(fit, contrast=myContrast)
    # get all of the DE genes and calculate Benjamini-Hochberg adjusted FDR
    tt <- topTags(qlf, n = Inf)
    # save in the list with the results for all the cell types
    de_per_cell_type[[cell_type]] <- tt$table
}

[1] "B_cells"
[1] "CD14_Monocytes"
[1] "CD4_T_cells"
[1] "CD8_T_cells"
[1] "Dendritic_cells"
[1] "FCGR3A_Monocytes"
[1] "NK_cells"

現(xiàn)在讓我們使用sc_toolbox包 (https://github.com/schillerlab/sc-toolbox) 中的 sc_toolbox.tools.de_res_to_anndata將其保存在原始adata對象的.uns中吨些，這會保存結(jié)果，就像使用scanpy創(chuàng)建的rank_genes_groups()函數(shù)一樣炒辉。像這樣存儲DE表的優(yōu)點(diǎn)是我們現(xiàn)在還可以使用標(biāo)準(zhǔn)的scanpy繪圖函數(shù)豪墅。我們將每種細(xì)胞類型的DEG表保存為.csv，因?yàn)槲覀兩院笤诩?xì)胞間通訊章節(jié)中將需要它們黔寇。

# get cell types that we ran the analysis for
cell_types = de_per_cell_type.keys()
# add the table to .uns for each cell type
for cell_type in cell_types:
    df = de_per_cell_type[cell_type]
    df["gene_symbol"] = df.index
    df["cell_type"] = cell_type
    sc_toolbox.tools.de_res_to_anndata(
        adata,
        df,
        groupby="cell_type",
        score_col="logCPM",
        pval_col="PValue",
        pval_adj_col="FDR",
        lfc_col="logFC",
        key_added="edgeR_" + cell_type,
    )
    df.to_csv(f"de_edgeR_{cell_type}.csv")

為了再次將表作為pandasDataFrame獲取偶器，我們使用sc.get.rank_genes_groups_df函數(shù)。

sc.get.rank_genes_groups_df(adata, group="CD14_Monocytes", key="edgeR_CD14_Monocytes")[
    :5
]

    names   scores  logfoldchanges  pvals   pvals_adj
0   FTH1    16.186098   -0.572272   0.001589    0.002712
1   MALAT1  16.025682   0.010049    0.877078    0.897245
2   B2M 15.468524   0.677266    0.0 0.0
3   TMSB4X  15.054697   -0.217131   0.004663    0.007331
4   FTL 14.984937   -0.041832   0.669951    0.710956

14.4.3 關(guān)于edgeR的tips

• 需要原始計(jì)數(shù)作為輸入
• 需要來自單細(xì)胞實(shí)驗(yàn)的pseudobulks偽重復(fù)
• 如果單細(xì)胞實(shí)驗(yàn)中有多個(gè)供體缝裤，并且用戶想要考慮患者的變異性状囱，我們建議為每個(gè)患者創(chuàng)建2或3個(gè)偽重復(fù)，并將患者信息納入設(shè)計(jì)矩陣中

14.5 單細(xì)胞特異算法

MAST使用兩部分廣義線性模型對單細(xì)胞基因表達(dá)進(jìn)行建模倘是。一部分模擬細(xì)胞中每個(gè)基因的離散表達(dá)率亭枷，而另一部分模擬條件連續(xù)表達(dá)水平。

MAST框架使用兩部分廣義線性模型對單細(xì)胞基因表達(dá)進(jìn)行建模搀崭。MAST的一個(gè)組件模擬細(xì)胞中每個(gè)基因的離散表達(dá)率叨粘，而另一個(gè)組件則模擬條件連續(xù)表達(dá)水平（以所表達(dá)的基因?yàn)闂l件）猾编。 欲了解更多詳細(xì)信息，請閱讀原始出版物[Finak et al., 2015]升敲。

MAST將標(biāo)準(zhǔn)化計(jì)數(shù)作為輸入答倡，因此我們首先采用“counts”層，然后執(zhí)行標(biāo)準(zhǔn)化步驟驴党。

adata.X = adata.layers["counts"].copy()
sc.pp.normalize_total(adata, target_sum=1e6)
sc.pp.log1p(adata)

我們定義了一個(gè)小輔助函數(shù)來處理R和Python之間的一些對象類型轉(zhuǎn)換問題瘪撇。我們還需要過濾掉每個(gè)亞群的少量細(xì)胞（本例中為3個(gè)）中表達(dá)的基因，因?yàn)槟Ｐ托枰軌蚬烙?jì)每個(gè)基因的方差港庄。

def prep_anndata(adata_):
    def fix_dtypes(adata_):
        df = pd.DataFrame(adata_.X.A, index=adata_.obs_names, columns=adata_.var_names)
        df = df.join(adata_.obs)
        return sc.AnnData(df[adata_.var_names], obs=df.drop(columns=adata_.var_names))

    adata_ = fix_dtypes(adata_)
    sc.pp.filter_genes(adata_, min_cells=3)
    return adata_

14.5.1 單組分析

與edgeR一樣倔既，可以使用完整數(shù)據(jù)集擬合模型，然后使用對比對每種感興趣的細(xì)胞類型執(zhí)行測試鹏氧，但在本篇中渤涌，我們僅顯示一種細(xì)胞類型（即CD14單核細(xì)胞）的MAST-RE管道，縮短運(yùn)行時(shí)間把还。

adata_mono = adata[adata.obs["cell_type"] == "CD14_Monocytes"].copy()
adata_mono

AnnData object with n_obs × n_vars = 5696 × 15701
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'n_genes', 'sample'
    var: 'name', 'n_cells'
    uns: 'edgeR_B_cells', 'edgeR_CD14_Monocytes', 'edgeR_CD4_T_cells', 'edgeR_CD8_T_cells', 'edgeR_Dendritic_cells', 'edgeR_FCGR3A_Monocytes', 'edgeR_NK_cells', 'log1p'
    obsm: 'X_pca', 'X_umap'
    layers: 'counts'

sc.pp.filter_genes(adata_mono, min_cells=3)
adata_mono

AnnData object with n_obs × n_vars = 5696 × 12268
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'n_genes', 'sample'
    var: 'name', 'n_cells'
    uns: 'edgeR_B_cells', 'edgeR_CD14_Monocytes', 'edgeR_CD4_T_cells', 'edgeR_CD8_T_cells', 'edgeR_Dendritic_cells', 'edgeR_FCGR3A_Monocytes', 'edgeR_NK_cells', 'log1p'
    obsm: 'X_pca', 'X_umap'
    layers: 'counts'

接下來我們?nèi)缟纤鲞^濾這兩個(gè)對象实蓬。

adata_mono = prep_anndata(adata_mono)
adata_mono

AnnData object with n_obs × n_vars = 5696 × 12268
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'n_genes', 'sample'
    var: 'n_cells'

與執(zhí)行多個(gè)統(tǒng)計(jì)測試時(shí)的情況一樣，必須使用Benjamini-Hochberg校正等方法對條件下DGE測試獲得的p值進(jìn)行校正吊履。

與edgeR分析類似安皱，我們定義一個(gè)用于分析的單獨(dú)函數(shù)：

• find_de_MAST_RE采用SingleCellExperiment對象作為輸入，并使用RE管道運(yùn)行 MAST艇炎。該函數(shù)的輸出是表（Python中的pandas DataFrame）酌伊，其中包含分析結(jié)果，例如每個(gè)基因的對數(shù)倍數(shù)變化冕臭、p值和FDR校正值腺晾。

adata_mono.obs["cell_type"] = [
    ct.replace(" ", "_") for ct in adata_mono.obs["cell_type"]
]
adata_mono.obs["cell_type"] = [
    ct.replace("+", "") for ct in adata_mono.obs["cell_type"]
]

%%R
find_de_MAST_RE <- function(adata_){
    # create a MAST object
    sca <- SceToSingleCellAssay(adata_, class = "SingleCellAssay")
    print("Dimensions before subsetting:")
    print(dim(sca))
    print("")
    # keep genes that are expressed in more than 10% of all cells
    sca <- sca[freq(sca)>0.1,]
    print("Dimensions after subsetting:")
    print(dim(sca))
    print("")
    # add a column to the data which contains scaled number of genes that are expressed in each cell
    cdr2 <- colSums(assay(sca)>0)
    colData(sca)$ngeneson <- scale(cdr2)
    # store the columns that we are interested in as factors
    label <- factor(colData(sca)$label)
    # set the reference level
    label <- relevel(label,"ctrl")
    colData(sca)$label <- label
    celltype <- factor(colData(sca)$cell_type)
    colData(sca)$celltype <- celltype
    # same for donors (which we need to model random effects)
    replicate <- factor(colData(sca)$replicate)
    colData(sca)$replicate <- replicate
    # create a group per condition-celltype combination
    colData(sca)$group <- paste0(colData(adata_)$label, ".", colData(adata_)$cell_type)
    colData(sca)$group <- factor(colData(sca)$group)
    # define and fit the model
    zlmCond <- zlm(formula = ~ngeneson + group + (1 | replicate), 
                   sca=sca, 
                   method='glmer', 
                   ebayes=F, 
                   strictConvergence=F,
                   fitArgsD=list(nAGQ = 0)) # to speed up calculations
    
    # perform likelihood-ratio test for the condition that we are interested in    
    summaryCond <- summary(zlmCond, doLRT='groupstim.CD14_Monocytes')
    # get the table with log-fold changes and p-values
    summaryDt <- summaryCond$datatable
    result <- merge(summaryDt[contrast=='groupstim.CD14_Monocytes' & component=='H',.(primerid, `Pr(>Chisq)`)], # p-values
                     summaryDt[contrast=='groupstim.CD14_Monocytes' & component=='logFC', .(primerid, coef)],
                     by='primerid') # logFC coefficients
    # MAST uses natural logarithm so we convert the coefficients to log2 base to be comparable to edgeR
    result[,coef:=result[,coef]/log(2)]
    # do multiple testing correction
    result[,FDR:=p.adjust(`Pr(>Chisq)`, 'fdr')]
    result = result[result$FDR<0.01,, drop=F]

    result <- stats::na.omit(as.data.frame(result))
    return(result)
}

我們運(yùn)行單核細(xì)胞的分析流程燕锥。

%%time
%%R -i adata_mono -o res
res <-find_de_MAST_RE(adata_mono)

[1] "Dimensions before subsetting:"
[1] 12268  5696
[1] ""
[1] "Dimensions after subsetting:"
[1] 1676 5696
[1] ""
CPU times: user 9min 35s, sys: 2.96 s, total: 9min 38s
Wall time: 9min 43s

看一下運(yùn)行結(jié)果辜贵。

res[:5]

    primerid    Pr(>Chisq)  coef    FDR
1   AAED1   8.410341e-19    0.709836    1.597106e-18
2   ABI1    1.401688e-05    0.312797    1.824921e-05
3   ABRACL  3.920183e-06    0.418028    5.201004e-06
4   ACADVL  2.121110e-26    -0.751305   4.656979e-26
5   ACOT9   3.424174e-151   2.921133    2.689503e-150

我們像以前一樣將結(jié)果存儲在.uns中。請注意归形，我們不需要score列托慨，因此我們只需將logFC傳遞到那里即可。

res["gene_symbol"] = res["primerid"]
res["cell_type"] = "CD14_Monocytes"
sc_toolbox.tools.de_res_to_anndata(
    adata,
    res,
    groupby="cell_type",
    score_col="coef",
    pval_col="Pr(>Chisq)",
    pval_adj_col="FDR",
    lfc_col="coef",
    key_added="MAST_CD14_Monocytes",
)

adata_copy = adata.copy()

14.6 可視化

我們可以使用熱圖和火山圖來可視化結(jié)果暇榴。熱圖的每一行對應(yīng)一個(gè)基因厚棵，每一列對應(yīng)一個(gè)單細(xì)胞。顏色越亮蔼紧，該基因在特定細(xì)胞中的表達(dá)就越高婆硬。由于我們只繪制DE基因，因此我們希望看到兩種條件之間表達(dá)的明顯差異奸例”蚍福火山圖通常用于可視化統(tǒng)計(jì)測試的結(jié)果向楼，它們在x軸上顯示表達(dá)變化（對數(shù)倍數(shù)變化），在y軸上顯示統(tǒng)計(jì)顯著性（FDR校正的p值）谐区。我們對FDR校正p值低于0.01且對數(shù)倍變化超過1.5的基因進(jìn)行顏色編碼湖蜕。

我們在繪制熱圖之前對數(shù)據(jù)進(jìn)行歸一化，以便更好地查看兩種條件之間表達(dá)的差異宋列。

adata.X = adata.layers["counts"].copy()
sc.pp.normalize_total(adata, target_sum=1e6)
sc.pp.log1p(adata)

接下來昭抒，我們?yōu)闊釄D定義一個(gè)輔助繪圖函數(shù)。

FDR = 0.01
LOG_FOLD_CHANGE = 1.5


def plot_heatmap(adata, group_key, group_name="cell_type", groupby="label"):
    cell_type = "_".join(group_key.split("_")[1:])
    res = sc.get.rank_genes_groups_df(adata, group=cell_type, key=group_key)
    res.index = res["names"].values
    res = res[
        (res["pvals_adj"] < FDR) & (abs(res["logfoldchanges"]) > LOG_FOLD_CHANGE)
    ].sort_values(by=["logfoldchanges"])
    print(f"Plotting {len(res)} genes...")
    markers = list(res.index)
    sc.pl.heatmap(
        adata[adata.obs[group_name] == cell_type].copy(),
        markers,
        groupby=groupby,
        swap_axes=True,
    )

最后炼杖，我們可以使用RE繪制edgeR和MAST的CD14+單核細(xì)胞的DEG熱圖灭返。

plot_heatmap(adata, "edgeR_CD14_Monocytes")

Plotting 303 genes...

plot_heatmap(adata, "MAST_CD14_Monocytes")

Plotting 436 genes...

我們觀察到，MAST根據(jù)我們給定的調(diào)整p值和對數(shù)倍變化的截止值識別出436個(gè)差異基因嘹叫，而edgeR識別出303個(gè)差異基因婆殿。

接下來，我們定義火山圖的輔助繪圖函數(shù)罩扇。

FDR = 0.01
LOG_FOLD_CHANGE = 1.5


def volcano_plot(adata, group_key, group_name="cell_type", groupby="label", title=None):
    cell_type = "_".join(group_key.split("_")[1:])
    result = sc.get.rank_genes_groups_df(adata, group=cell_type, key=group_key).copy()
    result["-logQ"] = -np.log(result["pvals"].astype("float"))
    lowqval_de = result.loc[abs(result["logfoldchanges"]) > LOG_FOLD_CHANGE]
    other_de = result.loc[abs(result["logfoldchanges"]) <= LOG_FOLD_CHANGE]

    fig, ax = plt.subplots()
    sns.regplot(
        x=other_de["logfoldchanges"],
        y=other_de["-logQ"],
        fit_reg=False,
        scatter_kws={"s": 6},
    )
    sns.regplot(
        x=lowqval_de["logfoldchanges"],
        y=lowqval_de["-logQ"],
        fit_reg=False,
        scatter_kws={"s": 6},
    )
    ax.set_xlabel("log2 FC")
    ax.set_ylabel("-log Q-value")

    if title is None:
        title = group_key.replace("_", " ")
    plt.title(title)
    plt.show()

volcano_plot(adata, "MAST_CD14_Monocytes")

volcano_plot(adata, "edgeR_CD14_Monocytes")

從熱圖婆芦，尤其是火山圖中，我們可以看到喂饥，edgeR識別出上調(diào)基因多于下調(diào)基因（刺激與對照）消约，而MAST識別出相似數(shù)量的上調(diào)和下調(diào)基因。