scater包简介

scater是一个优秀的单细胞转录组数据分析工具包，它可以对单细胞数据进行常规的质量控制，数据的标准化与归一化，以及数据的降维与可视化分析。它主要基于SingleCellExperiment类（来自SingleCellExperiment包）来进行操作处理，因此可以与其他许多Bioconductor包（如scran，batchelor和iSEE等）相互操作。

scater包主要含有以下特性：

  • Use of the SingleCellExperiment class as a data container for interoperability with a wide range of other Bioconductor packages;
  • Functions to import kallisto and Salmon results;
  • Simple calculation of many quality control metrics from the expression data;
  • Many tools for visualising scRNA-seq data, especially diagnostic plots for quality control;
  • Subsetting and many other methods for filtering out problematic cells and features;
  • Methods for identifying important experimental variables and normalising data ahead of downstream statistical analysis and modeling.

scater包的工作流程为：
使用scater包进行单细胞测序分析（一）：数据导入与SingleCellExperiment对象构建 - 图1

构建SingleCellExperiment对象

使用SingleCellExperiment函数导入单细胞转录组的基因表达矩阵构建一个SingleCellExperiment对象，该表达矩阵是一个行为基因，列为细胞的大型数据框。

SingleCellExperiment对象内容
使用scater包进行单细胞测序分析（一）：数据导入与SingleCellExperiment对象构建 - 图2

SingleCellExperiment对象常见操作
使用scater包进行单细胞测序分析（一）：数据导入与SingleCellExperiment对象构建 - 图3

  1. # 导入scater包
  2. library(scater)
  3. # 加载示例数据
  4. data("sc_example_counts")
  5. data("sc_example_cell_info")
  6. # 查看基因表达矩阵
  7. head(sc_example_counts)
  8.           Cell_001 Cell_002 Cell_003 Cell_004 Cell_005 Cell_006 Cell_007
  9. Gene_0001        0      123        2        0        0        0        0
  10. Gene_0002      575       65        3     1561     2311      160        2
  11. Gene_0003        0        0        0        0     1213        0        0
  12. Gene_0004        0        1        0        0        0       99      476
  13. Gene_0005        0        0       11        0        0        0        0
  14. Gene_0006        0        0        0        0        0        0      673
  15.           Cell_008 Cell_009 Cell_010 Cell_011 Cell_012 Cell_013 Cell_014
  16. Gene_0001       21        2        0     2624        1     1015        0
  17. Gene_0002        1        0        0        2        0     2710        0
  18. Gene_0003        1        0        0        2      178        0        0
  19. Gene_0004        0        1       66        0        1        0        1
  20. Gene_0005        0        1        0        0        2        2        0
  21. Gene_0006        0     3094        0        0      270        2        0
  22.           Cell_015 Cell_016 Cell_017 Cell_018 Cell_019 Cell_020 Cell_021
  23. Gene_0001        0        1       34        1        0        6        0
  24. Gene_0002        4        0      908      673      174      622     2085
  25. Gene_0003        0        0        0        0        1        0     3320
  26. Gene_0004        0      906      655     1020        1        0        0
  27. Gene_0005        0        0        0        2        0        0        3
  28. Gene_0006     1176        0        3        0        0        0        1
  30. # 查看样本信息
  31. head(sc_example_cell_info)
  32.              Cell Mutation_Status Cell_Cycle Treatment
  33. Cell_001 Cell_001        positive          S    treat1
  34. Cell_002 Cell_002        positive         G0    treat1
  35. Cell_003 Cell_003        negative         G1    treat1
  36. Cell_004 Cell_004        negative          S    treat1
  37. Cell_005 Cell_005        negative         G1    treat2
  38. Cell_006 Cell_006        negative         G0    treat1
  40. # 使用SingleCellExperiment函数构建SingleCellExperiment对象
  41. example_sce <- SingleCellExperiment(
  42.     assays = list(counts = sc_example_counts), 
  43.     colData = sc_example_cell_info
  44. )
  46. # 查看SingleCellExperiment对象
  47. example_sce
  48. class: SingleCellExperiment 
  49. dim: 2000 40 
  50. metadata(0):
  51. assays(1): counts
  52. rownames(2000): Gene_0001 Gene_0002 ... Gene_1999 Gene_2000
  53. rowData names(0):
  54. colnames(40): Cell_001 Cell_002 ... Cell_039 Cell_040
  55. colData names(4): Cell Mutation_Status Cell_Cycle Treatment
  56. reducedDimNames(0):
  57. spikeNames(0):
  59. View(example_sce)

使用scater包进行单细胞测序分析（一）：数据导入与SingleCellExperiment对象构建 - 图4

我们通常使用原始的count矩阵存储到SingleCellExperiment对象的“counts” Assay中，同时也可以使用counts函数提取SingleCellExperiment对象中的count表达矩阵。

  1. str(counts(example_sce))
  2.  int [1:2000, 1:40] 0 575 0 0 0 0 0 0 416 12 ...
  3.  - attr(*, "dimnames")=List of 2
  4.   ..$ : chr [1:2000] "Gene_0001" "Gene_0002" "Gene_0003" "Gene_0004" ...
  5.   ..$ : chr [1:40] "Cell_001" "Cell_002" "Cell_003" "Cell_004" ...
  7. head(counts(example_sce))

对于样本行和列的meta信息，我们也提供了一些常用函数来进行操作处理，如isSpike, sizeFactors, 和reducedDim等函数。

  1. # 添加一个新列的meta信息whee
  2. example_sce$whee <- sample(LETTERS, ncol(example_sce), replace=TRUE)
  4. # 使用colData函数查看列的meta信息
  5. colData(example_sce)
  6. DataFrame with 40 rows and 5 columns
  7.                 Cell Mutation_Status  Cell_Cycle   Treatment        whee
  8.          <character>     <character> <character> <character> <character>
  9. Cell_001    Cell_001        positive           S      treat1           N
  10. Cell_002    Cell_002        positive          G0      treat1           T
  11. Cell_003    Cell_003        negative          G1      treat1           Y
  12. Cell_004    Cell_004        negative           S      treat1           T
  13. Cell_005    Cell_005        negative          G1      treat2           C
  14. ...              ...             ...         ...         ...         ...
  15. Cell_036    Cell_036        negative          G0      treat1           Q
  16. Cell_037    Cell_037        negative          G0      treat1           X
  17. Cell_038    Cell_038        negative          G0      treat2           W
  18. Cell_039    Cell_039        negative          G1      treat1           B
  19. Cell_040    Cell_040        negative          G0      treat2           Z
  21. # 添加一个新行的meta信息
  22. rowData(example_sce)$stuff <- runif(nrow(example_sce))
  23. # 使用rowData函数查看行的meta信息
  24. rowData(example_sce)
  25. DataFrame with 2000 rows and 1 column
  26.                        stuff
  27.                    <numeric>
  28. Gene_0001  0.146899100858718
  29. Gene_0002  0.547358682611957
  30. Gene_0003  0.381470382912084
  31. Gene_0004 0.0698823253624141
  32. Gene_0005  0.577666614903137
  33. ...                      ...
  34. Gene_1996  0.810028552776203
  35. Gene_1997   0.92471176572144
  36. Gene_1998   0.73105761455372
  37. Gene_1999  0.496801204746589
  38. Gene_2000  0.135669085429981
  40. # 根据基因的表达过滤掉那些在所有细胞中表达量之和为0的基因
  41. keep_feature <- rowSums(counts(example_sce) > 0) > 0
  42. example_sce <- example_sce[keep_feature,]

对于原始的count表达矩阵，我们也提供了一些函数对其进行数据的归一化和标准化处理。如使用calculateCPM函数计算表达量的CPM(counts-per-million)值，其结果将会存储在SingleCellExperiment对象的“cpm” Assay中，可以通过cpm函数进行访问

  1. cpm(example_sce) <- calculateCPM(example_sce)
  2. head(cpm(example_sce))
  3.           Cell_001   Cell_002  Cell_003 Cell_004 Cell_005 Cell_006
  4. Gene_0001     0.00 749.529259  6.561271    0.000    0.000   0.0000
  5. Gene_0002  1344.85 396.092698  9.841906 5558.424 2826.476 923.5422
  6. Gene_0003     0.00   0.000000  0.000000    0.000 1483.563   0.0000
  7. Gene_0004     0.00   6.093734  0.000000    0.000    0.000 571.4418
  8. Gene_0005     0.00   0.000000 36.086989    0.000    0.000   0.0000
  9. Gene_0006     0.00   0.000000  0.000000    0.000    0.000   0.0000
  10.              Cell_007  Cell_008    Cell_009 Cell_010    Cell_011    Cell_012
  11. Gene_0001    0.000000 69.780424    2.698593   0.0000 9959.085768    5.882457
  12. Gene_0002    6.938109  3.322877    0.000000   0.0000    7.590767    0.000000
  13. Gene_0003    0.000000  3.322877    0.000000   0.0000    7.590767 1047.077301
  14. Gene_0004 1651.269847  0.000000    1.349296 168.5733    0.000000    5.882457
  15. Gene_0005    0.000000  0.000000    1.349296   0.0000    0.000000   11.764913
  16. Gene_0006 2334.673545  0.000000 4174.723091   0.0000    0.000000 1588.263322
  17.              Cell_013 Cell_014   Cell_015   Cell_016   Cell_017    Cell_018
  18. Gene_0001 2013.948828 0.000000    0.00000    2.64154  118.89733    2.069181
  19. Gene_0002 5377.144161 0.000000   11.56233    0.00000 3175.25816 1392.558811
  20. Gene_0003    0.000000 0.000000    0.00000    0.00000    0.00000    0.000000
  21. Gene_0004    0.000000 3.334801    0.00000 2393.23554 2290.52213 2110.564617
  22. Gene_0005    3.968372 0.000000    0.00000    0.00000    0.00000    4.138362
  23. Gene_0006    3.968372 0.000000 3399.32534    0.00000   10.49094    0.000000

同样的，我们也可以使用normalize函数进行数据的归一化处理，它将对原始的count矩阵进行一个log2的转换处理。This is done by dividing each count by its size factor (or scaled library size, if no size factors are defined), adding a pseudo-count and log-transforming. 归一化后的结果存储在”logcounts” Assay中，可以通过logcounts函数进行访问。

  1. # 使用normalize函数进行数据归一化
  2. example_sce <- normalize(example_sce)
  3. # 查看assay的信息
  4. assayNames(example_sce)
  5. [1] "counts"    "cpm"       "logcounts"
  7. head(logcounts(example_sce))

我们可以使用calcAverage函数计算基因的平均表达量

  1. head(calcAverage(example_sce))
  2.  Gene_0001  Gene_0002  Gene_0003  Gene_0004  Gene_0005  Gene_0006 
  3. 305.551749 325.719897 183.090462 162.143201   1.231123 187.167913

使用其他的方法导入基因表达矩阵

  • 对于CSV格式存储的基因表达矩阵，我们可以通过read.table()函数或data.table包中的fread()函数进行读取。
  • 对于一些大型的数据集，在读取的过程中会产生大量的缓存，需要较大的内存，因此我们可以通过Matrix包中的readSparseCounts()函数读取大型数据集，并将其存储为一个稀疏矩阵，可以有效减小系统的读取内存。
  • 对于来自10x Genomics产生的表达矩阵，我们可以通过DropletUtils包中的read10xCounts()函数进行读取，读取后它会自动生成一个SingleCellExperiment对象
  • 对于kallisto和Salmon等比对软件产生的基因表达矩阵，我们可以通过tximeta包中的 readSalmonResults()和readKallistoResults()函数进行读取。

参考来源：
http://www.bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/overview.html

使用scater包进行单细胞测序分析（一）：数据导入与SingleCellExperiment对象构建 - 图5
使用scater包进行单细胞测序分析（一）：数据导入与SingleCellExperiment对象构建 - 图6