20211109（周二）
1.今天準(zhǔn)備學(xué)習(xí)scanpy單細(xì)胞分析包，走在路上想，我用seurat和scanpy的意義在哪里，多學(xué)一個軟件嗎？
學(xué)習(xí)軟件后面的代碼邏輯 遠(yuǎn)遠(yuǎn)大于學(xué)會軟件游沿，后面抽時間一定要把seurat和scanpy的主要函數(shù)功能和代碼底層邏輯弄清楚，
而不是盲目跑程序肮砾，不要只是工具包的使用者诀黍；
2.回來問XXX，學(xué)習(xí)的方法是什么仗处？他說學(xué)習(xí)沒有捷徑而言眯勾，只有量變到質(zhì)變的轉(zhuǎn)化，你才能深入思考婆誓；

前面是之前的隨記吃环。于是就決定把seurat源碼學(xué)習(xí)一下，盡量不要只做調(diào)“黑匣子”的工具人洋幻。

閱讀seurat包的源碼并非易事郁轻，seurat包是面向?qū)ο蟮木幊趟季S，代碼中有大量的S3文留，S4基于泛型函數(shù)（generic function）的調(diào)用好唯。對于習(xí)慣面向過程編程的初學(xué)者，不易閱讀和理解厂庇。代碼中除了調(diào)用SeuratObject渠啊，sctransform等R包外输吏，seurat引用了Rcpp調(diào)用C++函數(shù)权旷，提高運算性能，不少函數(shù)還支持并行運算贯溅。就像做英語閱讀理解一樣拄氯，一半靠讀，一半靠猜它浅，猜不透就實操译柏。讀到后面發(fā)現(xiàn)底層都是數(shù)理統(tǒng)計和統(tǒng)計建模等知識，要不斷的查資料姐霍，一點點消化鄙麦，很難讀通透典唇，優(yōu)雅也只是個愿景。
也不得不佩服開發(fā)seurat包的大師們功底深厚胯府！在讀的過程中介衔，一些經(jīng)驗也可以慢慢記錄和積攢，這份記錄后面可能會修改骂因，前面理解的不到位炎咖。
seurat的github地址：https://github.com/satijalab/seurat

1. 如何方便的查閱源碼

我們一般去github網(wǎng)站seurat儲存庫查閱源碼，比如想查看FindMarkers函數(shù)的代碼寒波，在github搜索框中輸入“FindMarkers”乘盼，選擇僅在此儲存庫中搜索即可。

image.png

image.png

image.png

image.png

2. R語言與面向?qū)ο缶幊?/h4>

我們先對面向?qū)ο缶幊逃袀€粗略的認(rèn)識。
每種編程語言都有面向?qū)ο蠛兔嫦蜻^程的編程方式十减，非R語言獨創(chuàng)栈幸。編程語言也是相互學(xué)習(xí)相互借鑒，盡量彌補自身的不足帮辟，R語言運行慢就調(diào)用C++速址。R語言雖然被認(rèn)為是一種函數(shù)式語言, 但同樣支持面向?qū)ο缶幊獭?br> 其實無論是面向?qū)ο筮€是面向過程，都是我們在編程時解決問題的一種思維方式由驹。

面向過程（Procedure Oriented芍锚，簡稱PO）：當(dāng)解決一個問題的時候，我們一般會把所需要的步驟列出來蔓榄，然后通過計算機中的函數(shù)把這些步驟逐一實現(xiàn)并炮，這種過程化的敘事思維，就是面向過程思想甥郑。我們會把復(fù)用的代碼提煉出來逃魄，寫成公共函數(shù)，但沒有類的概念澜搅。
面向?qū)ο螅∣bject Oriented伍俘，簡稱OO）：面向?qū)ο缶幊讨杏袃蓚€重要的術(shù)語“類”和“對象”邪锌。類是對某個事物的概括定義，可以看做是一個藍(lán)圖癌瘾。對象則是對某個事物的具體實現(xiàn)秃流，可以看做是依照藍(lán)圖建立起的房子。當(dāng)解決一個問題的時候柳弄，面向?qū)ο髸咽挛锍橄蟪蓪ο蟮母拍畈罢停褪钦f這個問題里面有哪些對象，然后給對象賦一些屬性和方法碧注，然后讓每個對象去執(zhí)行自己的方法嚣伐，問題得到解決。
面向?qū)ο罂梢曰趯ο蟮墓残宰龀橄笃钾ぃ绦虻奶攸c是模塊化和封裝轩端，繼承、多態(tài)性的特性逝变，可以設(shè)計出低耦合的系統(tǒng)基茵，顯著的改善了軟件系統(tǒng)的可維護(hù)性，代碼易維護(hù)壳影、易復(fù)用拱层、易擴展。很多大型軟件開發(fā)選型上宴咧，大多會使用面向?qū)ο笳Z言編程根灯。

一個面向?qū)ο笙到y(tǒng)的核心是其實現(xiàn)的類 (class) 和方法(method)

• 類定義了對象共同擁有的某些屬性及其與其他類的關(guān)系
• 如果一個對象屬于某個類，則稱該對象是這個類的實例(instance)
• 方法是一種與特定類的對象關(guān)聯(lián)的函數(shù)

R語言的面向?qū)ο蟾渌鸍ava掺栅，php等語言不是很相似烙肺，它是基于泛型函數(shù)(generic function)的，而不是基于類層次結(jié)構(gòu)氧卧。
目前R中最熱門的OOP系統(tǒng)主要有四種：S3,S4,R6, Reference Class桃笙。

問題：什么是泛型函數(shù)？
R語言是一種函數(shù)式語言沙绝，假如需要編寫一個concat函數(shù)搏明；
如果我們的數(shù)據(jù)是數(shù)組，我們就編寫一個array_concat函數(shù)宿饱；
如果我們的數(shù)據(jù)是列表熏瞄，我們就編寫一個list_concat函數(shù)脚祟；
如果需求改變谬以，需要寫一個統(tǒng)一的函數(shù)，對數(shù)組和列表都能連接由桌，于是有了concat函數(shù)为黎。代碼如下：

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")
  }
}

但是邮丰，我們的需求一直在變化，如果新增不同的輸入數(shù)據(jù)類型铭乾，上面的函數(shù)有得重寫剪廉；每添加一種新的數(shù)據(jù)類型，concat函數(shù)都要修改炕檩，甚至重寫斗蒋，不便于代碼的維護(hù)，可讀性也很差笛质，代碼量很大泉沾，也不便于管理。于是妇押，就有了開發(fā)者想開發(fā)一個“可擴展”的函數(shù)跷究，在不改變原有代碼的情況下，對函數(shù)能添加新的數(shù)據(jù)類型和屬性敲霍。

R語言就有了泛型函數(shù)（Generic Function）的出現(xiàn)俊马。
泛型函數(shù)是一個函數(shù)族，其中的每個函數(shù)都有相似的功能肩杈，實現(xiàn)對不同對象使用不同的方法柴我。與一般函數(shù)相比，泛型函數(shù)多了一個判斷對象類并傳到指定函數(shù)的過程扩然。例如concat.array處理array對象屯换，而concat.list處理list對象。

3. 基于S3的面向?qū)ο缶幊?/h4>

S3是目前最容易使用的方式与学，特別是處理R包時彤悔，使代碼更易于理解和組織。

• S3實現(xiàn)了一種基于泛型函數(shù)的面向?qū)ο蠓绞?/li>
• 泛型函數(shù)是一種特殊的函數(shù), 其根據(jù)傳入對象的類型決定調(diào)用哪個具體的方法

S3函數(shù)的創(chuàng)建
S3對象組成：generic(generic FUN)+method(generic.class FUN）
通常用UseMethod()函數(shù)定義一個泛型函數(shù)的名稱索守，通過傳入?yún)?shù)的class屬性晕窑，來確定對應(yīng)的方法調(diào)用。

定義一個get_n_elements的泛型函數(shù)

• 用UseMethod()定義get_n_elements泛型函數(shù)
• 用get_n_elements.xxx的語法格式定義get_n_elements對象的行為
• get_n_elements.default是默認(rèn)行為

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

method(generic.class）函數(shù)卵佛，創(chuàng)建示例：

# 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調(diào)用時杨赤，先在methods中尋找對應(yīng)class，如果都沒有找到截汪，則會調(diào)用default方法疾牲。
get_n_elements.default <- function(x,...)
{
    length(unlist(x))
}

4.seurat包的S3泛型函數(shù)

切換到seurat包，我們可以看到大量的S3泛型函數(shù)衙解。

在seurat的NAMESPACE文件可以看到大量的S3泛型函數(shù)及方法阳柔，也可以通過R語句查看。

image.png

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

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

# 用pryr包中ftype()函數(shù)蚓峦，檢查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) 顯示泛型函數(shù)對特定對象方法的實現(xiàn)
getS3method("FindMarkers", "Seurat")

例如FindMarkers函數(shù)就是一個S3泛型函數(shù)舌剂；

# 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對象避消，先調(diào)用FindMarkers.Seurat()函數(shù)低滩，將FindMarkers的參數(shù)傳至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()代碼中提取歸一化的基因表達(dá)矩陣seurat_obj[['RNA']]@data岩喷，賦值給data.use吹埠，基因表達(dá)矩陣大小為36601 *8824 盖灸，36601個基因菲盾，8824個細(xì)胞竿屹。計算一些具體參數(shù)后，繼續(xù)傳遞給FindMarkers函數(shù)妇穴；
我們判斷class(data.use)爬虱，data.use為Assay，我們調(diào)用FindMarkers.Assay()腾它，我們將前面獲得的參數(shù)繼續(xù)傳遞給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：我們就這樣一層層查看傳入對象的類型，決定調(diào)用哪個具體的方法瞒滴，最后我們進(jìn)入到FindMarkers.default()這個函數(shù)曲梗，完成對FindMarkers函數(shù)源碼的解析。
我們熟知S3泛型函數(shù)的一些基本功能妓忍，就能順藤摸瓜虏两，一步步去解析seurat包的S3函數(shù)，也讓我們的理解變得更清晰世剖。