人妻福利在线视频,AV久艹在线

前言提到，在過去兩天的教程中，我們講解了使用omicverse進(jìn)行單細(xì)胞測序數(shù)據(jù)的質(zhì)控以及歸一化的一些思想。關(guān)于omicverse的使用文檔與安裝教程可以參考我們的readthedocs.

就是，本系列教程是我?guī)П究粕玫降?，所以概念?huì)盡可能地通俗，詳細(xì)，但對于急于求成的人，可能不是一個(gè)很好的教程。

1. 背景

我們現(xiàn)在有了一個(gè)經(jīng)過標(biāo)準(zhǔn)化的單細(xì)胞測序數(shù)據(jù)，其保留了細(xì)胞差異，同時(shí)消除了技術(shù)帶來的采樣誤差。通常來說，我們經(jīng)過上游的mapping（比對），我們會(huì)得到30,000到50,000個(gè)不等的基因。但是，我們只在第一步的質(zhì)控中，刪除了一小部分基因（少于20個(gè)細(xì)胞表達(dá)的基因）。但實(shí)際上，一個(gè)細(xì)胞表達(dá)的基因大約是3,000個(gè)左右。意味著我們的測序數(shù)據(jù)中的一大部分基因可能不具有與單細(xì)胞數(shù)據(jù)相關(guān)的生物學(xué)意義，這些基因大多數(shù)包含0計(jì)數(shù)，或者是在許多細(xì)胞中普遍出現(xiàn)的基因。

故在預(yù)處理環(huán)節(jié)，我們需要計(jì)算高可變基因（特征基因），排除那些不具有分析意義的基因，避免影響我們下游的建模以及生物學(xué)檢測。

特征選擇

在本章中，我們將介紹三種不同的特征基因選擇：基于基因離散度，基于基因歸一化方差以及基于基因的皮爾森殘差。在傳統(tǒng)的分析流程中，我們會(huì)采用基于基因離散度的方式去計(jì)算高可變基因，一般來說，我們首先確定了單細(xì)胞數(shù)據(jù)集中變異最大的一組基因。我們計(jì)算了所有單細(xì)胞中每個(gè)基因的平均值和離散度（方差/平均值），并根據(jù)基因的平均表達(dá)將其分為 20 個(gè)箱。然后，在每個(gè)箱內(nèi)，我們對箱內(nèi)所有基因的離散度進(jìn)行z歸一化，以識別表達(dá)值高度可變的基因。

2. 基于基因離散度

我們首先設(shè)置我們的環(huán)境。

import omicverse as ov
import scanpy as sc

ov.ov_plot_set()

接下來，我們加載此前分析所用到的數(shù)據(jù)集，并進(jìn)行移位對數(shù)標(biāo)準(zhǔn)化處理。

adata = sc.read(
    filename="s4d8_quality_control.h5ad",
    #backup_url="https://figshare.com/ndownloader/files/40014331",
)
#存儲(chǔ)原始數(shù)據(jù)以便后續(xù)還原
ov.utils.store_layers(adata,layers='counts')
adata.layers['counts']=adata.X.copy()

sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
adata

......The X of adata have been stored in counts
normalizing counts per cell
    finished (0:00:00)
AnnData object with n_obs × n_vars = 14814 × 20171
    obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'outlier', 'mt_outlier', 'scDblFinder_score', 'scDblFinder_class'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
    uns: 'layers_counts', 'log1p'
    layers: 'counts', 'soupX_counts'

接下來，我們調(diào)用scanpy包里的pp.highly_variable_genes函數(shù)來計(jì)算高可變基因，由于我們使用的是基于基因離散度的方法，故我們需要設(shè)置flavor='seurat'，該方法也是默認(rèn)方法?；诨螂x散度的方法尋找高變基因有兩個(gè)途徑：

指定目的高變基因數(shù)
指定離散度，平均值過濾高變基因

在這里，我們分別嘗試兩種不同的方法，首先是基于指定高變基因數(shù)的方法：

2.1 指定高可變基因數(shù)

由于是教程演示，所以我們不希望改變原來的anndata對象的.var，故我們設(shè)置inplace=False，旨在觀察和對比高可變基因

adata_dis_num=sc.pp.highly_variable_genes(
    adata,
    flavor="seurat",
    n_top_genes=2000,
    subset=False,
    inplace=False,
)
adata_dis_num

If you pass `n_top_genes`, all cutoffs are ignored.
extracting highly variable genes
    finished (0:00:00)

adata_dis_num['highly_variable'].value_counts()

False    18171
True      2000
Name: highly_variable, dtype: int64

我們發(fā)現(xiàn)，一共選擇了2000個(gè)高可變基因，這與我們最開始的分析目標(biāo)一致，至于選擇出來的高可變基因有什么用途，則留在下一章節(jié)討論

2.2 指定基因離散度與平均值閾值

除了指定高可變基因數(shù)量外，我們還可以通過基因的離散度與平均值的閾值，去計(jì)算高可變基因，但是，這種啟發(fā)式的方法是數(shù)據(jù)敏感的，閾值所計(jì)算出來的高可變基因可能并不是全部高可變。故基于閾值的方法在現(xiàn)在的分析中應(yīng)用較少。

adata_dis_cutoff=sc.pp.highly_variable_genes(
    adata,
    flavor="seurat",
    min_disp=0.5,
    min_mean=0.0125,
    max_mean=3,
    subset=False,
    inplace=False,
)
adata_dis_cutoff['highly_variable'].value_counts()

extracting highly variable genes
    finished (0:00:00)
False    17291
True      2880
Name: highly_variable, dtype: int64

通過簡單的計(jì)算我們發(fā)現(xiàn)，一共有2,880個(gè)高可變基因被選擇了出來，一般2000-3000內(nèi)的高可變基因數(shù)都是能接受的

因?yàn)榉歉呖勺兓蛟谙掠畏治鰰r(shí)會(huì)被過濾掉，我們會(huì)將歸一化值給保存進(jìn)raw文件，但這有個(gè)缺陷，我們將會(huì)永遠(yuǎn)失去非高可變基因的原始count。因此可以使用釋放函數(shù) ov.utils.retrieve_layers，與存放在adata.layers中的數(shù)據(jù)不同，retrieve_layers函數(shù)還可還原全基因的原始計(jì)數(shù)，我們以下將進(jìn)行一個(gè)簡單的對比實(shí)驗(yàn)。

#計(jì)算高可變基因（覆蓋，inplace=True）
adata.raw = adata
sc.pp.highly_variable_genes(
    adata,
    flavor="seurat",
)
#只保留高可變基因
adata = adata[:, adata.var.highly_variable]
print('shape: ',adata.shape)

adata_counts=adata.copy()
ov.utils.retrieve_layers(adata_counts,layers='counts')
print('normalize adata:',adata.X.max())
print('raw count adata:',adata_counts.X.max())
print('shape: ',adata_counts.shape)

extracting highly variable genes
    finished (0:00:00)
--> added
    'highly_variable', boolean vector (adata.var)
    'means', float vector (adata.var)
    'dispersions', float vector (adata.var)
    'dispersions_norm', float vector (adata.var)
shape:  (14814, 2880)
......The X of adata have been stored in raw
......The layers counts of adata have been retreved
normalize adata: 4.9448869840604806
raw count adata: 694.0
shape:  (14814, 2880)

#釋放所有基因
adata_counts=adata.raw.to_adata().copy()
ov.utils.retrieve_layers(adata_counts,layers='counts')
print('normalize adata:',adata.X.max())
print('raw count adata:',adata_counts.X.max())
print('shape: ',adata_counts.shape)

......The X of adata have been stored in raw
......The layers counts of adata have been retreved
normalize adata: 4.9448869840604806
raw count adata: 889.0
shape:  (14814, 20171)

可以發(fā)現(xiàn)，我們使用.raw還原基因名后，retrieve_layers即可自動(dòng)提取.raw中所包含的counts值。即數(shù)據(jù)的形狀為 (14814, 20171)

3. 基于基因歸一化方差

在過去的單細(xì)胞研究中人們發(fā)現(xiàn)，僅根據(jù)對數(shù)歸一化法，并使用離散度的高可變基因選擇方法無法解釋單細(xì)胞 RNA-seq 固有的均值-方差關(guān)系。故在seurat v3中，研究人員使用了方差穩(wěn)定變換來糾正這一點(diǎn)誤差。這意味著我們將不使用標(biāo)準(zhǔn)化后的數(shù)據(jù)來計(jì)算高可變基因。該方法的計(jì)算步驟如下：

我們首先計(jì)算每一個(gè)基因的平均值 $\overline{x_i}$ 與方差 $\sigma_i$ ，然后分別對平均值與方差進(jìn)行l(wèi)og對數(shù)化。然后我們使用2次多項(xiàng)式，將方差作為均值的函數(shù)，進(jìn)行局部多項(xiàng)式回歸。

$\sigma (x) = a x^2 +b x +c$

通過上述公式，我們可以獲得每一個(gè)基因的預(yù)期方差，然后我們進(jìn)行z轉(zhuǎn)換：

$z_{ij}= \frac {x_{ij}-\overline{x_i}} {\sigma (x_i)}$

其中 $z_{ij}$ 是細(xì)胞j中基因i的標(biāo)準(zhǔn)化值， $x_{ij}$ 是細(xì)胞j中基因i的原始值， $\overline{x_i}$ 是基因i的平均原始值，并且 $\sigma (x_i)$ 是從全局均值方差擬合得出的特征i的預(yù)期標(biāo)準(zhǔn)差。為了減少技術(shù)異常值的影響，我們將標(biāo)準(zhǔn)化值削減為最大值 $\sqrt {N}$ ，其中N是細(xì)胞總數(shù)。然后，對于每個(gè)基因，我們計(jì)算所有細(xì)胞標(biāo)準(zhǔn)化值的方差。該方差表示控制均值表達(dá)后單細(xì)胞離散度的度量，我們直接使用它對特征進(jìn)行排序。

在scanpy中，我們需要設(shè)定flavor="seurat_v3"以選擇基于基因歸一化方差的方法，并指定計(jì)數(shù)矩陣為未歸一化的矩陣，即layer='counts'，同時(shí)我們選擇標(biāo)準(zhǔn)化方差最高的 2,000 個(gè)基因作為“高度可變”。

adata_var_num=sc.pp.highly_variable_genes(
    adata,
    flavor="seurat_v3",
    layer='counts',
    n_top_genes=2000,
    subset=False,
    inplace=False,
)
adata_var_num['highly_variable'].value_counts()

If you pass `n_top_genes`, all cutoffs are ignored.
extracting highly variable genes
False    18171
True      2000
Name: highly_variable, dtype: int64

4. 基于基因皮爾森近似殘差

除了基于基因歸一化殘差以及離散度的方法外，我們還可以基于皮爾森近似殘差的思想來計(jì)算高可變基因，你一定注意到了我們此前的歸一化步驟中包含了皮爾森殘差的方法，但是在這里，我們依然是使用原始計(jì)數(shù)來計(jì)算高可變基因，這是scanpy的實(shí)現(xiàn)原理所導(dǎo)致的，即我們可以使用移位對數(shù)來獲得標(biāo)準(zhǔn)化數(shù)據(jù)，同時(shí)使用皮爾森殘差法來獲得高可變基因

針對沒有生物學(xué)變異的UMI或計(jì)數(shù)數(shù)據(jù)的常見建模假設(shè)是，每個(gè)基因 $g$ 在細(xì)胞 $c$ 中占據(jù)了總計(jì)數(shù) $n_c$ 的一定比例 $p_g$ [4, 6-10]。然后，觀察到的UMI計(jì)數(shù) $X_{c g}$ 被建模為泊松分布或負(fù)二項(xiàng)分布（NB），其期望值為 $\mu_{c g}=p_g n_c$ ，不考慮零膨脹：

$\begin{array}{r} X_{c g} \sim \text { Poisson }\left(\mu_{c g}\right) \text { or } \mathrm{NB}\left(\mu_{c g}, \theta\right), \\ \mu_{c g}=n_c p_g . \end{array}$

泊松模型有一個(gè)最大似然解，可以寫成閉式形式： $\hat{n}_c=\sum_g X_{c g}$ （測序深度）， $\hat{p}_g=\sum_c X_{c g} / \sum_c \hat{n}_c$ ，或者合并在一起：

$\hat{\mu}_{c g}=\frac{\sum_j X_{c j} \cdot \sum_i X_{i g}}{\sum_{i j} X_{i j}}$

對于負(fù)二項(xiàng)模型，這個(gè)解僅適用于近似情況。利用這個(gè)解，皮爾森近似殘差可表示為：

$Z_{c g}=\frac{X_{c g}-\hat{\mu}_{c g}}{\sqrt{\hat{\mu}_{c g}+\hat{\mu}_{c g}^2 / \theta}},$

其中 $\mu_{c g}+\mu_{c g}^2 / \theta$ 是NB的方差， $\theta \rightarrow \infty$ 時(shí)得到泊松極限。皮爾森近似殘差的方差，除了一個(gè)常數(shù)外，等于皮爾森 $\chi^2$ 擬合優(yōu)度統(tǒng)計(jì)量，用來量化每個(gè)基因與這個(gè)常量表達(dá)模型的偏差程度。正如Aedin Culhane指出的那樣，在泊松模型下對皮爾森近似殘差進(jìn)行奇異值分解被稱為對應(yīng)分析，這是一種歷史悠久的方法。

Hafemeister和Satija建議在高變異基因（HVG）篩選中使用來自相關(guān)NB回歸模型的皮爾森近似殘差，也可作為下游處理的數(shù)據(jù)轉(zhuǎn)換方法。與此同時(shí)，Townes等人建議使用來自上述泊松模型的偏差殘差進(jìn)行HVG篩選

在scanpy中，我們將使用experimental.pp.highly_variable_genes來計(jì)算基于皮爾森殘差的高可變基因

adata_pearson_num=sc.experimental.pp.highly_variable_genes(
    adata, 
    flavor="pearson_residuals",
    layer='counts',
    n_top_genes=2000,
    subset=False,
    inplace=False,
)
adata_pearson_num['highly_variable'].value_counts()

extracting highly variable genes
False    18171
True      2000
Name: highly_variable, dtype: int64

adata_dis_num

現(xiàn)在，我們已經(jīng)獲得了3種不同的特征基因選擇方法所得到的高可變基因，我們可以對比一下基因之間的重疊率，至于功能的研究則不在本章的考慮范圍內(nèi)。我們先將基因名賦予此前的結(jié)果

adata_dis_num.index=adata.var_names.copy()
adata_dis_cutoff.index=adata.var_names.copy()
adata_var_num.index=adata.var_names.copy()
adata_pearson_num.index=adata.var_names.copy()

import matplotlib.pyplot as plt
from matplotlib_venn import venn3

# 三個(gè)列表的元素
list1 = set(adata_dis_num.loc[adata_dis_num['highly_variable']==True].index.tolist())
list2 = set(adata_var_num.loc[adata_var_num['highly_variable']==True].index.tolist())
list3 = set(adata_pearson_num.loc[adata_pearson_num['highly_variable']==True].index.tolist())

# 繪制 Venn 圖
venn = venn3([list1, list2, list3], set_labels=('Dis', 'Var', 'Pearson'))


# 顯示圖形
plt.title("Venn Diagram of Three HVGs")
plt.show()

三種不同的HVGs的Venn圖

有意思的事情出現(xiàn)了，我們發(fā)現(xiàn)三種不同方法所找到的高可變基因（HVGs）僅有656個(gè)是相同的，這意味著不同的方法所尋找到的高可變基因可能會(huì)影響你的下游分析的誤差。在我們的分析中，我們推薦使用皮爾森殘差法來獲得高可變基因。

在omicverse中，我們可以直接使用ov.pp.preprocess完成預(yù)處理步驟，需要注意的是，當(dāng)omicverse的版本小于1.4.13時(shí)，mode的參數(shù)只能設(shè)置為scanpy或pearson，而在更高的版本中，normalize|HVGs：我們使用 | 來控制預(yù)處理步驟，| 前用于歸一化步驟，可以是 shiftlog 或 pearson，| 后用于高度可變基因計(jì)算步驟，可以是 pearson 或 seurat。我們的默認(rèn)值是 shiftlog|pearson。

adata = sc.read(
    filename="s4d8_quality_control.h5ad",
    #backup_url="https://figshare.com/ndownloader/files/40014331",
)
#存儲(chǔ)原始數(shù)據(jù)以便后續(xù)還原
ov.utils.store_layers(adata,layers='counts')
adata.layers['counts']=adata.X.copy()

adata=ov.pp.preprocess(adata,
                       mode='shiftlog|pearson',
                       n_HVGs=2000,)
adata

......The X of adata have been stored in counts
Begin robust gene identification
After filtration, 20171/20171 genes are kept. Among 20171 genes, 20171 genes are robust.
End of robust gene identification.
Begin size normalization: shiftlog and HVGs selection pearson
normalizing counts per cell The following highly-expressed genes are not considered during normalization factor computation:
[]
    finished (0:00:00)
extracting highly variable genes
--> added
    'highly_variable', boolean vector (adata.var)
    'highly_variable_rank', float vector (adata.var)
    'highly_variable_nbatches', int vector (adata.var)
    'highly_variable_intersection', boolean vector (adata.var)
    'means', float vector (adata.var)
    'variances', float vector (adata.var)
    'residual_variances', float vector (adata.var)
End of size normalization: shiftlog and HVGs selection pearson



AnnData object with n_obs × n_vars = 14814 × 20171
    obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'outlier', 'mt_outlier', 'scDblFinder_score', 'scDblFinder_class'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells', 'percent_cells', 'robust', 'mean', 'var', 'residual_variances', 'highly_variable_rank', 'highly_variable_features'
    uns: 'layers_counts', 'log1p', 'hvg'
    layers: 'counts', 'soupX_counts'

注意，與scanpy不同，omicverse計(jì)算高可變基因后，將保存為var['highly_variable_features']，而在scanpy中，高可變將保存為var['highly_variable']

adata.var.head()

adata.write_h5ad('s4d8_preprocess.h5ad')

5. 思考

我們?yōu)槭裁匆x擇高可變基因/特征？
不同的高可變基因/特征最后只獲得了656個(gè)基因?yàn)榻患@是為什么？
選擇2000個(gè)高可變基因還是選擇3000個(gè)高可變基因，你認(rèn)為區(qū)別是什么？

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

單細(xì)胞測序最好的教程（三）：特征基因選擇

單細(xì)胞測序最好的教程（三）：特征基因選擇

1. 背景

2. 基于基因離散度

2.1 指定高可變基因數(shù)

2.2 指定基因離散度與平均值閾值

3. 基于基因歸一化方差

4. 基于基因皮爾森近似殘差

5. 思考

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

單細(xì)胞測序最好的教程（三）：特征基因選擇

1. 背景

2. 基于基因離散度

2.1 指定高可變基因數(shù)

2.2 指定基因離散度與平均值閾值

3. 基于基因歸一化方差

4. 基于基因皮爾森近似殘差

5. 思考

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av