天子呼來不上船苫拍，
自稱臣是菜鳥團瓜富。
在這里，和國際同行一起學習單細胞數(shù)據(jù)分析纽什。
我今天并沒有敲滿100行代碼措嵌，以至于我在寫這篇文章的時候都有點不好意思了。要讀/寫長代碼的想法腦海里輪回了很久了芦缰，哪來的時間呢铅匹？寫吧。
我們?yōu)槭裁匆獙W習長代碼?
對大學畢業(yè)才開始學習一門語言的人來講饺藤，學習基本概念并不是一件難事包斑。Perl語言，聽過吧涕俗，有那么多的操作符和參數(shù)罗丰，在不和C語言來比較的前提下，這并不是一個順從人類本能的語言再姑。
知道一門語言的基本概念是一回事萌抵，能寫出一手好程序是另一回事。
Seurat Weekly 目前已經(jīng)進行到第十四期了，之前我們讀了它的很多幫助文檔绍填，這一次我們來試著讀一讀它的源代碼霎桅。借此機會呢，也回答兩個有趣的問題讨永。
首先滔驶，我們把Seurat的源代碼clone 到 本地便于我們查找和閱讀。
git clone https://github.com/satijalab/seurat.git
然后我們就要在結構上看看Seurat的框架
tree -L 1
.
|-- CODE_OF_CONDUCT.md
|-- DESCRIPTION
|-- LICENSE
|-- NAMESPACE
|-- NEWS.md
|-- R
|-- README.md
|-- _pkgdown.yaml
|-- appveyor.yml
|-- azure-pipelines.yml
|-- cran-comments.md
|-- data
|-- index.md
|-- inst
|-- man
|-- seurat.Rproj
|-- src
|-- tests
|-- travis_setup.sh
`-- vignettes
這時候就發(fā)揮程序員最能讀文檔的功能了卿闹，先把眼前看到的文檔讀一遍揭糕，至少要cat 一番，知道它是如何架構的锻霎。然后是找出除了配置文件之外最重要的一個文件著角。也許在之前我們非常想看vignettes，也就是教程旋恼，但是教程我們已經(jīng)寫的太多了吏口，先把它往后排一點。這里我們先看 R 這個文件夾：
.
|-- RcppExports.R
|-- clustering.R
|-- convenience.R
|-- data.R
|-- differential_expression.R
|-- dimensional_reduction.R
|-- generics.R
|-- integration.R
|-- mixscape.R
|-- objects.R
|-- preprocessing.R
|-- reexports.R
|-- tree.R
|-- utilities.R
|-- visualization.R
|-- zfRenderSeurat.old
`-- zzz.R
我們也看到這里都是點R的文件冰更，也就是R腳本锨侯，這時候我們就需要挨個來看了。當然了冬殃，如果你是比較喜歡畫圖囚痴，可以先學visualization.R,這里有實打實的可視化腳本，腳本太長审葬，我們就不再展示了深滚。不管從哪里開始最終都是要全看的，不是說要背下來涣觉，而是學習人家寫代碼的風格和一些小的技巧痴荐。然而，這里可能回憶道一個閱讀障礙官册，那就是RcppExports.R 這個腳本生兆，和一般的R代碼不同，里面幾乎都是這個格式的代碼：
RunModularityClusteringCpp <- function(SNN, modularityFunction, resolution, algorithm, nRandomStarts, nIterations, randomSeed, printOutput, edgefilename) {
.Call('_Seurat_RunModularityClusteringCpp', PACKAGE = 'Seurat', SNN, modularityFunction, resolution, algorithm, nRandomStarts, nIterations, randomSeed, printOutput, edgefilename)
}
RunUMISampling <- function(data, sample_val, upsample = FALSE, display_progress = TRUE) {
.Call('_Seurat_RunUMISampling', PACKAGE = 'Seurat', data, sample_val, upsample, display_progress)
}
RunUMISamplingPerCell <- function(data, sample_val, upsample = FALSE, display_progress = TRUE) {
.Call('_Seurat_RunUMISamplingPerCell', PACKAGE = 'Seurat', data, sample_val, upsample, display_progress)
}
于是我們就要找到Rcpp到底在哪放著膝宁。在src的目錄下放著：
|-- Makevars
|-- ModularityOptimizer.cpp
|-- ModularityOptimizer.h
|-- RModularityOptimizer.cpp
|-- RcppExports.cpp
|-- data_manipulation.cpp
|-- data_manipulation.h
|-- fast_NN_dist.cpp
|-- integration.cpp
|-- integration.h
|-- snn.cpp
|-- snn.h
`-- valid_pointer.c
想要看懂這些代碼除了四級詞匯還需要一些C語言的基礎知識鸦难。其實真看的話也并沒有那么難，一行一行讀就是了员淫，《紅樓夢》都看了合蔽。
#include "ModularityOptimizer.h"
#include <algorithm>
#include <exception>
#include <functional>
#include <numeric>
#include <stdexcept>
using namespace ModularityOptimizer;
using namespace std::chrono;
JavaRandom::JavaRandom(uint64_t seed) {
setSeed(seed);
}
void JavaRandom::setSeed(uint64_t seed) {
this->seed = (seed ^ uint64_t(0x5DEECE66D)) & ((uint64_t(1) << 48) - 1);
}
int JavaRandom::next(int bits) {
// Only 31 bits ever used.
seed = (seed * uint64_t(0x5DEECE66D) + uint64_t(0xB)) & ((uint64_t(1) << 48) - 1);
return static_cast<int>(seed >> (48 - bits));
}
要把這部分讀完估計是要發(fā)一段時間的，下面我們看開一個比較輕松的模塊vignettes是大家比較熟悉的教程介返，一個教程一個md拴事。
.
|-- archive.Rmd
|-- archive.yaml
|-- assets
|-- atacseq_integration_vignette.Rmd
|-- cell_cycle_vignette.Rmd
|-- conversion_vignette.Rmd
|-- de_vignette.Rmd
|-- dim_reduction_vignette.Rmd
|-- essential_commands.Rmd
|-- extensions.Rmd
|-- future_vignette.Rmd
|-- get_started.Rmd
|-- hashing_vignette.Rmd
|-- install.Rmd
|-- integration_introduction.Rmd
|-- integration_large_datasets.Rmd
|-- integration_mapping.Rmd
|-- integration_rpca.Rmd
|-- interaction_vignette.Rmd
|-- merge_vignette.Rmd
|-- mixscape_vignette.Rmd
|-- multimodal_reference_mapping.Rmd
|-- multimodal_vignette.Rmd
|-- pbmc3k_tutorial.Rmd
|-- sctransform_vignette.Rmd
|-- spatial_vignette.Rmd
|-- v4_changes.Rmd
|-- vignettes.yaml
|-- visualization_vignette.Rmd
`-- weighted_nearest_neighbor_analysis.Rmd
看完本文你也許會納悶：為什么總有人在找教程沃斤，源代碼里面不是滿滿的教程嗎？而且是Rmd的拿出來那是可以直接跑的啊刃宵，到處找人要什么pipeline衡瓶，也值當?shù)摹?/p>
在這里，我們提出兩個問題牲证，大家可以在留言區(qū)回答：
- 為什么Seurat的函數(shù)運行完返回的還是Seurat對象哮针，它讀入的和輸出的都一樣，好奇怪从隆。
- 為什么Seurat的聚類結果中編號是從零開始的，而且缭裆，細胞數(shù)量依次遞減
希望您從源碼中給出你的看法键闺。
