R在讀取和處理數(shù)據(jù)的過程中會將所有的變量和占用都儲存在RAM當(dāng)中，這樣一來，對于海量的單細(xì)胞RNA-seq數(shù)據(jù)（尤其是超過250k的細(xì)胞量），即使在服務(wù)器當(dāng)中運行，Seurat、metacell、monocle這一類的R包的使用還是會產(chǎn)生內(nèi)存不足的問題。

但是最近我發(fā)現(xiàn)了一個基于python的單細(xì)胞基因表達(dá)分析包scanpy，能夠很好地在我這個僅4G內(nèi)存的小破機(jī)上分析378k的細(xì)胞，并且功能豐富程度不亞于Seurat。它包含了數(shù)據(jù)預(yù)處理、可視化、聚類、偽時間分析和軌跡推斷、差異表達(dá)分析。根據(jù)官網(wǎng)描述，scanpy能夠有效處理超過1,000,000個細(xì)胞的數(shù)據(jù)集。

Scanpy – Single-Cell Analysis in Python：https://scanpy.readthedocs.io/en/latest/

安裝與數(shù)據(jù)下載

安裝好python3之后，終端運行：

pip install scanpy

若安裝過程出現(xiàn)問題，請參考:
https://scanpy.readthedocs.io/en/latest/installation.html

骨髓單細(xì)胞轉(zhuǎn)錄組測序數(shù)據(jù)下載地址：
https://preview.data.humancellatlas.org

Step0, 讀取數(shù)據(jù)

運行python

import numpy as np
import pandas as pd
import scanpy as sc
# 可以直接讀取10Xgenomics的.h5格式數(shù)據(jù)
adata = sc.read_10x_h5("/Users/shinianyike/Desktop/ica_bone_marrow_h5.h5"
               , genome=None, gex_only=True)
adata.var_names_make_unique()

查看數(shù)據(jù)：

adata

    AnnData object with n_obs × n_vars = 378000 × 33694 
        var: 'gene_ids'

共378000個細(xì)胞，33694個基因。

Step1, 數(shù)據(jù)預(yù)處理

這一步目的將數(shù)據(jù)進(jìn)行篩選和過濾，一些測序質(zhì)量差的數(shù)據(jù)會讓后續(xù)的分析產(chǎn)生噪音和干擾，因此我們需要將它們?nèi)コ?/p>

展示在所有的細(xì)胞當(dāng)中表達(dá)占比最高的20個基因。

sc.pl.highest_expr_genes(adata, n_top=20)

表達(dá)水平前20的基因.png

基礎(chǔ)過濾：
去除表達(dá)基因200以下的細(xì)胞；
去除在3個細(xì)胞以下表達(dá)的基因。

sc.pp.filter_cells(adata, min_genes=200)   # 去除表達(dá)基因200以下的細(xì)胞
sc.pp.filter_genes(adata, min_cells=3)     # 去除在3個細(xì)胞以下表達(dá)的基因

adata

    AnnData object with n_obs × n_vars = 335618 × 24888 
        obs: 'n_genes'
        var: 'gene_ids', 'n_cells'

共335618個細(xì)胞，24888個基因。(可以發(fā)現(xiàn)細(xì)胞跟基因數(shù)量都減少了)

質(zhì)量控制：

在測序后，有很大比例是低質(zhì)量的細(xì)胞，可能是因為細(xì)胞破碎造成的胞質(zhì)RNA流失，由于線粒體比單個的轉(zhuǎn)錄分子要大得多，不容易在破碎的細(xì)胞膜中泄漏出去，這樣一來就導(dǎo)致測序結(jié)果顯示線粒體基因的比例在細(xì)胞內(nèi)占比異常高。這一步質(zhì)量控制的目的就是將這些低質(zhì)量的細(xì)胞去除掉。

計算線粒體基因占所有基因的比例:

mito_genes = adata.var_names.str.startswith('MT-')
# for each cell compute fraction of counts in mito genes vs. all genes
# the `.A1` is only necessary as X is sparse (to transform to a dense array after summing)
adata.obs['percent_mito'] = np.sum(
    adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1
# add the total counts per cell as observations-annotation to adata
adata.obs['n_counts'] = adata.X.sum(axis=1).A1

作小提琴圖，查看線粒體基因占比分布:

sc.pl.violin(adata, ['n_genes', 'n_counts', 'percent_mito'],
             jitter=0.4, multi_panel=True)

細(xì)胞表達(dá)的基因種數(shù)、基因數(shù)、線粒體基因占比

由于數(shù)據(jù)點實在太多，小提琴已被數(shù)據(jù)點覆蓋，無法顯示出來。

這里還可以做一個散點圖，也可以直觀地顯示出一些異常分布的數(shù)據(jù)點。

sc.pl.scatter(adata, x='n_counts', y='percent_mito')
sc.pl.scatter(adata, x='n_counts', y='n_genes')

adata

    AnnData object with n_obs × n_vars = 335618 × 24888 
        obs: 'n_genes', 'percent_mito', 'n_counts'
        var: 'gene_ids', 'n_cells'

335618個細(xì)胞，24888個基因；
下面這一步進(jìn)行真正的過濾，根據(jù)上面做的圖，去除基因數(shù)目過多（大于等于4000）和線粒體基因占比過多（大于等于0.3）的細(xì)胞：

adata = adata[adata.obs['n_genes'] < 4000, :]
adata = adata[adata.obs['percent_mito'] < 0.3, :]

過濾后查看剩下多少細(xì)胞和基因。

adata

    View of AnnData object with n_obs × n_vars = 328435 × 24888 
        obs: 'n_genes', 'percent_mito', 'n_counts'
        var: 'gene_ids', 'n_cells'

328435個細(xì)胞，24888個基因。

數(shù)據(jù)標(biāo)準(zhǔn)化

在繪圖之前，還要進(jìn)行一步數(shù)據(jù)標(biāo)準(zhǔn)化，將表達(dá)量用對數(shù)計算一遍，這樣有利于繪圖和展示。

sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)
sc.pp.log1p(adata)

adata.raw = adata # 儲存標(biāo)準(zhǔn)化后的AnnaData Object

識別差異表達(dá)基因

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

sc.pl.highly_variable_genes(adata)

將保守的基因去除，留下差異表達(dá)的基因用于后續(xù)分析。

adata = adata[:, adata.var['highly_variable']]

adata

    View of AnnData object with n_obs × n_vars = 328435 × 1372 
        obs: 'n_genes', 'percent_mito', 'n_counts'
        var: 'gene_ids', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'

328435個細(xì)胞，1372個基因。

回歸每個細(xì)胞總計數(shù)和線粒體基因表達(dá)百分比的影響。將數(shù)據(jù)放縮到方差為1。單細(xì)胞數(shù)據(jù)集可能包含“不感興趣”的變異來源。這不僅包括技術(shù)噪音，還包括批次效應(yīng)，甚至包括生物變異來源（細(xì)胞周期階段）。正如(Buettner, et al NBT，2015)中所建議的那樣，從分析中回歸這些信號可以改善下游維數(shù)減少和聚類。
這一步高內(nèi)存需求預(yù)警，建議清理電腦緩存，關(guān)閉后臺不使用的應(yīng)用。

sc.pp.regress_out(adata, ['n_counts', 'percent_mito'])

    /Users/shinianyike/anaconda3/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
      from pandas.core import datetools
Scale each gene to unit variance. Clip values exceeding standard deviation 10.

sc.pp.scale(adata, max_value=10)

Step2, 主成分分析

主成分分析是一種將數(shù)據(jù)降維的分析方法，是考察多個變量間相關(guān)性一種多元統(tǒng)計方法，研究如何通過少數(shù)幾個主成分來揭示多個變量間的內(nèi)部結(jié)構(gòu)，即從原始變量中導(dǎo)出少數(shù)幾個主成分，使它們盡可能多地保留原始變量的信息，且彼此間互不相關(guān).通常數(shù)學(xué)上的處理就是將原來P個指標(biāo)作線性組合，作為新的綜合指標(biāo)。

sc.tl.pca(adata, svd_solver='arpack') # PCA分析
sc.pl.pca(adata, color='CST3') #繪圖

作碎石圖，觀測主成分的質(zhì)量。這個圖用于選擇后續(xù)應(yīng)該使用多少個PC，用于計算細(xì)胞間的相鄰距離。

sc.pl.pca_variance_ratio(adata, log=True)

在這里先將數(shù)據(jù)保存，便于后續(xù)操作的更改。

adata.write("pca_results.h5ad")

Step3， 聚類分析

計算細(xì)胞間的距離：
這里的參數(shù)就先按照默認(rèn)值設(shè)定：

sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

參數(shù)說明：
n_neighbors指的是每個點的鄰近點的數(shù)量，據(jù)評論區(qū)@小光amateur 所說neighbors的個數(shù)越多，聚類數(shù)會越少。
pc的數(shù)量依賴于上面所做的碎石圖，一般是選在拐點處的的主成分，只需要一個粗略值，不同的pc數(shù)量所產(chǎn)生的聚類形狀也不同。我后來更改為PC=16，效果比下圖要好一些。

將距離嵌入圖中：

sc.tl.umap(adata)
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])

sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=False)

聚類

sc.tl.louvain(adata)

sc.pl.umap(adata, color=['louvain'])

output_51_0.png

這里得到了29類細(xì)胞，每個顏色代表一種。
將數(shù)據(jù)保存。

adata.write("umap.h5ad")

尋找marker基因

marker基因通常是細(xì)胞表面抗原，用于定義出該細(xì)胞的細(xì)胞類型。
為了定義每個簇屬于什么細(xì)胞，根據(jù)基因的差異表達(dá)水平，將每個簇排名前25的基因?qū)С觥?/p>

sc.tl.rank_genes_groups(adata, 'louvain', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

下一步的工作是找出每一個簇的marker基因?qū)?yīng)的細(xì)胞類型，這主要依靠一些數(shù)據(jù)庫或生物學(xué)的相關(guān)背景知識。

文章已發(fā)布到微信公眾號：百味科研芝士，歡迎關(guān)注。