【在公粽號回復(fù)20210105就能看到示例數(shù)據(jù)了】
前面我們已經(jīng)學(xué)習(xí)了單細(xì)胞轉(zhuǎn)錄組分析的：使用Cell Ranger得到表達(dá)矩陣和doublet檢測张峰，今天我們開始Seurat標(biāo)準(zhǔn)流程的學(xué)習(xí)锐墙。這一部分的內(nèi)容，網(wǎng)上有很多帖子，基本上都是把Seurat官網(wǎng)PBMC的例子重復(fù)一遍，這回我換一個數(shù)據(jù)集，細(xì)胞類型更多，同時也會加入一些實際分析中很有用的技巧轩娶。
1. 導(dǎo)入數(shù)據(jù)儿奶，創(chuàng)建Seurat對象
library(Seurat)
library(tidyverse)
testdf=read.table("test_20210105.txt",header = T,row.names = 1)
test.seu=CreateSeuratObject(counts = testdf)
看一下長什么樣子
> test.seu
An object of class Seurat
33538 features across 6746 samples within 1 assay
Active assay: RNA (33538 features)
#①1 assay表示有一套數(shù)據(jù)框往，假如用Seurat里面的函數(shù)去過批次效應(yīng)，
#這里會有2個assay闯捎，另外一個是去批次（整合）之后的椰弊；
#②cell hashing的tag表達(dá)矩陣生成Seurat對象，這時的assay為"HTO"瓤鼻，不叫"RNA"秉版；
#其他情況類似
測試數(shù)據(jù)有33538個基因，6746個細(xì)胞茬祷。除此之外清焕，還要關(guān)注一下另外兩層信息：test.seu@meta.data這個數(shù)據(jù)框用來存儲元數(shù)據(jù)，每一個細(xì)胞都有多個屬性祭犯；test.seu[["RNA"]]@counts這個稀疏矩陣用來存儲原始UMI表達(dá)矩陣秸妥。
> head(test.seu@meta.data)
orig.ident nCount_RNA nFeature_RNA
A_AAACCCAAGGGTCACA A 3714 1151
A_AAACCCAAGTATAACG A 1855 816
A_AAACCCAGTCTCTCAC A 1530 823
A_AAACCCAGTGAGTCAG A 11145 1087
A_AAACCCAGTGGCACTC A 2289 834
A_AAACGAAAGCCAGAGT A 3714 990
#我這里CB的前面人為加上了樣本來源，用下劃線連接沃粗，orig.ident是自動識別得到的
> test.seu[["RNA"]]@counts[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
A_AAACCCAAGGGTCACA A_AAACCCAAGTATAACG A_AAACCCAGTCTCTCAC A_AAACCCAGTGAGTCAG
MIR1302-2HG . . . .
FAM138A . . . .
OR4F5 . . . .
AL627309.1 . . . .
2. 簡單過濾
接下來粥惧，我們根據(jù)每個細(xì)胞內(nèi)部線粒體基因表達(dá)占比、檢測到的基因數(shù)最盅、檢測的UMI總數(shù)這三個方面來對細(xì)胞進(jìn)行簡單的過濾突雪。
先計算細(xì)胞內(nèi)線粒體基因表達(dá)占比，類似的核糖體基因(大多為RP開頭)也能這樣計算涡贱，還要注意不要將線粒體基因的MT-寫成了MT咏删，不然就把別的基因也算進(jìn)去了：
test.seu[["percent.mt"]] <- PercentageFeatureSet(test.seu, pattern = "^MT-") #正則表達(dá)式，表示以MT-開頭问词；test.seu[["percent.mt"]]這種寫法會在meta.data矩陣加上一列
這里我已經(jīng)根據(jù)預(yù)先設(shè)定好的閾值過濾了督函，代碼如下
test.seu <- subset(test.seu, subset = nCount_RNA > 1000 &
nFeature_RNA < 5000 &
percent.mt < 30 &
nFeature_RNA > 600)
過濾之后的數(shù)值分布如下，用到VlnPlot()函數(shù)戏售，該繪圖函數(shù)里面的feature參數(shù)可以是meta.data矩陣的某一列侨核，也可以是某一個基因，很多文章都用這種圖展示marker gene
VlnPlot(test.seu,features = c("nCount_RNA", "nFeature_RNA", "percent.mt"),pt.size = 0)
nFeature_RNA/nCount_RNA不能太泄嘣帧（空液滴）搓译，不能太大（doublet、測序技術(shù)限制）锋喜， 而且閾值設(shè)定要綜合多個樣本來看些己，像下面這樣
一般在CD45陰性的細(xì)胞中percent.mt的閾值大一些豌鸡，50%也看過幾次了
3. LogNormalize，消除文庫大小的影響
如何標(biāo)準(zhǔn)化：LogNormalize: Feature counts for each cell are divided by the total counts for that cell and multiplied by the scale.factor. This is then natural-log transformed using log1p.（先相除段标，再求對數(shù)）
test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)
標(biāo)準(zhǔn)化之后的矩陣存儲在test.seu[["RNA"]]@data
> test.seu[["RNA"]]@data[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
A_AAACCCAAGGGTCACA A_AAACCCAAGTATAACG A_AAACCCAGTCTCTCAC A_AAACCCAGTGAGTCAG
MIR1302-2HG . . . .
FAM138A . . . .
OR4F5 . . . .
AL627309.1 . . . .
4. 找Variable基因
因為單細(xì)胞表達(dá)矩陣很稀疏（很多0）涯冠，選high variable基因的目的可以找到包含信息最多的基因（很多基因的表達(dá)差不多都是0），同時極大提升軟件運(yùn)行速度
test.seu <- FindVariableFeatures(test.seu, selection.method = "vst", nfeatures = 2000)
這些基因存儲在VariableFeatures(test.seu)逼庞，有時候可能需要人為指定high variable基因蛇更，可以這樣：
VariableFeatures(test.seu)="specific genes"
5. scale表達(dá)矩陣
（基于前面得到的data矩陣）
這一步之后，所有基因的表達(dá)值的分布就差不多了赛糟，不然表達(dá)值不在一個數(shù)量級派任，對后續(xù)降維聚類影響挺大。新的矩陣存儲在test.seu[["RNA"]]@scale.data里面璧南。
test.seu <- ScaleData(test.seu, features = rownames(test.seu))
默認(rèn)只對上一步選出來的基因scale掌逛，這里調(diào)整為所有基因，是為了方便以后畫熱圖（畫熱圖一般會用scale之后的z-score）
6. 降維聚類
（基于前面得到的high variable基因的scale矩陣）
test.seu <- RunPCA(test.seu, npcs = 50, verbose = FALSE)
test.seu <- FindNeighbors(test.seu, dims = 1:30)
test.seu <- FindClusters(test.seu, resolution = 0.5)
test.seu <- RunUMAP(test.seu, dims = 1:30)
test.seu <- RunTSNE(test.seu, dims = 1:30)
Run開頭的函數(shù)降維司倚，F(xiàn)ind開頭的函數(shù)聚類豆混，一般就這幾步，相對固定动知。PCA將原來2000維的數(shù)據(jù)降到50維皿伺，dims參數(shù)表示使用多少個主成分（一般20左右就可以了，多幾個少幾個對結(jié)果影響不大）拍柒，resolution參數(shù)表達(dá)聚類的分辨率心傀，這個值大于0，一般都是在0-1范圍里面調(diào)整拆讯，越大得到的cluster越多脂男，這個值可以反復(fù)調(diào)整，并不會改變降維的結(jié)果（也就是tsne种呐、umap圖的二維坐標(biāo)）宰翅。
這一步之后的數(shù)據(jù)是這樣的
> test.seu
An object of class Seurat
33538 features across 6746 samples within 1 assay
Active assay: RNA (33538 features)
3 dimensional reductions calculated: pca, umap, tsne
# 幾種降維方式都會呈現(xiàn)出來
聚類之后test.seu@meta.data多了兩列，RNA_snn_res.0.5記錄了你用的分辨率爽室，最終的聚類結(jié)果保存在seurat_clusters中
> head(test.seu@meta.data)
orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5
A_AAACCCAAGGGTCACA A 3714 1151 9.585353 8
A_AAACCCAAGTATAACG A 1855 816 12.776280 0
A_AAACCCAGTCTCTCAC A 1530 823 14.248366 12
A_AAACCCAGTGAGTCAG A 11145 1087 2.853297 4
A_AAACCCAGTGGCACTC A 2289 834 15.640017 3
A_AAACGAAAGCCAGAGT A 3714 990 5.654281 0
seurat_clusters
A_AAACCCAAGGGTCACA 8
A_AAACCCAAGTATAACG 0
A_AAACCCAGTCTCTCAC 12
A_AAACCCAGTGAGTCAG 4
A_AAACCCAGTGGCACTC 3
A_AAACGAAAGCCAGAGT 0
7. tsne/umap展示結(jié)果
library(cowplot)
test.seu$patient=str_replace(test.seu$orig.ident,"_.*$","")
p1 <- DimPlot(test.seu, reduction = "tsne", group.by = "patient", pt.size=0.5)
p2 <- DimPlot(test.seu, reduction = "tsne", group.by = "ident", pt.size=0.5, label = TRUE,repel = TRUE) #后面兩個參數(shù)用來添加文本標(biāo)簽
p3 <- DimPlot(test.seu, reduction = "umap", group.by = "patient", pt.size=0.5)
p4 <- DimPlot(test.seu, reduction = "umap", group.by = "ident", pt.size=0.5, label = TRUE,repel = TRUE)
fig_tsne <- plot_grid(p1, p2, labels = c('patient','ident'),align = "v",ncol = 2)
ggsave(filename = "tsne.pdf", plot = fig_tsne, device = 'pdf', width = 30, height = 15, units = 'cm')
fig_umap <- plot_grid(p3, p4, labels = c('patient','ident'),align = "v",ncol = 2)
ggsave(filename = "umap.pdf", plot = fig_umap, device = 'pdf', width = 30, height = 15, units = 'cm')
ident表示每個細(xì)胞的標(biāo)簽汁讼，聚類之后就是聚類的結(jié)果，在一些特定場景可以更換阔墩。
在umap圖中嘿架，cluster之間的距離更明顯
從上面的圖可以看出不同樣本其實是有批次效應(yīng)的，下一講我會介紹兩種去批次效應(yīng)的方法啸箫。
因水平有限耸彪，有錯誤的地方，歡迎批評指正忘苛！
