20211109（周二）
1.今天准备学习scanpy单细胞分析包，走在路上想，我用seurat和scanpy的意义在哪里，多学一个软件吗？
学习软件后面的代码逻辑 远远大于学会软件，后面抽时间一定要把seurat和scanpy的主要函数功能和代码底层逻辑弄清楚，
而不是盲目跑程序，不要只是工具包的使用者；
2.回来问XXX，学习的方法是什么？他说学习没有捷径而言，只有量变到质变的转化，你才能深入思考；

前面是之前的随记。于是就决定把seurat源码学习一下，尽量不要只做调“黑匣子”的工具人。

阅读seurat包的源码并非易事，seurat包是面向对象的编程思维，代码中有大量的S3，S4基于泛型函数（generic function）的调用。对于习惯面向过程编程的初学者，不易阅读和理解。代码中除了调用SeuratObject，sctransform等R包外，seurat引用了Rcpp调用C++函数，提高运算性能，不少函数还支持并行运算。就像做英语阅读理解一样，一半靠读，一半靠猜，猜不透就实操。读到后面发现底层都是数理统计和统计建模等知识，要不断的查资料，一点点消化，很难读通透，优雅也只是个愿景。
也不得不佩服开发seurat包的大师们功底深厚！在读的过程中，一些经验也可以慢慢记录和积攒，这份记录后面可能会修改，前面理解的不到位。
seurat的github地址：https://github.com/satijalab/seurat

1. 如何方便的查阅源码

我们一般去github网站seurat储存库查阅源码，比如想查看FindMarkers函数的代码，在github搜索框中输入“FindMarkers”，选择仅在此储存库中搜索即可。

image.png

image.png

image.png

image.png

2. R语言与面向对象编程

我们先对面向对象编程有个粗略的认识。
每种编程语言都有面向对象和面向过程的编程方式，非R语言独创。编程语言也是相互学习相互借鉴，尽量弥补自身的不足，R语言运行慢就调用C++。R语言虽然被认为是一种函数式语言, 但同样支持面向对象编程。
其实无论是面向对象还是面向过程，都是我们在编程时解决问题的一种思维方式。

面向过程（Procedure Oriented，简称PO）：当解决一个问题的时候，我们一般会把所需要的步骤列出来，然后通过计算机中的函数把这些步骤逐一实现，这种过程化的叙事思维，就是面向过程思想。我们会把复用的代码提炼出来，写成公共函数，但没有类的概念。
面向对象（Object Oriented，简称OO）：面向对象编程中有两个重要的术语“类”和“对象”。类是对某个事物的概括定义，可以看做是一个蓝图。对象则是对某个事物的具体实现，可以看做是依照蓝图建立起的房子。当解决一个问题的时候，面向对象会把事物抽象成对象的概念，就是说这个问题里面有哪些对象，然后给对象赋一些属性和方法，然后让每个对象去执行自己的方法，问题得到解决。
面向对象可以基于对象的共性做抽象，程序的特点是模块化和封装，继承、多态性的特性，可以设计出低耦合的系统，显著的改善了软件系统的可维护性，代码易维护、易复用、易扩展。很多大型软件开发选型上，大多会使用面向对象语言编程。

一个面向对象系统的核心是其实现的类 (class) 和方法(method)

• 类定义了对象共同拥有的某些属性及其与其他类的关系
• 如果一个对象属于某个类，则称该对象是这个类的实例(instance)
• 方法是一种与特定类的对象关联的函数

R语言的面向对象跟其他Java，php等语言不是很相似，它是基于泛型函数(generic function)的，而不是基于类层次结构。
目前R中最热门的OOP系统主要有四种：S3,S4,R6, Reference Class。

问题：什么是泛型函数？
R语言是一种函数式语言，假如需要编写一个concat函数；
如果我们的数据是数组，我们就编写一个array_concat函数；
如果我们的数据是列表，我们就编写一个list_concat函数；
如果需求改变，需要写一个统一的函数，对数组和列表都能连接，于是有了concat函数。代码如下：

array_concat <- function(x,y) { 
  array(c(x, y))
}

list_concat <- function(x,y) {
  list(c(unlist(x), unlist(y)))
}

concat <- function(x, y) {
  if (is.array(x) && is.array(y)){
    return(array_concat(x, y))
  } else if (is.list(x) && is.list(y)) {
    return(list_concat(x, y))
  } else {
    stop("not supported type")
  }
}

但是，我们的需求一直在变化，如果新增不同的输入数据类型，上面的函数有得重写；每添加一种新的数据类型，concat函数都要修改，甚至重写，不便于代码的维护，可读性也很差，代码量很大，也不便于管理。于是，就有了开发者想开发一个“可扩展”的函数，在不改变原有代码的情况下，对函数能添加新的数据类型和属性。

R语言就有了泛型函数（Generic Function）的出现。
泛型函数是一个函数族，其中的每个函数都有相似的功能，实现对不同对象使用不同的方法。与一般函数相比，泛型函数多了一个判断对象类并传到指定函数的过程。例如concat.array处理array对象，而concat.list处理list对象。

3. 基于S3的面向对象编程

S3是目前最容易使用的方式，特别是处理R包时，使代码更易于理解和组织。

• S3实现了一种基于泛型函数的面向对象方式
• 泛型函数是一种特殊的函数, 其根据传入对象的类型决定调用哪个具体的方法

S3函数的创建
S3对象组成：generic(generic FUN)+method(generic.class FUN）
通常用UseMethod()函数定义一个泛型函数的名称，通过传入参数的class属性，来确定对应的方法调用。

定义一个get_n_elements的泛型函数

• 用UseMethod()定义get_n_elements泛型函数
• 用get_n_elements.xxx的语法格式定义get_n_elements对象的行为
• get_n_elements.default是默认行为

get_n_elements <- function(x,...)
{
    UseMethod("get_n_elements")
}

method(generic.class）函数，创建示例：

# Create a data.frame method for get_n_elements
get_n_elements.data.frame <- function(x, ...)
{
    nrow(x) * ncol(x) # or prod(dim(x))
}
 
# Create a default method for get_n_elements
# 在使用UseMethod调用时，先在methods中寻找对应class，如果都没有找到，则会调用default方法。
get_n_elements.default <- function(x,...)
{
    length(unlist(x))
}

4.seurat包的S3泛型函数

切换到seurat包，我们可以看到大量的S3泛型函数。

在seurat的NAMESPACE文件可以看到大量的S3泛型函数及方法，也可以通过R语句查看。

image.png

# List Methods for S3 Generic Functions or Classes
methods(class="Seurat")

# S3 class
.S3methods(class="Seurat")

# 用pryr包中ftype()函数，检查FindMarkers类型
pryr::ftype(FindMarkers)
# [1] "s3"      "generic"

methods(Seurat::FindMarkers)
# [1] FindMarkers.Assay*    FindMarkers.default*  FindMarkers.DimReduc* FindMarkers.Seurat*  

pryr::is_s3_generic("FindMarkers")
# [1] TRUE

# getS3method(f, class) 显示泛型函数对特定对象方法的实现
getS3method("FindMarkers", "Seurat")

例如FindMarkers函数就是一个S3泛型函数；

# NAMESPACE
S3method(FindMarkers,Assay)
S3method(FindMarkers,DimReduc)
S3method(FindMarkers,Seurat)
S3method(FindMarkers,default)

FindMarkers <- function(object, ...) {
  UseMethod(generic = 'FindMarkers', object = object)
}

# differential_expression.R
FindMarkers.default <- function(){...}
FindMarkers.Assay <- function(){...}
FindMarkers.DimReduc <- function(){...}
FindMarkers.Seurat <- function(){...}

案例：我们有1个seurat对象seurat_obj，查找cluster0和cluster1的差异基因。

DefaultAssay(seurat_obj) <- "RNA" 
DEG_wilcox <- FindMarkers(seurat_obj, ident.1 = 0, ident.2= 1,  test.use = "wilcox", min.pct = 0.5, logfc.threshold = 0.5, only.pos = T) 
head(DEG_wilcox) 
#           p_val avg_log2FC pct.1 pct.2 p_val_adj 
# RBP7         0   1.698551 0.726 0.004         0 
# PGD          0   1.670035 0.823 0.046         0 
# AGTRAP       0   1.488062 0.833 0.136         0 
# TNFRSF1B     0   1.696910 0.858 0.119         0 
# EFHD2        0   2.097743 0.917 0.079         0 
# CDA          0   1.577383 0.760 0.001         0 
dim(DEG_wilcox) 
# [1] 789   5

class(seurat_obj)
# [1] "Seurat"
# attr(,"package")
# [1] "SeuratObject"

问题：我们该如何一层层去看FindMarkers代码？
说明：FindMarkers的具体解析会在《seurat-FindAllMarkers()源码解析》说明，我们只摘出重要的代码：
step1：通过判断class(seurat_obj)，我们知道seurat_obj是seurat对象，先调用FindMarkers.Seurat()函数，将FindMarkers的参数传至FindMarkers.Seurat。

object <- seurat_obj 
ident.1 = 0 
ident.2= 1 
test.use = "wilcox" 
min.pct = 0.5 
logfc.threshold = 0.5 
only.pos = T

step2：FindMarkers.Seurat()代码中提取归一化的基因表达矩阵seurat_obj[['RNA']]@data，赋值给data.use，基因表达矩阵大小为36601 *8824 ，36601个基因，8824个细胞。计算一些具体参数后，继续传递给FindMarkers函数；
我们判断class(data.use)，data.use为Assay，我们调用FindMarkers.Assay()，我们将前面获得的参数继续传递给FindMarkers.Assay。

# select which data to use 
assay <- assay %||% DefaultAssay(object = object) 
assay <- "RNA" 
data.use <- object[[assay]] 
dim(data.use)  
# [1] 36601  8824

...
  de.results <- FindMarkers( 
    object = data.use, 
    slot = slot, 
    cells.1 = cells$cells.1, 
    cells.2 = cells$cells.2, 
    features = features, 
    logfc.threshold = logfc.threshold, 
    test.use = test.use, 
    min.pct = min.pct, 
    min.diff.pct = min.diff.pct, 
    verbose = verbose, 
    only.pos = only.pos, 
    max.cells.per.ident = max.cells.per.ident, 
    random.seed = random.seed, 
    latent.vars = latent.vars, 
    min.cells.feature = min.cells.feature, 
    min.cells.group = min.cells.group, 
    pseudocount.use = pseudocount.use, 
    mean.fxn = mean.fxn, 
    base = base, 
    fc.name = fc.name, 
    densify = densify, 
    ... 
  ) 
  return(de.results)

class(data.use)
# [1] "Assay" 
# attr(,"package") 
# [1] "SeuratObject"

step3：我们就这样一层层查看传入对象的类型，决定调用哪个具体的方法，最后我们进入到FindMarkers.default()这个函数，完成对FindMarkers函数源码的解析。
我们熟知S3泛型函数的一些基本功能，就能顺藤摸瓜，一步步去解析seurat包的S3函数，也让我们的理解变得更清晰。