1.下載和整理數(shù)據(jù)
示例https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE231920
需要3個(gè)文件：
GSM7306054_sample1_barcodes.tsv.gz
GSM7306054_sample1_features.tsv.gz
GSM7306054_sample1_matrix.mtx.gz
下下來(lái)是個(gè)壓縮包
1.1 解包文件
untar("GSE231920_RAW.tar",exdir = "input")
1.2 單細(xì)胞文件組織的要求
這三個(gè)文件是10X的標(biāo)準(zhǔn)文件，需要放在同一個(gè)文件夾里，并且名字是固定的，不能有前綴。
“barcodes.tsv.gz”：存儲(chǔ)的是barcodes，相當(dāng)于細(xì)胞的編號(hào)，是表達(dá)矩陣的列名。
“features.tsv.gz”：存儲(chǔ)的是基因名稱，是表達(dá)矩陣的行名，也可以是”genes.tsv.gz”。
“matrix.mtx.gz”:存儲(chǔ)的是每個(gè)位置的數(shù)值，是表達(dá)矩陣的內(nèi)容，僅存儲(chǔ)了非零的數(shù)值。
1.3 修改文件名稱
library(stringr)
nn = str_remove(dir("input/"),"GSM7306054_sample1_")
nn
## [1] "barcodes.tsv.gz" "features.tsv.gz" "matrix.mtx.gz"
file.rename(paste0("input/",dir("input/")),
paste0("input/",nn))
## [1] TRUE TRUE TRUE
dir("input/")
## [1] "barcodes.tsv.gz" "features.tsv.gz" "matrix.mtx.gz"
2.讀取并創(chuàng)建Seurat對(duì)象
2.1讀取文件
rm(list = ls())
#清空右上角的所有變量，方便反復(fù)調(diào)試代碼
library(Seurat)
library(patchwork)
library(tidyverse)
ct = Read10X("input/")
#讀取標(biāo)準(zhǔn)文件，接受的參數(shù)是文件夾名稱。文件夾里的三個(gè)文件合到一起才是完整的單細(xì)胞表達(dá)矩陣。
dim(ct)
## [1] 33538 10218
兩個(gè)數(shù)字分別是行數(shù)(基因數(shù))和列數(shù)(細(xì)胞數(shù))。
2.2 稀疏矩陣
class(ct)
## [1] "dgCMatrix"
## attr(,"package")
## [1] "Matrix"
ct[c("CD3D","TCL1A","MS4A1"), 1:30]
## 3 x 30 sparse Matrix of class "dgCMatrix"
##
## CD3D . . 5 . . 1 . . 2 . 3 . . . 4 . . 3 . . . . . . . 5 3 . . .
## TCL1A . . . . . . . . . . . . 3 . . . . . . . . 1 1 . . . . . . .
## MS4A1 4 . . . . . . . . 4 . . 1 . 1 . 2 . . . . 2 4 7 5 . . . 1 .
稀疏矩陣是存儲(chǔ)0值比較多的數(shù)據(jù)用的，用“.”表示0，可以節(jié)省空間。單細(xì)胞矩陣?yán)锩嬗写罅康?值。
2.3 創(chuàng)建Seurat對(duì)象
seu.obj <- CreateSeuratObject(counts = ct,
min.cells = 3, #一個(gè)基因至少要在3個(gè)細(xì)胞里有表達(dá)，才被保留
min.features = 200) #一個(gè)細(xì)胞里至少要200個(gè)基因有表達(dá)，才被保留
dim(seu.obj)
## [1] 20648 10105
2.4 細(xì)胞抽樣
!!!這一步是為了為了節(jié)省計(jì)算資源，我們減一下細(xì)胞的數(shù)量，實(shí)戰(zhàn)時(shí)不能減!!!!
set.seed(1234) #set.seed(1234)是設(shè)計(jì)隨機(jī)種子，作用是讓后面的隨機(jī)抽樣變得可重復(fù)，只要多次運(yùn)行時(shí)，setseed的括號(hào)里的數(shù)字沒(méi)變，抽樣的結(jié)果就不會(huì)變。
seu.obj = subset(seu.obj,downsample = 3000)
2.5 對(duì)象
廣義的“對(duì)象”是Rstudio右上角所有的數(shù)據(jù)，向量數(shù)據(jù)框矩陣列表都是對(duì)象。
狹義的“對(duì)象”是由R包的作者定義的，以固定模式組織的數(shù)據(jù)，里面的結(jié)構(gòu)都是規(guī)定好的。
可以用class查看
class(seu.obj)
## [1] "Seurat"
## attr(,"package")
## [1] "SeuratObject"
#意思是：這是一個(gè)Seurat對(duì)象，出自于SeuratObject這個(gè)包。
提取它的子集，有@和$兩個(gè)符號(hào),具體該用哪一個(gè)可以試試，一般第一層次都是@，后面的就不一定了
2.6 探索Seurat對(duì)象的meta信息
seu.obj@meta.data %>% head()
## orig.ident nCount_RNA nFeature_RNA
## AAACCCACAGTCGGTC-1 SeuratProject 4243 1256
## AAACGAAAGAATCGCG-1 SeuratProject 7307 2577
## AAAGAACAGCTTACGT-1 SeuratProject 8154 1881
## AAAGAACAGGTTCATC-1 SeuratProject 8223 2182
## AAAGAACAGTCTGTAC-1 SeuratProject 3884 1377
## AAAGAACTCCACCTCA-1 SeuratProject 3997 1307
或
seu.obj[[]] %>% head()
## orig.ident nCount_RNA nFeature_RNA
## AAACCCACAGTCGGTC-1 SeuratProject 4243 1256
## AAACGAAAGAATCGCG-1 SeuratProject 7307 2577
## AAAGAACAGCTTACGT-1 SeuratProject 8154 1881
## AAAGAACAGGTTCATC-1 SeuratProject 8223 2182
## AAAGAACAGTCTGTAC-1 SeuratProject 3884 1377
## AAAGAACTCCACCTCA-1 SeuratProject 3997 1307
行名是barcodes,也就是表達(dá)矩陣?yán)锩娴牧忻?，相?dāng)于細(xì)胞的編號(hào)。
orig.ident是細(xì)胞的原始分類，如果是單樣本數(shù)據(jù)，在前面CreateSeuratObject時(shí)，如果指定project參數(shù)，就會(huì)顯示在這一列，整列都是同一個(gè)單詞，因?yàn)槲覀兦懊鏇](méi)指定，所以顯示默認(rèn)值SeuratProject。
nCount_RNA是每個(gè)細(xì)胞中所有基因的表達(dá)量之和，也就是表達(dá)矩陣的列和。
nFeature_RNA是每個(gè)細(xì)胞中表達(dá)量不為零的基因的數(shù)量。
exp = seu.obj@assays$RNA$counts #這是表達(dá)矩陣
exp[1:4,1:4]
## 4 x 4 sparse Matrix of class "dgCMatrix"
## AAACCCACAGTCGGTC-1 AAACGAAAGAATCGCG-1 AAAGAACAGCTTACGT-1
## AL627309.1 . . .
## AL627309.3 . . .
## AL669831.5 . . .
## FAM87B . . .
## AAAGAACAGGTTCATC-1
## AL627309.1 .
## AL627309.3 .
## AL669831.5 .
## FAM87B .
sum(exp[,1])
## [1] 4243
table(exp[,1]!=0) #統(tǒng)計(jì)TRUE和FALSE的個(gè)數(shù)
##
## FALSE TRUE
## 19392 1256
3.質(zhì)控
3.1 質(zhì)控的指標(biāo)及原因
這里是對(duì)細(xì)胞進(jìn)行的質(zhì)控，指標(biāo)是：（1）線粒體基因含量不能過(guò)高；（2）nFeature_RNA 不能過(guò)高或過(guò)低
nFeature_RNA是每個(gè)細(xì)胞中檢測(cè)到的基因數(shù)量。nCount_RNA是細(xì)胞內(nèi)檢測(cè)到的分子總數(shù)。nFeature_RNA過(guò)低，表示該細(xì)胞可能已死/將死或是空液滴。太高的nCount_RNA和/或nFeature_RNA表明“細(xì)胞”實(shí)際上可以是兩個(gè)或多個(gè)細(xì)胞。結(jié)合線粒體基因count數(shù)除去異常值，即可除去大多數(shù)雙峰/死細(xì)胞/空液滴，因此它們過(guò)濾是常見(jiàn)的預(yù)處理步驟。 參考自：https://www.biostars.org/p/407036/
3.2計(jì)算線粒體基因比例
人類的線粒體基因都是以MT-開(kāi)頭
seu.obj[["percent.mt"]] <- PercentageFeatureSet(seu.obj, pattern = "^MT-")
head(seu.obj@meta.data)
## orig.ident nCount_RNA nFeature_RNA percent.mt
## AAACCCACAGTCGGTC-1 SeuratProject 4243 1256 6.292717
## AAACGAAAGAATCGCG-1 SeuratProject 7307 2577 2.572875
## AAAGAACAGCTTACGT-1 SeuratProject 8154 1881 4.083885
## AAAGAACAGGTTCATC-1 SeuratProject 8223 2182 4.377964
## AAAGAACAGTCTGTAC-1 SeuratProject 3884 1377 4.763131
## AAAGAACTCCACCTCA-1 SeuratProject 3997 1307 3.402552
VlnPlot(seu.obj,
features = c("nFeature_RNA",
"nCount_RNA",
"percent.mt"),
ncol = 3,pt.size = 0.5)
小提琴的寬窄代表對(duì)應(yīng)縱坐標(biāo)的細(xì)胞數(shù)量多少。如果以后遇到點(diǎn)的數(shù)量太多，把小提琴都給蓋上的情況，那就pt.size = 0，大小為0就是不畫(huà)點(diǎn)。
根據(jù)上面的小提琴圖里可以看到所有細(xì)胞的三個(gè)指標(biāo)的分布情況，我們過(guò)濾就是把三個(gè)指標(biāo)比大部隊(duì)數(shù)值大的點(diǎn)給過(guò)濾掉。例如這個(gè)數(shù)據(jù)過(guò)濾標(biāo)準(zhǔn)可以是：
seu.obj = subset(seu.obj,
nFeature_RNA < 6000 &
nCount_RNA < 30000 &
percent.mt < 18)
標(biāo)準(zhǔn)不是唯一的，大點(diǎn)兒小點(diǎn)兒?jiǎn)栴}不大。如果沒(méi)有尖尖的線，不過(guò)濾也可以。有的公共數(shù)據(jù)拿到時(shí)就已經(jīng)是過(guò)濾好的。尖尖部分對(duì)應(yīng)的縱坐標(biāo)數(shù)值太大，屬于離群值，就不要了。 一般過(guò)濾掉尖尖部分就可以了，不要過(guò)濾的太狠，標(biāo)準(zhǔn)見(jiàn)下：
過(guò)濾后
4.降維聚類分群
4.1 理解降維這件事
對(duì)于表達(dá)矩陣來(lái)說(shuō)，一個(gè)基因就是一個(gè)維度，有幾萬(wàn)個(gè)維度，太多了，需要減少，稱之為降維。
PCA是主成分分析，完成線性的降維，會(huì)把幾萬(wàn)個(gè)維度轉(zhuǎn)換為幾十個(gè)新的維度，而UMAP和tSNE是非線性的降維，線性降維不夠徹底，整點(diǎn)高級(jí)的，進(jìn)一步降維。
UMAP 和tSNE二選一就行，沒(méi)有必須用哪個(gè)的說(shuō)法，只能說(shuō)UMAP用的多一些，圖一般會(huì)緊湊一些。
4.2 file.exists——存在即跳過(guò)
file.exists("filename")
#存在的文件名會(huì)返回TRUE，而不存在與工作目錄下的文件名則會(huì)返回FALSE
因此下面的代碼是，結(jié)合if語(yǔ)句實(shí)現(xiàn)了分情況討論：
如果工作目錄下不存在”obj.Rdata”這個(gè)文件，則運(yùn)行大括號(hào)里的代碼(最后一句是save，所以運(yùn)行完也就保存了”obj.Rdata”)； 如果工作目錄下存在”obj.Rdata”這個(gè)文件(說(shuō)明大括號(hào)里的代碼已經(jīng)運(yùn)行過(guò)了)，則不運(yùn)行大括號(hào)里的代碼，也就跳過(guò)了這個(gè)耗時(shí)比較長(zhǎng)的限速步驟
總之這個(gè)技巧可以避免多次重復(fù)運(yùn)行限速步驟，非常實(shí)用。但需要注意?，一旦輸入數(shù)據(jù)有改動(dòng)，比如說(shuō)前面的過(guò)濾標(biāo)準(zhǔn)改了，那么就要?jiǎng)h除工作目錄下的obj.Rdata文件，才能重新運(yùn)行這段代碼。
f = "obj.Rdata"
if(!file.exists(f)){
seu.obj = seu.obj %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData(features = rownames(.)) %>%
RunPCA(features = VariableFeatures(.)) %>%
FindNeighbors(dims = 1:15) %>% #選擇前15個(gè)維度，15是個(gè)比較折中的值，選擇的數(shù)量越多計(jì)算量越大，而且信息越冗余；選的太少又不具有代表性。一般是10-20
FindClusters(resolution = 0.5) %>% #分群的分辨率是0.5，這個(gè)參數(shù)的選擇范圍是0.1-1.5之間，數(shù)值越大分出來(lái)的群越多，數(shù)值越小分出來(lái)的點(diǎn)越少
RunUMAP(dims = 1:15) %>%
RunTSNE(dims = 1:15)
save(seu.obj,file = f)
}
load(f)
ElbowPlot(seu.obj)
p1 <- DimPlot(seu.obj, reduction = "umap",label = T)+NoLegend();p1
標(biāo)記#的兩個(gè)地方可能會(huì)需要改動(dòng)
一個(gè)是dims = 1:15,代表選擇前15個(gè)維度，15是個(gè)比較折中的值，選擇的數(shù)量越多計(jì)算量越大，而且信息越冗余；選的太少又不具有代表性。一般是10-20，根據(jù)ElbowPlot來(lái)選擇拐點(diǎn)的橫坐標(biāo)(拐點(diǎn)就是在這個(gè)點(diǎn)之前縱坐標(biāo)下降比較快，在這個(gè)點(diǎn)之后縱坐標(biāo)下降比較慢)。不是很重要，只要拐點(diǎn)在15之前，選15就一點(diǎn)問(wèn)題都沒(méi)有。
一個(gè)是resolution = 0.5，代表分群的分辨率是0.5，這個(gè)參數(shù)的選擇范圍是0.1-1.5之間，數(shù)值越大分出來(lái)的群越多，數(shù)值越小分出來(lái)的點(diǎn)越少。0.5也是一個(gè)比較折中的值，如果后面注釋不困難，就不用改動(dòng)。
如果決定要改動(dòng)那么就不能只改代碼，要把obj.Rdata從工作目錄下刪掉，再重新運(yùn)行修改后的代碼。
