单细胞kallisto传统流程

发表于2021-12-04更新于2022-10-21

kallisto/bustool/scanpy流程

通常，10X的上游都是cellranger，但是这玩意消耗大，速度慢，而且定制化程度很高，cellranger虽然是集成的STAR的流程，但是下游结果文件固定，参数也很多。个人用起来太重量级了。因此，有了某些大佬，决定针对scRNA的实际需求，改进了kallisto工具，并结合高端的bustool，因此就产生了这个上游kallisto+bustools，下游scanpy+其他高级可视化的workflow。

Node1.data

10X官网下载scRNA的测序数据即可，只需要输入一些基本信息，即可随意下载，根据自己的电脑配置下载，我这里下载的是1K PBMC的，大了放不下，最近动漫下多了。下载数据没什么特别的，但是需要注意，文库版本是哪个版本的，10X的V2文库和V3文库差别很大，而且这个参数后面需要。

Node2.Ref

kallisto是把fastq与转录组参考相对比，因此需要的并不是传统的转录组下载基因组的参考，而是需要去ensembl下载转录组的文件，注意，转录组，转录组，转录组，几个文件注意版本对应分别为：

1
2
3
Homo_sapiens.GRCh38.cdna.all.fa.gz #这个是cDNA的，也就是转录本参考
Homo_sapiens.GRCh38.104.gtf.gz #注意要去转录本的文件夹里下载转录本的gtf
Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz #这个是基因组的

借助一个叫kb-python的warpper，可以很容易的使用kallisto和bustool的联合，这里贴出官方的参数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
$ kb ref
usage: kb ref [-h] [--tmp TMP] [--keep-tmp] [--verbose] -i INDEX -g T2G -f1
              FASTA [-f2 FASTA] [-c1 T2C] [-c2 T2C] [-n N]
              [-d {human,mouse,linnarsson}] [-k K]
              [--workflow {standard,lamanno,nucleus,kite}] [--lamanno]
              [--overwrite]
              fasta gtf [feature]

Build a kallisto index and transcript-to-gene mapping

positional arguments:
  fasta                 Genomic FASTA file(s), comma-delimited
  gtf                   Reference GTF file(s), comma-delimited
  feature               [`kite` workflow only] Path to TSV containing barcodes
                        and feature names.

optional arguments:
  -h, --help            Show this help message and exit
  --tmp TMP             Override default temporary directory
  --keep-tmp            Do not delete the tmp directory
  --verbose             Print debugging information
  -n N                  Number of files to split the index into. If this
                        option is specified, the FASTA that is normally used
                        to create an index is split into `N` approximately-
                        equal parts. Each of these FASTAs are indexed
                        separately.
  -d {human,mouse,linnarsson}
                        Download a pre-built kallisto index (along with all
                        necessary files) instead of building it locally
  -k K                  Use this option to override the k-mer length of the
                        index. Usually, the k-mer length automatically
                        calculated by `kb` provides the best results.
  --workflow {standard,lamanno,nucleus,kite}
                        Type of workflow to prepare files for. Use `lamanno`
                        for RNA velocity based on La Manno et al. 2018 logic.
                        Use `nucleus` for RNA velocity on single-nucleus RNA-
                        seq reads. Use `kite` for feature barcoding. (default:
                        standard)
  --lamanno             Deprecated. Use `--workflow lamanno` instead.
  --overwrite           Overwrite existing kallisto index

required arguments:
  -i INDEX              Path to the kallisto index to be constructed. If `-n`
                        is also specified, this is the prefix for the n
                        indices to construct.
  -g T2G                Path to transcript-to-gene mapping to be generated
  -f1 FASTA             [Optional with -d] Path to the cDNA FASTA (lamanno,
                        nucleus) or mismatch FASTA (kite) to be generated

required arguments for `lamanno` and `nucleus` workflows:
  -f2 FASTA             Path to the intron FASTA to be generated
  -c1 T2C               Path to generate cDNA transcripts-to-capture
  -c2 T2C               Path to generate intron transcripts-to-capture

特别注意，如果在一开始加了-d参数，那么会从一个在线存储里面下载预先已经建立好的index，然而不出所料是下载不了的。根据后面的教程，如果要走standard的单细胞流程，需要2个文件，一个idx文件，就是索引，和一个t2g文件，transcripts-to-capture，就是一个基因注释文件，大概长这样，理论上这东西直接自己提取也不是不行……

1
2
3
4
5
6
7
8
9
10
11
12
13
14
NST00000456328.2    ENSG00000223972.5    DDX11L1
ENST00000450305.2    ENSG00000223972.5    DDX11L1
ENST00000488147.1    ENSG00000227232.5    WASH7P
ENST00000619216.1    ENSG00000278267.1    MIR6859-1
ENST00000473358.1    ENSG00000243485.5    MIR1302-2HG
ENST00000469289.1    ENSG00000243485.5    MIR1302-2HG
ENST00000607096.1    ENSG00000284332.1    MIR1302-2
ENST00000417324.1    ENSG00000237613.2    FAM138A
ENST00000461467.1    ENSG00000237613.2    FAM138A
ENST00000606857.1    ENSG00000268020.3    OR4G4P
ENST00000642116.1    ENSG00000240361.2    OR4G11P
ENST00000492842.2    ENSG00000240361.2    OR4G11P
ENST00000641515.2    ENSG00000186092.6    OR4F5
ENST00000335137.4    ENSG00000186092.6    OR4F5

还要注意，里面有个workflow参数，是用来指定对应流程的，这个教程先不考虑其他，直接走最经典的单细胞转录组流程，那么这个参数可以不动，默认就是standard，但是，如果准备走最早的velocity流程，那么需要使用lamanno，如果需要走单细胞核测序的RNA velocity（就是scVelo），需要使用nucleus，而要走10X的Feature Barcoding（测表面蛋白）的话，则需要使用kite。

这里我们走经典转录组，而且不需要太复杂，因此脚本可以为stdindex.sh

1
2
3
4
5
6
7
8
9
#!/usr/bin

transindexpath=../refdata/standard/trans104.idx
g2fpath=../refdata/standard/g2f.txt
transfastapath="../refdata/standard/cDNA.fa
genomepath=../refdata/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gtfpath=../refdata/Homo_sapiens.GRCh38.104.gtf.gz

kb ref -i $transindexpath -g $g2fpath -f1 $transfastapath $fastapath $gtfpath

运行大概持续了15min，速度还是相当快的，而且结果和直接下载得到的cdna.all.fa和idex文件基本一致，当然因为版本原因大小有一点区别的还是。

这里需要注意，如果我们下载了cdna.all.fa文件，最后输出的时候如果同文件夹下存在的话，会跳过这一步的输出。最后idx大小大概是3G多，g2f文件10多M，cDNA文件400多M。

所以，这个就是经典scRNA-seq流程的准备文件了，下面就是准备进入分析流程了。

Node3.count

这一步仍然是利用kb-python来进行比对和计数，由于kb-python的封装，这一步之后就能直接出所有的下游需要的文件了。以及…其实kb-python本身也就只有2个外部接口罢了。

这里也给出kb count的参数，可以看到，这里多了很多就是了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
usage: kb count [-h] [--tmp TMP] [--keep-tmp] [--verbose] -i INDEX -g T2G -x
                TECHNOLOGY [-o OUT] [-w WHITELIST] [-t THREADS] [-m MEMORY]
                [--workflow {standard,lamanno,nucleus,kite,kite:10xFB}]
                [--mm | --tcc] [--filter [{bustools}]] [-c1 T2C] [-c2 T2C]
                [--overwrite] [--dry-run] [--lamanno | --nucleus]
                [--loom | --h5ad]
                fastqs [fastqs ...]

Generate count matrices from a set of single-cell FASTQ files. Run `kb --list`
to view single-cell technology information.

positional arguments:
  fastqs                FASTQ files

optional arguments:
  -h, --help            Show this help message and exit
  --tmp TMP             Override default temporary directory
  --keep-tmp            Do not delete the tmp directory
  --verbose             Print debugging information
  -o OUT                Path to output directory (default: current directory)
  -w WHITELIST          Path to file of whitelisted barcodes to correct to. If
                        not provided and bustools supports the technology, a
                        pre-packaged whitelist is used. If not, the bustools
                        whitelist command is used. (`kb --list` to view
                        whitelists)
  -t THREADS            Number of threads to use (default: 8)
  -m MEMORY             Maximum memory used (default: 4G)
  --workflow {standard,lamanno,nucleus,kite,kite:10xFB}
                        Type of workflow. Use `lamanno` for RNA velocity based
                        on La Manno et al. 2018 logic. Use `nucleus` for RNA
                        velocity on single-nucleus RNA-seq reads. Use `kite`
                        for feature barcoding. Use `kite:10xFB` for 10x
                        Genomics Feature Barcoding technology. (default:
                        standard)
  --mm                  Include reads that pseudoalign to multiple genes.
  --tcc                 Generate a TCC matrix instead of a gene count matrix.
  --filter [{bustools}]
                        Produce a filtered gene count matrix (default:
                        bustools)
  --overwrite           Overwrite existing output.bus file
  --dry-run             Dry run
  --lamanno             Deprecated. Use `--workflow lamanno` instead.
  --nucleus             Deprecated. Use `--workflow nucleus` instead.
  --loom                Generate loom file from count matrix
  --h5ad                Generate h5ad file from count matrix

required arguments:
  -i INDEX              Path to kallisto index/indices, comma-delimited
  -g T2G                Path to transcript-to-gene mapping
  -x TECHNOLOGY         Single-cell technology used (`kb --list` to view)

required arguments for `lamanno` and `nucleus` workflows:
  -c1 T2C               Path to cDNA transcripts-to-capture
  -c2 T2C               Path to intron transcripts-to-captured

同样，可以使用stdcount.sh来进行这一步:

1
2
3
4
5
6
#!/usr/bin

fastqpath=../refdata/standard/trans104.idx
g2fpath=../refdata/standard/g2f.txt

kb count --h5ad -i $fastqpath -g $g2fpath -x 10xv3 -o output --filter bustools -t 2 data/*_001.fastq.gz

这里其实挺奇怪的，如果用过cellranger就能知道，用cellranger进行每个lane比对的时候是需要I1，R1和R2三个文件才能进行比对的，但是其实kallisto或者单独的STAR solo以及其他一些工具在做这一步的时候反而不需要I1那个文件，只需要R1和R2就行。

先看看这三个文件是什么：

“The source material consists of 6 FASTQ files split into two sequencing lanes L001 and L002, each with three reads of R1 (barcodes), R2 (cDNA sequences), I1 (illumina lane info)”

这三个文件的来源是下机的原始文件经过mkfastq或者illumina的bcl2fastq出来的，那么，需要R1和R2就很容易理解，barcode和cDNA自然是最需要的。但是为什么可以不要lane info？而且，上面的stdcount.sh会一并将lane info文件进入比对流程，但是仍然能够正常比对，尝试去掉lane info文件后再比对，也能比对出来。

但是，如果加入I1文件的话，比对的结果只有600M左右，但不加上I1文件则与官网给出的结果差不多，有1.3G，首先正确性上来说，肯定是不加I1文件是正确的结果，可以先看看这个所谓的lane info是什么：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
head pbmc_1k_v3_S1_L001_I1_001.fastq

@A00228:279:HFWFVDMXX:1:1101:8486:1000 1:N:0:NCATTACT
NCATTACT
+
#FFFFFFF
@A00228:279:HFWFVDMXX:1:1101:10782:1000 1:N:0:NCATTACT
NCATTACT
+
#FFFFFFF
@A00228:279:HFWFVDMXX:1:1101:12626:1000 1:N:0:NCATTACT
NCATTACT
+
#FFFFFFF

tail pbmc_1k_v3_S1_L001_I1_001.fastq
@A00228:279:HFWFVDMXX:1:2488:9670:37059 1:N:0:GGCAATGG
GGCAATGG
+
FFFFFFFF
@A00228:279:HFWFVDMXX:1:2488:13467:37059 1:N:0:GGCAATGG
GGCAATGG
+
FFFFFFFF

首先能看到，基本上全部是一样的东西，而且一个特征是，这个里面是不带有任何序列信息的，唯一长得像的就是比如NCATTACT，GGCAATGG这种了，那么这个8碱基的东西是啥呢。

根据10X V3的建库的官方说明以及提供的Single_Index_Kit_T_Set_A.csv这个文件，这个8碱基的东西应该是所谓的sample_index，一般就是使用的上面那个文件里面提供的官方的sample index序列，当然，自己指定也行，这个东西的作用很简答，就是用来区分样本用的。

根据下游的分析，大致能够判断，将sample index一起放进去比对，对上游比对结果影响不是特别大，当然，是有影响的，比如使用bustools inspect查看比对结果时能发现，将index一起放进比对的结果里面mean reads per barcode少了不少，甚至已经低于V3的标准（2W reads per barcode），理论上结果应当以不加入sample index的结果为准。

而且，从最后聚类的结果也能明显发现，不加sample index的时候，聚类结果和scanpy的官网示例基本一致，都是聚类为0-7，但是加了sample index则变成了0-9，可以认为加了sample index进比对流程的话，会影响比对出来的gene的count数。

总之，最后count的结果，就是output文件夹里的东西，大概有这些：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
output
├── 10x_version3_whitelist.txt
├── counts_filtered
│   ├── adata.h5ad
│   ├── cells_x_genes.barcodes.txt
│   ├── cells_x_genes.genes.txt
│   └── cells_x_genes.mtx
├── counts_unfiltered
│   ├── adata.h5ad
│   ├── cells_x_genes.barcodes.txt
│   ├── cells_x_genes.genes.txt
│   └── cells_x_genes.mtx
├── filter_barcodes.txt
├── inspect.json
├── kb_info.json
├── matrix.ec
├── output.bus
├── output.filtered.bus
├── output.unfiltered.bus
├── run_info.json
└── transcripts.txt

2 directories, 18 files

其中bus文件是比对后通过bustools压缩的文件，而标准的10X流程则只需要counts_unfiltered/filtered中的h5ad文件或者经典三文件即可。

Node4.DownStream-Data Preprocess

那么，继续往下，以不加sample index的结果为准，下游就能进入scanpy或者scater的标准QC过程了，这里以scanpy为主：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import numpy as np
import pandas as pd
import scanpy as sc
import anndata
from sklearn.decomposition import TruncatedSVD
import matplotlib
import matplotlib.pyplot as plt

sc.settings.verbosity = 3
sc.settings.set_figure_params(dpi=80)


adata = anndata.read_h5ad("output/counts_unfiltered/adata.h5ad")
adata.var["gene_id"] = adata.var.index.values

t2g = pd.read_table("../refdata/standard/g2f.txt", header=None, names=["tid", "gene_id", "gene_name","version_name","chr","start","end","strand"], sep="\t")
t2g = t2g.dropna(axis=0)
t2g.index = t2g.gene_id
t2g = t2g.loc[~t2g.index.duplicated(keep='first')]

adata.var["gene_name"] = adata.var.gene_id.map(t2g["gene_name"])
adata.var = adata.var.fillna("unknowngene")
adata.var.index = adata.var["gene_name"].tolist()

adata.var_names_make_unique()

这里有2个超级坑，其一，是如果使用最新版本的基因组和gtf文件，可能出现有个ensembl-id没有对应gene_name的情况，这很正常，但是进入anndata会导致后续无法分析，而且scannpy又没有提供去除Na或者空值之类的接口，因此，简单起见，如果使用了版本较高的基因组和gtf文件，务必检查前面ref的时候的那个g2f.txt里面又没有缺失或者空值，如果有，务必对anndata里面的var进行处理，我这里用的是一个比较粗暴的办法，直接把缺失设置为unknowngene，后面能够正常进行。

1
adata.var = adata.var.fillna("unknowngene")

其二，在设置var的index的时候，必须如下设置：

1
adata.var.index = adata.var["gene_name"].tolist()

官网教程上要么是给的干净数据，要么是之前版本的接口，没有对index对象做tolist处理，这样会导致anndata的var.index设置错误，进一步导致后续部分步骤和所有的结果无法保存，而且一旦session关闭，这次处理过后的h5ad文件将无法再被读入python，等同于要重新跑。

继续，读入h5ad之后，就构成了一个叫anndata的对象，这个对象跟R里面的# SummarizedExperiment/Singlecellexperiment对象特别像。根据官方文档，anndata至少有3个部分

anndata.var：主要是基因级别信息，比如上面的gene_name这些。

anndata.obs：主要是cell级别的信息，比如barcode，group这些

anndata.X：核心数据，用于存储稀疏矩阵

比如在最后分析结束后，adata里面有这些东西：

1
2
3
4
5
6
7
8
9
adata

AnnData object with n_obs × n_vars = 1175 × 25834
    obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain', 'leiden'
    var: 'gene_name', 'gene_id', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'
    uns: 'log1p', 'hvg', 'pca', 'neighbors', 'louvain', 'paga', 'louvain_sizes', 'umap', 'leiden', 'leiden_colors', 'louvain_colors'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    obsp: 'distances', 'connectivities'

可以看到，还有一些没说的属性，比如uns，这里存放了一些不能被结构话的数据，比如每个group绘图的时候用什么颜色映射，paga的路径之类的，里面大概长这样:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
adata.uns

OverloadedDict, wrapping:
    OrderedDict([('log1p', {'base': None}), ('hvg', {'flavor': 'seurat'}), ('pca', {'params': {'zero_center': True, 'use_highly_variable': True}, 'variance': array([66.33609  , 48.00583  , 43.741863 , 24.337816 , 11.7938385,
        9.654164 ,  8.945425 ,  8.060145 ,  7.9558783,  7.075499 ,
        6.960134 ,  6.713339 ,  6.5481   ,  6.3729434,  6.224372 ,
        6.066858 ,  5.907557 ,  5.595799 ,  5.498731 ,  5.339343 ,
        5.264697 ,  5.085398 ,  5.0702558,  4.9665728,  4.9449844,
        4.8837957,  4.7512255,  4.7258067,  4.6373615,  4.589517 ,
        4.484859 ,  4.4121437,  4.3468637,  4.294298 ,  4.2609587,
        4.2075477,  4.201293 ,  4.138263 ,  4.093854 ,  4.079861 ,
        4.035178 ,  4.023861 ,  3.9789526,  3.9368029,  3.9152558,
        3.8964837,  3.8508425,  3.766862 ,  3.7391431,  3.7227194],
      dtype=float32), 'variance_ratio': array([0.04367469, 0.03160632, 0.02879899, 0.01602365, 0.00776489,
       0.00635616, 0.00588953, 0.00530668, 0.00523803, 0.0046584 ,
       0.00458245, 0.00441996, 0.00431117, 0.00419585, 0.00409803,
       0.00399433, 0.00388945, 0.00368419, 0.00362028, 0.00351534,
       0.0034662 , 0.00334815, 0.00333818, 0.00326992, 0.0032557 ,
       0.00321542, 0.00312814, 0.0031114 , 0.00305317, 0.00302167,
       0.00295276, 0.00290489, 0.00286191, 0.0028273 , 0.00280535,
       0.00277019, 0.00276607, 0.00272457, 0.00269533, 0.00268612,
       0.0026567 , 0.00264925, 0.00261968, 0.00259193, 0.00257775,
       0.00256539, 0.00253534, 0.00248005, 0.0024618 , 0.00245098],
      dtype=float32)}), ('neighbors', {'connectivities_key': 'connectivities', 'distances_key': 'distances', 'params': {'n_neighbors': 10, 'method': 'umap', 'random_state': 0, 'metric': 'euclidean', 'n_pcs': 10}}), ('louvain', {'params': {'resolution': 0.7, 'random_state': 42}}), ('paga', {'connectivities': <6x6 sparse matrix of type '<class 'numpy.float64'>'
    with 14 stored elements in Compressed Sparse Row format>, 'connectivities_tree': <6x6 sparse matrix of type '<class 'numpy.float64'>'
    with 5 stored elements in Compressed Sparse Row format>, 'groups': 'louvain', 'pos': array([[  9.89749122, -11.9940879 ],
       [ 12.33582657,  -6.49087334],
       [ 16.42411062,  -6.62426792],
       [ 12.64145549,  -8.17883229],
       [  8.90823793, -14.36047967],
       [  7.95605852, -13.65646286]])}), ('louvain_sizes', array([407, 297, 188, 108, 102,  73])), ('umap', {'params': {'a': 0.5830300205483709, 'b': 1.334166992455648}}), ('leiden', {'params': {'resolution': 1, 'random_state': 0, 'n_iterations': -1}}), ('leiden_colors', ['#1f77b4', '#ff7f0e', '#279e68', '#d62728', '#aa40fc', '#8c564b', '#e377c2', '#b5bd61', '#17becf', '#aec7e8', '#ffbb78', '#98df8a', '#ff9896', '#c5b0d5', '#c49c94']), ('louvain_colors', ['#1f77b4', '#ff7f0e', '#279e68', '#d62728', '#aa40fc', '#8c564b'])])
With overloaded keys:
    ['neighbors'].

而adata.var则是这个样子:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
adata.var


     gene_name     gene_id     n_cells     highly_variable     means     dispersions     dispersions_norm     mean     std
TNFRSF4     TNFRSF4     ENSG00000186827.11     59     True     0.106400     1.206044     2.019484     0.054150     0.252805
TNFRSF18     TNFRSF18     ENSG00000186891.14     39     True     0.068430     1.109991     1.729919     0.035445     0.202188
ATAD3B     ATAD3B     ENSG00000160072.20     139     False     0.133877     0.300382     -0.710786     0.090129     0.260384
THAP3     THAP3     ENSG00000041988.15     75     False     0.078003     0.520083     -0.048459     0.049195     0.203547
unknowngene     unknowngene     ENSG00000260179.1     11     False     0.011523     0.469297     -0.201563     0.007132     0.078769
...     ...     ...     ...     ...     ...     ...     ...     ...     ...
unknowngene-20032     unknowngene     ENSG00000276256.1     7     False     0.009656     0.626575     0.272579     0.005590     0.074509
unknowngene-20035     unknowngene     ENSG00000273748.1     8     False     0.006963     0.207106     -0.991981     0.004621     0.059163
unknowngene-20038     unknowngene     ENSG00000278817.1     15     True     0.019199     1.074168     1.621924     0.010459     0.103165
unknowngene-20043     unknowngene     ENSG00000278384.1     61     False     0.067741     0.598364     0.187531     0.041752     0.191105
unknowngene-20052     unknowngene     ENSG00000271254.6     8     True     0.011055     1.205207     2.016962     0.005692     0.078602

25834 rows × 9 columns

adata.obs则是：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
adata.obs


n_genes     percent_mito     n_counts     louvain     leiden
barcode                     
AAACCCAAGGAGAGTA     3137     0.112912     9131.0     1     4
AAACGCTTCAGCCCAG     2497     0.081725     6424.0     2     3
AAAGAACAGACGACTG     2059     0.062348     4940.0     4     10
AAAGAACCAATGGCAG     1540     0.068987     3218.0     4     10
AAAGAACGTCTGCAAT     2484     0.066250     7366.0     0     8
...     ...     ...     ...     ...     ...
TTTCCTCTCTCTTGCG     3834     0.116246     12723.0     3     11
TTTGATCTCTTTGGAG     2567     0.086051     8065.0     0     8
TTTGGTTAGTAACCTC     2012     0.083509     5233.0     0     0
TTTGGTTGTAGAATAC     3365     0.121682     10585.0     1     1
TTTGTTGCAATTAGGA     2001     0.085224     5069.0     2     3

1175 rows × 5 columns

adata.X则是：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
adata.X


array([[-0.214197  , -0.17530552, -0.34613955, ..., -0.10138597,
        -0.21847485, -0.0724164 ],
       [-0.214197  , -0.17530552, -0.34613955, ..., -0.10138597,
        -0.21847485, -0.0724164 ],
       [-0.214197  , -0.17530552, -0.34613955, ..., -0.10138597,
        -0.21847485, -0.0724164 ],
       ...,
       [-0.214197  , -0.17530552, -0.34613955, ..., -0.10138597,
        -0.21847485, -0.0724164 ],
       [-0.214197  , -0.17530552,  2.2082584 , ..., -0.10138597,
        -0.21847485, -0.0724164 ],
       [-0.214197  , -0.17530552, -0.34613955, ..., -0.10138597,
        -0.21847485, -0.0724164 ]], dtype=float32)

从形式也能看到，var和obs是用DataFrame对象存储的，因此，各种pandas的方法都能直接在上面应用，进行一些统计，而X则是np的array，对其进行操作时则需要用numpy的操作方式。而其他非结构化的数据更多的是代码运行中需要的，对研究来说没有特殊意义。此外，知道了这个之后，我们其实就能够手动向var和obs里面添加信息，比如添加gene的type，染色体的正负链，亦或者是样本的临床信息，病理信息，这样的组合型对象非常灵活，对后续分析帮助很大。

下面继续QC的过程，首先，由于这算是比对后的数据，因此需要先在UMI，barcode级别进行一下质控，本质类似于cellranger跑出来的那个html文件，当然，上面的bustools的部分也完成了一点点质控的部分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
fig, ax = plt.subplots(figsize=(7, 7))

x = np.asarray(adata.X.sum(axis=1))[:,0]
y = np.asarray(np.sum(adata.X>0, axis=1))[:,0]

ax.scatter(x, y, color="green", alpha=0.25)
ax.set_xlabel("UMI Counts")
ax.set_ylabel("Genes Detected")
ax.set_xscale('log')
ax.set_yscale('log')

ax.set_xlim(1)
ax.set_ylim(1)


plt.show()

这里会出来一个这个“Create a plot showing genes detected as a function of UMI counts.”就是一张横轴为UMI数量，纵轴为基因数量的散点图，如果这张图的趋势基本趋于缓和，说明基本覆盖度已经达到饱和了，在增加测序深度也不会有更多的基因覆盖，当然在1K的PBMC数据里，其实数据基本分布于对角线，缓和趋势不是特别明显，示例数据嘛不强求。

然后就是根据UMI的分布估计细胞数目，这个图应该是最为常见的：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
expected_num_cells =  1200 #@param {type:"integer"}
knee = np.sort((np.array(adata.X.sum(axis=1))).flatten())[::-1]

fig, ax = plt.subplots(figsize=(10, 7))

ax.loglog(knee, range(len(knee)), linewidth=5, color="g")
ax.axvline(x=knee[expected_num_cells], linewidth=3, color="k")
ax.axhline(y=expected_num_cells, linewidth=3, color="k")

ax.set_xlabel("UMI Counts")
ax.set_ylabel("Set of Barcodes")

plt.grid(True, which="both")
plt.show()

其实这个图可以和上面bustools inspect的结果里面的预估唯一barcode数进行对照来看，如果偏移得很离谱说明数据有问题，可能是index有问题，或者barcode文件有问题，或者就是测序质量特别差导致UMI分布不均匀。不过一般问题不大。

然后需要看看UMI对应最多的基因是什么（等同于表达量最高的基因）

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

这里调用的是scanpy的计算高表达基因的功能，本质上就是计算基因在UMI上的分配情况，可以直接理解为基因的表达量。这里官方给了说明，通常，最容易出现高表达的就是MALAT1，MT基因，RP基因，actin基因，血红蛋白基因，假基因或者ERCC，这是正常的现象。但是当这些基因类型表达高得离谱得时候，比如top30全是某一类基因的时候，就分别与几种异常情况对应了：其中MT基因特别多，说明太多代表细胞存在死亡或者细胞质泄漏的问题，如果ERCC较多说明建库时ERCC加多了技术误差就没法矫正，假基因太多说明比对存在问题或者基因组版本用错了。对于这些基因是否需要先彻底去除后进行分析，目前没有绝对的定论，可以通常是针对细胞进行去除，但是对这些基因本身不会直接去掉。

在上面对UMI和barcode的质控没有问题后，进入常见单细胞最熟悉的环节，开始对细胞和基因进行操作，准备进行下游分析了。

1
2
3
4
5
# Removes cells with less than 1070 umi counts
adata = adata[np.asarray(adata.X.sum(axis=1)).reshape(-1) > 1070]

# Removes genes with 0 umi counts（也就是有基因，但是没有转录本）
adata = adata[:, np.asarray(adata.X.sum(axis=0)).reshape(-1) > 0]

原文这里是说是根据阈值进行过滤，它说是过滤低于1170umi count的细胞，代码是中关于这个axis=i，我看到的最经典的解释是：下标的变化方向；axis=0就是下标沿着第0轴递归操作，axis=1就是下标沿着第1轴递归操作。

第二个操作的目的很容易懂，就是将UMI为0的基因去除，等同于去掉没有测到的基因。但是第一个操作很有意思，它去掉了低于UMI=1070的细胞，但是这个阈值是从哪里来的呢。
顺便注意，经过下面两个操作，现在adata本身已经被改变了。因此，可以先将原本的矩阵进行一次保存到一个anndata内置的raw属性里面，这样就能在改动矩阵数据的同时保留原本的表达数据。

接下来就开始正式进入“传统”的表达矩阵类似级别的分析了

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

这里就跟seurat非常像了，先筛选细胞里表达不到200个基因的，再筛选在小于3个细胞中表达的基因说白了就是一个是筛低质量细胞，一个筛低质量基因。

1
2
3
4
mito_genes = adata.var_names.str.startswith('MT-')
adata.obs['percent_mito'] = np.sum(adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1
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)

查看质量绘图，对应seurat里的对应的图。
n_gene代表细胞内的基因数，可以看到基本在2000-4000（其实貌似也就一般……）
nCount可以理解为测序深度（文库大小），后续跟文库normalization有关。
percent.mt就是线粒体基因含量。

再根据上面的数据进行一次过滤

1
2
mask = np.logical_or((adata.obs.n_genes < 6500).values, (adata.obs.n_genes > 200).values, (adata.obs.percent_mito < 0.2).values)
adata = adata[mask, :]

然后进行文库归一化，之前我们看到了，测序深度这个问题，并不是在所有细胞中都是一样的（也就是reads不是一样多，基因覆盖度也不一样）。这种情况导致根据UMI弄到的基因并不可比，所以需要矫正，对应的就是Normalize counts，把文库都归一化道10000（也就是seurat里的sizefactor）

1
2
# normalize counts in each cell to be equal
sc.pp.normalize_total(adata, target_sum=10**4)

然后取表达值的log压缩表达量的范围

1
sc.pp.log1p(adata)

再来选择高变基因，跟Seurat里一样

1
2
3
4
5
sc.pp.highly_variable_genes(adata, min_mean=0.01, max_mean=8, min_disp=1, n_top_genes=2000, flavor="seurat", n_bins=20)
sc.pl.highly_variable_genes(adata)

#此外，虽说提供了“消除”指定“因素”影响的功能，类似于seurat里消除细胞周期的影响。
# sc.pp.regress_out(adata, ['n_counts', 'percent_mito'])

然后进行scale:

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

其实就是将数据变为0位均值，10为方差的正态分布，才能进行pca。所以为啥要选方差为10…..

Node5.DownStream-Clustering

一切准备就绪，上游QC，数据变换都执行完了，准备开始进入聚类分群了

首先是PCA

1
2
3
4
5
# We perform PCA on just the highly variable genes
sc.tl.pca(adata, svd_solver='arpack', use_highly_variable=True)
sc.pl.pca(adata, color='CST3')
sc.pl.pca_variance_ratio(adata, log=True)
sc.pp.neighbors(adata, n_neighbors=20, n_pcs=10)

这里的参数也就svd的方法可以改动，其他基本没什么可以改的，注意，如果在之前实际上能看到和示例并不完全相同，应该就是前面flavor的问题，但其实这里效果也挺好的，能看到分得非常开就是了。

此外，scanpy官方教程都说了

often, a rough estimate of the number of PCs does fine.

然后就可以执行Cluster了，这里有两个路径，其一是：

louvain+paga+umap；

1
2
3
4
sc.tl.louvain(adata)
sc.tl.paga(adata)
sc.pl.paga(adata, plot=False)
sc.tl.umap(adata, init_pos='paga')

第二则是：

ledian+umap

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

理论上，louvain和ledian不会有太大的区别，一般二者聚类的结果都是基本一致的。

然后就是挑选marker基因进行注释了。这一步这里不讲，这一次主要是上游区别与cellranger的一个流程汇总，注释放在可视化的部分里面专门讲，大概是下下次，下一次就就做一下cellranger的流程，二者用同一份数据对比看看有什么区别。