劉小澤寫于19.8.9-第三單元第五講：利用scRNAseq包學(xué)習(xí)scater
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

復(fù)習(xí)

我們之前在第14篇筆記中學(xué)習(xí)了scRNAseq包的數(shù)據(jù)囤官，知道了它可以分成4組：pluripotent stem cells 分化而成的 neural progenitor cells (NPC非迹，神經(jīng)前體細(xì)胞) 四苇，還有 GW16（radial glia，放射狀膠質(zhì)細(xì)胞） 右遭、GW21（newborn neuron，新生兒神經(jīng)元） 阿宅、GW21+3(maturing neuron接谨，成熟神經(jīng)元) ，其中NPC和其他組區(qū)別較大芽唇。并且使用常規(guī)的操作進(jìn)行了初步的探索顾画，但是畢竟是單細(xì)胞數(shù)據(jù)，就會(huì)有專業(yè)的包對(duì)它進(jìn)行處理，于是這次就來介紹第一款：scater包(https://www.bioconductor.org/packages/release/bioc/html/scater.html)

首先來看官網(wǎng)教程

親測(cè)亲雪，6個(gè)小時(shí)就能學(xué)完Scater官方提供的三個(gè)教程

這個(gè)包是EMBL和劍橋大學(xué)發(fā)布的，是為分析單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)而開發(fā)疚膊，它包含了一些特性：

• 它需要利用SingleCellExperiment這個(gè)對(duì)象义辕，就是這個(gè)東西（來自Bioconductor-workshop）：
  
• 可以導(dǎo)入非比對(duì)工具 kallisto and Salmon 得到的定量結(jié)果
• 計(jì)算了大量的QC指標(biāo)，方便過濾
• 可視化方面做得不錯(cuò)寓盗，設(shè)計(jì)了大量的函數(shù)（尤其針對(duì)質(zhì)控）灌砖，并且功能如其名
• 從2017年7月，scater包就改變了整體架構(gòu)傀蚌，從之前的SCESet對(duì)象更改成了更多人都在用的SingleCellExperiment 基显，使用Bioconductor 3.6 (2017.10發(fā)布)安裝的包，都會(huì)是SingleCellExperiment 對(duì)象善炫。如果要更迭也不難撩幽，官方給了許多解決方案，例如：toSingleCellExperiment函數(shù)箩艺、updateSCESet函數(shù)

假設(shè)現(xiàn)在有一個(gè)包含表達(dá)量信息的矩陣窜醉，其中包含的feature信息可以是gene, exon, region, etc.

第一步 創(chuàng)建一個(gè)SingleCellExperiment對(duì)象 (官網(wǎng) 24 May 2019)

• 需要注意的是，官方友情提示艺谆，在導(dǎo)入對(duì)象之前榨惰，最好是將表達(dá)量數(shù)據(jù)存為矩陣；
• 如果是較大的數(shù)據(jù)集静汤，官方建議使用chunk-by-chunk的方法琅催，參考Matrix 包，然后使用readSparseCounts函數(shù)虫给，有效減少內(nèi)存使用量（因?yàn)樗梢圆粚⒋罅康?表達(dá)量放進(jìn)內(nèi)存）藤抡；
• 如果是導(dǎo)入10X的數(shù)據(jù)，使用DropletUtils 包的read10xCounts函數(shù)即可抹估，它會(huì)自動(dòng)生成一個(gè)SingleCellExperiment對(duì)象
• 對(duì)于非比對(duì)定量工具杰捂，scater也提供了readSalmonResults、readKallistoResults 支持兩款軟件棋蚌，它的背后利用的是tximport

rm(list = ls()) 
Sys.setenv(R_MAX_NUM_DLLS=999) ##在R3.3版本中嫁佳，只能有100個(gè)固定的動(dòng)態(tài)庫(kù)限制，到了3.4版本以后谷暮，就能夠使用Sys.setenv(R_MAX_NUM_DLLS=xxx)進(jìn)行設(shè)置蒿往，而這個(gè)數(shù)字根據(jù)個(gè)人情況設(shè)定
options(stringsAsFactors = F) 

# 使用包自帶的測(cè)試數(shù)據(jù)進(jìn)行操作
library(scater)
data("sc_example_counts")
data("sc_example_cell_info")

> dim(sc_example_counts)
[1] 2000   40
> sc_example_counts[1:3,1:3] # 基因表達(dá)矩陣
          Cell_001 Cell_002 Cell_003
Gene_0001        0      123        2
Gene_0002      575       65        3
Gene_0003        0        0        0

> sc_example_cell_info[1:3,1:3] # 樣本注釋信息
             Cell Mutation_Status Cell_Cycle
Cell_001 Cell_001        positive          S
Cell_002 Cell_002        positive         G0
Cell_003 Cell_003        negative         G1

# 接下來就是構(gòu)建對(duì)象(日后只需要復(fù)制粘貼替換即可)
example_sce <- SingleCellExperiment(
  assays = list(counts = sc_example_counts), 
  colData = sc_example_cell_info
)

> example_sce
class: SingleCellExperiment 
dim: 2000 40 
metadata(0):
assays(1): counts
rownames(2000): Gene_0001 Gene_0002 ... Gene_1999 Gene_2000
rowData names(0):
colnames(40): Cell_001 Cell_002 ... Cell_039 Cell_040
colData names(4): Cell Mutation_Status Cell_Cycle Treatment
reducedDimNames(0):
spikeNames(0):

注意到上面構(gòu)建對(duì)象時(shí)使用了counts = sc_example_counts這么一個(gè)定義，官方也推薦湿弦，使用counts作為導(dǎo)入表達(dá)矩陣的名稱瓤漏，這樣會(huì)方面下面的counts函數(shù)提取；另外還支持exprs蔬充、tmp蝶俱、cpm、fpkm這樣的輸入名稱

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

調(diào)用或修改行或列的metadata比較方便：

# 默認(rèn)調(diào)用/修改 列饥漫，所以example_sce$whee就是新增一列metadata
example_sce$whee <- sample(LETTERS, ncol(example_sce), replace=TRUE)
colData(example_sce)

# 如果對(duì)行新增一行metadata(注意這里rowData和原來的矩陣沒有關(guān)系榨呆，它操作的是一些注釋信息)
rowData(example_sce)$stuff <- runif(nrow(example_sce))
rowData(example_sce)

除此以外，還有一些比較復(fù)雜的函數(shù)：例如isSpike 對(duì)spike-in操作庸队，sizeFactors 是進(jìn)行標(biāo)準(zhǔn)化時(shí)對(duì)細(xì)胞文庫(kù)大小計(jì)算的結(jié)果积蜻、reducedDim 對(duì)降維結(jié)果（reduced dimensionality results）操作

另外，對(duì)這個(gè)對(duì)象取子集也是很方便的彻消，例如要過濾掉在所有細(xì)胞中都不表達(dá)的基因：

# 過濾不表達(dá)基因
keep_feature <- rowSums(counts(example_sce) > 0) > 0
example_sce <- example_sce[keep_feature,]

第二步 (可選)計(jì)算一堆表達(dá)統(tǒng)計(jì)值 (官網(wǎng) 24 May 2019)

如果要計(jì)算CPM值竿拆，之前一直使用log2(edgeR::cpm(dat)+1)進(jìn)行計(jì)算，這個(gè)包自己做了一個(gè)函數(shù)：

# 計(jì)算的CPM值存到example_sce對(duì)象的標(biāo)準(zhǔn)命名（cpm）中去
cpm(example_sce) <- calculateCPM(example_sce)

另外還可以提供歸一化：normalize函數(shù)宾尚，它計(jì)算得到：log2-transformed normalized expression values

# 總體計(jì)算方法是：dividing each count by its size factor (or scaled library size, if no size factors are defined), adding a pseudo-count and log-transforming (翻譯一下：將每個(gè)count值除以size factor丙笋，記得之前edgeR進(jìn)行標(biāo)準(zhǔn)化就計(jì)算了這么一個(gè)值，它就是為了均衡各個(gè)樣本文庫(kù)差異煌贴；如果沒有size factor不见，也可以對(duì)文庫(kù)大小進(jìn)行歸一化)，接著加一個(gè)值(例如1崔步，為了不讓log為難)稳吮，最后log計(jì)算
example_sce <- normalize(example_sce) # 結(jié)果保存在logcounts中
assayNames(example_sce)
## [1] "counts"    "cpm"       "logcounts"

注意：表達(dá)矩陣的標(biāo)準(zhǔn)命名中，exprs和logcounts是同義詞井濒，它是為了和老版本的scater函數(shù)兼容：

identical(exprs(example_sce), logcounts(example_sce))
## [1] TRUE

另外灶似，我們也可以根據(jù)需要?jiǎng)?chuàng)建一個(gè)和原始count矩陣同樣維度的新矩陣，存儲(chǔ)在assay中

# 比如創(chuàng)建一個(gè)判斷的矩陣瑞你，看看原來count矩陣中的每個(gè)值是不是都大于0酪惭，結(jié)果是一堆的邏輯值
assay(example_sce, "is_expr") <- counts(example_sce)>0

還有，calcAverage函數(shù)可以計(jì)算樣本歸一化以后者甲，各個(gè)基因的平均表達(dá)量(如果樣本還沒進(jìn)行標(biāo)準(zhǔn)化春感，那么它首先會(huì)計(jì)算size factor)

head(calcAverage(example_sce))
##  Gene_0001  Gene_0002  Gene_0003  Gene_0004  Gene_0005  Gene_0006 
## 305.551749 325.719897 183.090462 162.143201   1.231123 187.167913

第三步 數(shù)據(jù)可視化(官網(wǎng) 24 May 2019)

重點(diǎn)包含這幾方面：

• plotExpression ：畫出一個(gè)或多個(gè)基因在細(xì)胞中的表達(dá)量水平
• plotReducedDim：(計(jì)算或)繪制降維后的坐標(biāo)
• 其他的QC圖

# 創(chuàng)建并標(biāo)準(zhǔn)化
library(scater)
data("sc_example_counts")
data("sc_example_cell_info")
example_sce <- SingleCellExperiment(
    assays = list(counts = sc_example_counts),
    colData = sc_example_cell_info
) 
example_sce <- normalize(example_sce)
example_sce
## class: SingleCellExperiment 
## dim: 2000 40 
## metadata(1): log.exprs.offset
## assays(2): counts logcounts
## rownames(2000): Gene_0001 Gene_0002 ... Gene_1999 Gene_2000
## rowData names(0):
## colnames(40): Cell_001 Cell_002 ... Cell_039 Cell_040
## colData names(4): Cell Mutation_Status Cell_Cycle Treatment
## reducedDimNames(0):
## spikeNames(0):

先來繪制基因表達(dá)相關(guān)的圖：

通過這種方式，可以方便檢查差異分析虏缸、擬時(shí)序分析等等的feature結(jié)果鲫懒；它默認(rèn)使用標(biāo)準(zhǔn)化后的logcounts值，當(dāng)然刽辙，可以定義exprs_values參數(shù)來修改

# 最簡(jiǎn)單的圖
# 默認(rèn)使用標(biāo)準(zhǔn)化后的logcounts值
plotExpression(example_sce, rownames(example_sce)[1:6])

# 增加分組信息：定義x是一個(gè)離散型分組變量
plotExpression(example_sce, rownames(example_sce)[1:6],
    x = "Mutation_Status", exprs_values = "logcounts") 

# 查看繪制的x這個(gè)metadata
> colData(example_sce)$Mutation_Status
 [1] "positive" "positive" "negative" "negative" "negative" "negative" "positive" "positive"
 [9] "negative" "positive" "negative" "negative" "positive" "negative" "negative" "negative"
[17] "positive" "negative" "negative" "negative" "positive" "positive" "negative" "positive"
[25] "negative" "positive" "positive" "negative" "positive" "negative" "negative" "positive"
[33] "positive" "negative" "positive" "negative" "negative" "negative" "negative" "negative"

不難看到窥岩，它是將細(xì)胞先根據(jù)x = "Mutation_Status"進(jìn)行了分組，然后對(duì)6個(gè)基因繪制了表達(dá)分布圖

這個(gè)x參數(shù)的設(shè)置很講究：它的英文含義是 covariate to be shown on the x-axis宰缤，定義了x軸上的協(xié)變量颂翼。簡(jiǎn)單理解晃洒，就是x軸上按照什么來定義，如果x是一個(gè)分類的離散型變量(比如這里的positive朦乏、negative)球及，那么x軸就是為了分組，結(jié)果就是小提琴圖呻疹；如果x是一個(gè)連續(xù)的變量(比如下面??要演示的某個(gè)基因表達(dá)量)吃引，那么x軸就是為了看數(shù)值的變化，結(jié)果就是散點(diǎn)圖

# 定義x是一個(gè)連續(xù)型變量
plotExpression(example_sce, rownames(example_sce)[1:6],
    x = "Gene_0001")

可以自定義顏色诲宇、形狀际歼、大小的區(qū)分惶翻，例如：

plotExpression(example_sce, rownames(example_sce)[1:6],
               colour_by = "Cell_Cycle", shape_by = "Mutation_Status", 
               size_by = "Gene_0002")
# 利用兩個(gè)metadata：Cell_Cycle(區(qū)分顏色)姑蓝、Mutation_Status(區(qū)分形狀)
# 利用一個(gè)表達(dá)量指標(biāo)：Gene_0002(區(qū)分大小)

上面這張圖可以再加上之前說的x參數(shù)，還是按照Mutation_Status進(jìn)行x軸分組

# 添加中位線吕粗、x軸分組
plotExpression(example_sce, rownames(example_sce)[7:12],
    x = "Mutation_Status", exprs_values = "counts", 
    colour = "Cell_Cycle", show_median = TRUE, 
    xlab = "Mutation Status", log = TRUE)

再來繪制降維相關(guān)的圖：

SingleCellExperiment對(duì)象中包含了reducedDims接口纺荧，其中存儲(chǔ)了細(xì)胞降維后的坐標(biāo)，可以用reducedDim颅筋、reducedDims函數(shù)獲取

關(guān)于這兩個(gè)函數(shù)的不同：使用?reducedDim就能獲得
For reducedDim, a numeric matrix is returned containing coordinates for cells (rows) and dimensions (columns).

For reducedDims, a named SimpleList of matrices is returned, with one matrix for each type of dimensionality reduction method.

降維首先利用PCA方法：

# runPCA結(jié)果保存在sce對(duì)象的PCA中宙暇。默認(rèn)情況下，runPCA會(huì)根據(jù)500個(gè)變化差異最顯著的feature的log-count值進(jìn)行計(jì)算议泵，當(dāng)然這個(gè)數(shù)量可以通過ntop參數(shù)修改占贫。
example_sce <- runPCA(example_sce) 

> reducedDimNames(example_sce)
[1] "PCA"

> example_sce
class: SingleCellExperiment 
dim: 2000 40 
metadata(1): log.exprs.offset
assays(2): counts logcounts
rownames(2000): Gene_0001 Gene_0002 ... Gene_1999 Gene_2000
rowData names(0):
colnames(40): Cell_001 Cell_002 ... Cell_039 Cell_040
colData names(4): Cell Mutation_Status Cell_Cycle Treatment
reducedDimNames(1): PCA
spikeNames(0):

注：runPCA源代碼在此：https://rdrr.io/bioc/scater/src/R/runPCA.R

任何的降維結(jié)果，都能用plotReducedDim函數(shù)作圖

plotReducedDim(example_sce, use_dimred = "PCA", 
    colour_by = "Treatment", shape_by = "Mutation_Status")

然后可以利用表達(dá)量添加顏色先口、大小分組

# 看特定基因在PCA過程中起到的作用
plotReducedDim(example_sce, use_dimred = "PCA", 
    colour_by = "Gene_1000", size_by = "Gene_0500")

除了使用plotReducedDim型奥，還能使用plotPCA自己去生成，但前提還是要使用PCA的計(jì)算結(jié)果碉京；如果檢測(cè)不到PCA的計(jì)算坐標(biāo)厢汹，這個(gè)函數(shù)會(huì)自己runPCA計(jì)算一遍。盡管如此谐宙，還是推薦先進(jìn)行計(jì)算烫葬，再作圖。因?yàn)橛袝r(shí)候我們需要利用一個(gè)數(shù)據(jù)集做出多張不同的圖(就像上面??一樣)凡蜻，但是每做一張圖這個(gè)函數(shù)都要運(yùn)行一遍太費(fèi)時(shí)間搭综，如果已經(jīng)計(jì)算好，那么它就能直接調(diào)用划栓，十分方便

# 最簡(jiǎn)單的plotPCA
plotPCA(example_sce)

它也可以像plotReducedDim一樣设凹，定義顏色、大小

plotPCA(example_sce, colour_by = "Gene_0001", size_by = "Gene_1000")

另外我們可以自己選擇進(jìn)行PCA的數(shù)據(jù)茅姜，例如使用自己的feature_control（例如使用的ERCC spike-in ）闪朱，來看看數(shù)據(jù)中是否存在技術(shù)誤差而導(dǎo)致差異

# 默認(rèn)情況下月匣，runPCA會(huì)根據(jù)500個(gè)變化差異最顯著的feature。這里定義 feature_set可以覆蓋默認(rèn)設(shè)置
# 看官方描述：eature_set Character vector of row names, a logical vector or a numeric vector of indices indicating a set of features to use for PCA. This will override any \code{ntop} argument if specified.

example_sce2 <- runPCA(example_sce, 
    feature_set = rowData(example_sce)$is_feature_control)
plotPCA(example_sce2)

還可以繪制多個(gè)主成分：

example_sce <- runPCA(example_sce, ncomponents=20)
# 繪制4個(gè)成分
plotPCA(example_sce, ncomponents = 4, colour_by = "Treatment",
        shape_by = "Mutation_Status")

接著使用t-SNE(t-distributed stochastic neighbour embedding)降維：**

關(guān)于tsne這個(gè)流行的算法奋姿，有必要了解一下：

• tsne的作者Laurens強(qiáng)調(diào)锄开，可以通過t-SNE的可視化圖提出一些假設(shè)，但是不要用t-SNE來得出一些結(jié)論称诗，想要驗(yàn)證你的想法萍悴，最好用一些其他的辦法。
• t-SNE中集群之間的距離并不表示相似度 寓免，同一個(gè)數(shù)據(jù)上運(yùn)行t-SNE算法多次癣诱，很有可能得到多個(gè)不同“形態(tài)”的集群。但話說回來袜香，真正有差異的群體之間撕予，不管怎么變換形態(tài)，它們還是有差別
• 關(guān)于perplexity的使用：(默認(rèn)值是30) 如果忽視了perplexity帶來的影響蜈首，有的時(shí)候遇到t-SNE可視化效果不好時(shí)实抡，對(duì)于問題無從下手。perplexity表示了近鄰的數(shù)量欢策，例如設(shè)perplexity為2吆寨，那么就很有可能得到很多兩個(gè)一對(duì)的小集群。
  
• 有的時(shí)候會(huì)出現(xiàn)同一集群被分為兩半的情況踩寇，但群間的距離并不能說明什么啄清，解決這個(gè)問題，只需要跑多次找出效果最好的就可以了

引用自： https://bindog.github.io/blog/2018/07/31/t-sne-tips/
很好的tsne可視化：https://distill.pub/2016/misread-tsne/

和PCA類似俺孙，先runTSNE辣卒，再plotTSNE。另外注意鼠冕，為了重復(fù)結(jié)果要設(shè)置隨機(jī)種子添寺，因?yàn)閠sne每次映射的坐標(biāo)結(jié)果都不同。官方強(qiáng)烈建議懈费，使用不同的隨機(jī)種子和perplexity數(shù)值出圖

# Perplexity of 10 just chosen here arbitrarily. 
set.seed(1000)
example_sce <- runTSNE(example_sce, perplexity=10)
plotTSNE(example_sce, colour_by = "Gene_0001", size_by = "Gene_1000")

還可以使用diffusion maps降維：

example_sce <- runDiffusionMap(example_sce)
plotDiffusionMap(example_sce, colour_by = "Gene_0001", size_by = "Gene_1000")

第四步 質(zhì)控（官網(wǎng) 24 May 2019)

作為關(guān)鍵的一步计露，值得關(guān)注！

質(zhì)控主要包含以下三步：

• 質(zhì)控并過濾細(xì)胞
• 質(zhì)控并過濾feature(基因)
• 對(duì)實(shí)驗(yàn)的一些變量進(jìn)行質(zhì)控

首先登場(chǎng)的是calculateQCMetrics函數(shù)

它對(duì)每個(gè)細(xì)胞和feature信息計(jì)算了大量的統(tǒng)計(jì)指標(biāo)憎乙，分別存在colData和rowData中票罐，這個(gè)函數(shù)默認(rèn)對(duì)count值進(jìn)行計(jì)算，但是也可以通過參數(shù)exprs_values進(jìn)行修改

example_sce <- calculateQCMetrics(example_sce)

> colnames(colData(example_sce)) # 樣本質(zhì)控
 [1] "Cell"                           "Mutation_Status"               
 [3] "Cell_Cycle"                     "Treatment"                     
 [5] "is_cell_control"                "total_features_by_counts"      
 [7] "log10_total_features_by_counts" "total_counts"                  
 [9] "log10_total_counts"             "pct_counts_in_top_50_features" 
[11] "pct_counts_in_top_100_features" "pct_counts_in_top_200_features"
[13] "pct_counts_in_top_500_features"

> colnames(rowData(example_sce)) # feature質(zhì)控
[1] "is_feature_control"    "mean_counts"           "log10_mean_counts"     "n_cells_by_counts"    
[5] "pct_dropout_by_counts" "total_counts"          "log10_total_counts"   

另外可以此時(shí)設(shè)置對(duì)照：對(duì)feature（基因）信息來講泞边，可以添加ERCC该押、線粒體基因等信息，對(duì)細(xì)胞來講阵谚，可以添加empty wells蚕礼、 visually damaged cells等信息烟具，然后接下來計(jì)算的QC指標(biāo)就會(huì)包含這些信息

example_sce <- calculateQCMetrics(example_sce, 
    feature_controls = list(ERCC = 1:20, mito = 500:1000),
    cell_controls = list(empty = 1:5, damaged = 31:40))

all_col_qc <- colnames(colData(example_sce))
all_col_qc <- all_col_qc[grep("ERCC", all_col_qc)]
# 取出的就是ERCC的指標(biāo)
> all_col_qc
[1] "total_features_by_counts_ERCC"       "log10_total_features_by_counts_ERCC"
[3] "total_counts_ERCC"                   "log10_total_counts_ERCC"            
[5] "pct_counts_ERCC"                     "pct_counts_in_top_50_features_ERCC" 
[7] "pct_counts_in_top_100_features_ERCC" "pct_counts_in_top_200_features_ERCC"
[9] "pct_counts_in_top_500_features_ERCC"

如果定義了control組，那么計(jì)算的結(jié)果還會(huì)有一類：就是將control匯總在一起計(jì)算的feature_control和cell_control奠蹬；非control的匯總在一起得到endogenous和non_control（這幾個(gè)名稱在我們命名數(shù)據(jù)集時(shí)應(yīng)該避開）

關(guān)于細(xì)胞質(zhì)控的幾個(gè)重點(diǎn)關(guān)注結(jié)果

• total_counts：每個(gè)細(xì)胞中總表達(dá)量(文庫(kù)大小)

• total_features_by_counts：細(xì)胞中超過規(guī)定閾值(默認(rèn)是0)的feature數(shù)量

• pct_counts_X：屬于feature control組(也就是這里的X)的表達(dá)量占比

  如果使用的不是counts值而是其他矩陣朝聋，比如tpm，那么結(jié)果中也會(huì)將counts替換為tpm

關(guān)于feature質(zhì)控的幾個(gè)重點(diǎn)關(guān)注結(jié)果

• mean_counts：基因/feature的平均表達(dá)量
• pct_dropout_by_counts：基因在細(xì)胞中表達(dá)量為0的細(xì)胞數(shù)占比（該基因丟失率）
• pct_counts_Y：屬于cell control組的表達(dá)量占比
• log10_mean_counts：歸一化 log10 scale
• n_cells_by_counts：多少個(gè)細(xì)胞表達(dá)了該基因

繪圖函數(shù)--檢測(cè)高表達(dá)基因(plotHighestExprs)

默認(rèn)顯示前50個(gè)基因囤躁。圖中每一行表示一個(gè)基因冀痕，每個(gè)線條(bar)表示這個(gè)基因在不同細(xì)胞的表達(dá)量(可以想象成把基因表達(dá)量的箱線圖轉(zhuǎn)了一下)。圓圈是每個(gè)基因表達(dá)量的中位數(shù)狸演，并且在圖中經(jīng)過了排序言蛇。

plotHighestExprs(example_sce, exprs_values = "counts")

我們期待這個(gè)圖會(huì)給出符合常理的結(jié)果，比如：線粒體基因宵距、actin腊尚、ribosomal protein、MALAT1消玄，另外或許有少量的ERCC spike-in在這里有顯示跟伏。ERCC數(shù)量少還可以丢胚，但是如果在top50基因中翩瓜，ERCC占比過多，就表示加入了太高濃度的外源RNA携龟，高ERCC含量與低質(zhì)量數(shù)據(jù)相關(guān)兔跌，通常是排除的標(biāo)準(zhǔn)。除此以外峡蟋，還有可能會(huì)見到一些假基因或預(yù)測(cè)的基因名坟桅，這表明可能比對(duì)過程出現(xiàn)了問題

繪圖函數(shù)—表達(dá)頻率比均值(plotExprsFreqVsMean)

表達(dá)頻率指的是非0表達(dá)量的細(xì)胞數(shù)量(當(dāng)然這個(gè)表達(dá)閾值也可以自己定義)；理論上蕊蝗，對(duì)于大多數(shù)基因來說仅乓，它應(yīng)該和平均表達(dá)量正相關(guān)

plotExprsFreqVsMean(example_sce)

但其實(shí)圖中也能看出一些異常：

• 例如表達(dá)量均值很低，但覆蓋的細(xì)胞數(shù)量由非常高（圖中左上部）蓬戚，可能的原因是：比對(duì)到一些假基因（pseudo-genes ）（假基因是原來的能翻譯成蛋白的基因經(jīng)過各種突變導(dǎo)致喪失功能的基因）夸楣，然后造成了一個(gè)假象

  一般來說，結(jié)尾是P1等字眼的都是假基因（如：KRAS-->KRASP1）子漩；根據(jù)https://www.genenames.org/download/statistics-and-files/數(shù)據(jù)豫喧，目前存在13311個(gè)假基因

• 又如，表達(dá)量均值很高幢泼，但細(xì)胞數(shù)量卻很少（圖中右下部）紧显，可能因?yàn)?strong>PCR擴(kuò)增偏好性（或者存在罕見的細(xì)胞群）

繪圖函數(shù)—總feature表達(dá)量 vs feature control占比(plotColData)

這個(gè)函數(shù)是對(duì)細(xì)胞進(jìn)行質(zhì)控的，原理很簡(jiǎn)單缕棵，就是看看細(xì)胞中實(shí)際表達(dá)量高還是feature control（例如ERCC）的占比高孵班。如果實(shí)際表達(dá)量高涉兽，feature control占比小的話，就說明細(xì)胞質(zhì)量不錯(cuò)篙程；相反花椭，就說明存在blank and failed cells的情況。我們最后就是看散點(diǎn)主要集中在哪邊（右下部效果較好）

plotColData(example_sce, x = "total_features_by_counts",
    y = "pct_counts_feature_control", colour = "Mutation_Status") +
    theme(legend.position = "top") +
    stat_smooth(method = "lm", se = FALSE, size = 2, fullrange = TRUE)

繪圖函數(shù)—表達(dá)量累計(jì)貢獻(xiàn)圖(plotScater)

這個(gè)函數(shù)先在表達(dá)量最高的基因中選一部分（默認(rèn)是500個(gè)）房午，然后從高到低累加矿辽，看看它們對(duì)每個(gè)細(xì)胞文庫(kù)的貢獻(xiàn)值如何。它將不同細(xì)胞的不同基因表達(dá)分布繪制出來郭厌，就像芯片數(shù)據(jù)或bulk轉(zhuǎn)錄組數(shù)據(jù)每個(gè)樣本做一個(gè)箱線圖袋倔，來看基因表達(dá)量分布。但是由于單細(xì)胞數(shù)據(jù)樣本太多折柠，沒辦法全畫出箱線圖宾娜。

為了看不同細(xì)胞的表達(dá)量差異，可以利用colData中的數(shù)據(jù)將細(xì)胞分類扇售；

默認(rèn)使用counts值前塔，但是只要在assays存在的表達(dá)矩陣，就可以在exprs_values參數(shù)中自定義

plotScater(example_sce, block1 = "Mutation_Status", block2 = "Treatment",
     colour_by = "Cell_Cycle", nfeatures = 300, exprs_values = "counts")

這個(gè)圖在處理不同批次細(xì)胞時(shí)就可以很清楚地看到它們之間的差異

繪圖函數(shù)—QC結(jié)果兩兩比較

# 比較feature的屬性
plotRowData(example_sce, x = "n_cells_by_counts", y = "mean_counts")
# 比較樣本的屬性
p1 <- plotColData(example_sce, x = "total_counts", 
    y = "total_features_by_counts")
p2 <- plotColData(example_sce, x = "pct_counts_feature_control",
    y = "total_features_by_counts")
p3 <- plotColData(example_sce, x = "pct_counts_feature_control",
    y = "pct_counts_in_top_50_features")
# 最后可以組合
multiplot(p1, p2, p3, cols = 3)

第五步-1 對(duì)細(xì)胞過濾（官網(wǎng) 24 May 2019)

對(duì)列取子集

列就是細(xì)胞樣本承冰，最簡(jiǎn)單的方法就是直接按照類似數(shù)據(jù)框的操作：

example_sce <- example_sce[,1:40]

另外华弓，還提供了filter函數(shù)（受到dplyr包的同名函數(shù)啟發(fā)）

filter(example_sce, Treatment == "treat1")

設(shè)置條件來過濾

例如設(shè)置細(xì)胞的總表達(dá)量不能低于100,000，總表達(dá)基因數(shù)不能少于500

keep.total <- example_sce$total_counts > 1e5
keep.n <- example_sce$total_features_by_counts > 500
filtered <- example_sce[,keep.total & keep.n]
dim(filtered)
## [1] 2000   37

除此以外困乒，還有一種更靈活的方式：isOutlier函數(shù)寂屏，它設(shè)置的是距離中位數(shù)差幾個(gè)MAD值（絕對(duì)中位差 median absolute deviations），超過設(shè)定值的被認(rèn)為是離群點(diǎn)（outlier）娜搂，它們就會(huì)被舍去迁霎。例如設(shè)置過小的離群點(diǎn)為：log count值低于中位數(shù)的3個(gè)MAD值。

可以看到百宇，這個(gè)函數(shù)沒有“一刀切”考廉，可以根據(jù)測(cè)序深度、spike-in的影響携御、細(xì)胞類型等進(jìn)行調(diào)整昌粤。
如何設(shè)置閾值是一門學(xué)問，官方也推薦看看：https://bioconductor.org/packages/release/workflows/html/simpleSingleCell.html

keep.total <- isOutlier(example_sce$total_counts, nmads=3, 
    type="lower", log=TRUE)
filtered <- example_sce[,keep.total]

根據(jù)QC結(jié)果來過濾

比如用一個(gè)PCA圖：

example_sce <- runPCA(example_sce, use_coldata = TRUE,
    detect_outliers = TRUE)
plotReducedDim(example_sce, use_dimred="PCA_coldata")

結(jié)果就會(huì)在example_sce對(duì)象中增加一個(gè)outlier的接口因痛，我們可以這樣看：

summary(example_sce$outlier)
##    Mode   FALSE 
## logical      40

第五步-2 對(duì)feature過濾（官網(wǎng) 24 May 2019)

最簡(jiǎn)單的方法是用自帶的函數(shù)（當(dāng)然婚苹，也可以在構(gòu)建對(duì)象前就過濾好）：

keep_feature <- nexprs(example_sce, byrow=TRUE) >= 4
example_sce <- example_sce[keep_feature,]
dim(example_sce)
## [1] 1753   40

第六步 探索不同實(shí)驗(yàn)因素對(duì)表達(dá)量的影響

# 一般是要先標(biāo)準(zhǔn)化
example_sce <- normalize(example_sce)
plotExplanatoryVariables(example_sce)

結(jié)果中的每條線都表示一個(gè)影響因子，當(dāng)然我們也可以選擇一部分

plotExplanatoryVariables(example_sce,
    variables = c("total_features_by_counts", "total_counts",
        "Mutation_Status", "Treatment", "Cell_Cycle"))

結(jié)果看到：因?yàn)槭菧y(cè)試小數(shù)據(jù)的緣故鸵膏，total_counts 和 total_features_by_counts 對(duì)基因表達(dá)變化差異貢獻(xiàn)很高（接近10%）膊升，真實(shí)數(shù)據(jù)中這兩個(gè)占比應(yīng)該在1-5%之間

第七步 去除技術(shù)誤差

Scaling normalization

可選的有normalize函數(shù)，官方還強(qiáng)烈推薦了 友包scran的computeSumFactors谭企、computeSpikeFactors函數(shù)

Batch correction

Unlike scaling biases, these are usually constant across all cells in a given batch but different for each gene.

可以用limma的removeBatchEffect函數(shù)

library(limma)
batch <- rep(1:2, each=20)
corrected <- removeBatchEffect(logcounts(example_sce), block=batch)
assay(example_sce, "corrected_logcounts") <- corrected

另外還推薦了：scran的mnnCorrect函數(shù)（之前也介紹過：http://www.reibang.com/p/b7f6a5efed85）

好廓译，上面的官網(wǎng)教程結(jié)束评肆，下面繼續(xù)對(duì)scRNAseq這個(gè)包中的四組細(xì)胞類型進(jìn)行scater操作

回到scRNA數(shù)據(jù)

首先載入數(shù)據(jù)

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm)

# 得到RSEM矩陣
assay(fluidigm)  <- assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)

> ct[1:4,1:4]
         SRR1275356 SRR1274090 SRR1275251 SRR1275287
A1BG              0          0          0          0
A1BG-AS1          0          0          0          0
A1CF              0          0          0          0
A2M               0          0          0         33
> table(rowSums(ct)==0)
FALSE  TRUE 
17096  9159 

# 樣本注釋信息
pheno_data <- as.data.frame(colData(fluidigm))
# 然后做scater的數(shù)據(jù)
sce <- SingleCellExperiment(
    assays = list(counts = ct), 
    colData = pheno_data
    )
> sce
class: SingleCellExperiment 
dim: 26255 130 
metadata(0):
assays(1): counts
rownames(26255): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
rowData names(0):
colnames(130): SRR1275356 SRR1274090 ... SRR1275366 SRR1275261
colData names(28): NREADS NALIGNED ... Cluster1 Cluster2
reducedDimNames(0):
spikeNames(0):

簡(jiǎn)單質(zhì)控

sce <- calculateQCMetrics(sce)

然后是過濾

基因?qū)用?/h5>

# 進(jìn)行一個(gè)標(biāo)準(zhǔn)化
exprs(sce) <- log2(calculateCPM(sce) + 1)

genes <- rownames(rowData(sce))
genes[grepl('^MT-',genes)] # 檢測(cè)線粒體基因
genes[grepl('^ERCC-',genes)] # 檢測(cè)ERCC
# 假入有線粒體基因或ERCC，就可以加入feature control組中
sce <- calculateQCMetrics(sce, 
        feature_controls = list(ERCC = grep('^ERCC',genes)))

colnames(rowData(sce)) # 查看信息
# 只過濾那些在所有細(xì)胞都沒有表達(dá)的基因
keep_feature <- rowSums(exprs(sce) > 0) > 0
# table(keep_feature)
sce <- sce[keep_feature,]
# sce

細(xì)胞層面

# 細(xì)胞質(zhì)控屬性非常多
colnames(colData(sce)) # 查看信息

可視化同上

學(xué)習(xí)SC3包

library(SC3) # BiocManager::install('SC3')
sce <- sc3_estimate_k(sce) # 預(yù)估亞群(24個(gè))
metadata(sce)$sc3$k_estimation
rowData(sce)$feature_symbol <- rownames(rowData(sce))

# 接下來正式運(yùn)行非区，kn參數(shù)表示的預(yù)估聚類數(shù)
# 我們這里自定義為4組(因?yàn)橐阎嬲褪?組瓜挽，實(shí)際上需要探索)
kn <- 4 # 還可以設(shè)置3、5看看結(jié)果
sc3_cluster <- "sc3_4_clusters"
Sys.time()
sce <- sc3(sce, ks = kn, biology = TRUE) # 運(yùn)行會(huì)很慢
Sys.time()

# 最后進(jìn)行可視化--比較先驗(yàn)分類和SC3的聚類的一致性
sc3_plot_consensus(sce, k = kn, show_pdata = c("Biological_Condition",sc3_cluster))
# 最后進(jìn)行可視化--表達(dá)量信息
sc3_plot_expression(sce, k = kn, show_pdata =  c("Biological_Condition",sc3_cluster))
# 最后進(jìn)行可視化--可能的marker基因
sc3_plot_markers(sce, k = kn, show_pdata =  c("Biological_Condition",sc3_cluster))
# PCA上展示SC3的聚類結(jié)果
plotPCA(sce, colour_by =  sc3_cluster )