1.下載和整理數(shù)據(jù)

示例https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE231920
需要3個文件：
GSM7306054_sample1_barcodes.tsv.gz
GSM7306054_sample1_features.tsv.gz
GSM7306054_sample1_matrix.mtx.gz
下下來是個壓縮包

1.1 解包文件

untar("GSE231920_RAW.tar",exdir = "input")

1.2 單細胞文件組織的要求

這三個文件是10X的標(biāo)準(zhǔn)文件，需要放在同一個文件夾里，并且名字是固定的拳亿，不能有前綴。
“barcodes.tsv.gz”：存儲的是barcodes谈截，相當(dāng)于細胞的編號，是表達矩陣的列名。
“features.tsv.gz”：存儲的是基因名稱适贸，是表達矩陣的行名黎侈，也可以是”genes.tsv.gz”察署。
“matrix.mtx.gz”:存儲的是每個位置的數(shù)值，是表達矩陣的內(nèi)容峻汉，僅存儲了非零的數(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對象

2.1讀取文件

rm(list = ls()) 
#清空右上角的所有變量，方便反復(fù)調(diào)試代碼
library(Seurat)
library(patchwork)
library(tidyverse)
ct = Read10X("input/")
#讀取標(biāo)準(zhǔn)文件休吠，接受的參數(shù)是文件夾名稱扳埂。文件夾里的三個文件合到一起才是完整的單細胞表達矩陣。
dim(ct)
## [1] 33538 10218

兩個數(shù)字分別是行數(shù)(基因數(shù))和列數(shù)(細胞數(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 .

稀疏矩陣是存儲0值比較多的數(shù)據(jù)用的阳懂，用“.”表示0，可以節(jié)省空間。單細胞矩陣?yán)锩嬗写罅康?值希太。

2.3 創(chuàng)建Seurat對象

seu.obj <- CreateSeuratObject(counts = ct, 
                           min.cells = 3,         #一個基因至少要在3個細胞里有表達克饶，才被保留
                           min.features = 200)    #一個細胞里至少要200個基因有表達，才被保留
dim(seu.obj)
## [1] 20648 10105

2.4 細胞抽樣

!!!這一步是為了為了節(jié)省計算資源誊辉，我們減一下細胞的數(shù)量矾湃，實戰(zhàn)時不能減!!!!

set.seed(1234)  #set.seed(1234)是設(shè)計隨機種子，作用是讓后面的隨機抽樣變得可重復(fù)堕澄，只要多次運行時邀跃，setseed的括號里的數(shù)字沒變，抽樣的結(jié)果就不會變蛙紫。
seu.obj = subset(seu.obj,downsample = 3000)

2.5 對象

廣義的“對象”是Rstudio右上角所有的數(shù)據(jù)拍屑，向量數(shù)據(jù)框矩陣列表都是對象。
狹義的“對象”是由R包的作者定義的坑傅，以固定模式組織的數(shù)據(jù)僵驰，里面的結(jié)構(gòu)都是規(guī)定好的。
可以用class查看

class(seu.obj)
## [1] "Seurat"
## attr(,"package")
## [1] "SeuratObject"

#意思是：這是一個Seurat對象唁毒，出自于SeuratObject這個包蒜茴。

提取它的子集，有@和$兩個符號,具體該用哪一個可以試試浆西，一般第一層次都是@粉私，后面的就不一定了

2.6 探索Seurat對象的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,也就是表達矩陣?yán)锩娴牧忻喈?dāng)于細胞的編號近零。
orig.ident是細胞的原始分類诺核，如果是單樣本數(shù)據(jù)，在前面CreateSeuratObject時久信，如果指定project參數(shù)窖杀，就會顯示在這一列，整列都是同一個單詞裙士，因為我們前面沒指定入客，所以顯示默認(rèn)值SeuratProject。
nCount_RNA是每個細胞中所有基因的表達量之和潮售，也就是表達矩陣的列和。
nFeature_RNA是每個細胞中表達量不為零的基因的數(shù)量锅风。

exp = seu.obj@assays$RNA$counts #這是表達矩陣
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)計TRUE和FALSE的個數(shù)
## 
## FALSE  TRUE 
## 19392  1256

3.質(zhì)控

3.1 質(zhì)控的指標(biāo)及原因

這里是對細胞進行的質(zhì)控酥诽，指標(biāo)是：（1）線粒體基因含量不能過高；（2）nFeature_RNA 不能過高或過低

nFeature_RNA是每個細胞中檢測到的基因數(shù)量皱埠。nCount_RNA是細胞內(nèi)檢測到的分子總數(shù)肮帐。nFeature_RNA過低，表示該細胞可能已死/將死或是空液滴。太高的nCount_RNA和/或nFeature_RNA表明“細胞”實際上可以是兩個或多個細胞训枢。結(jié)合線粒體基因count數(shù)除去異常值托修，即可除去大多數(shù)雙峰/死細胞/空液滴，因此它們過濾是常見的預(yù)處理步驟恒界。 參考自：https://www.biostars.org/p/407036/

3.2計算線粒體基因比例

人類的線粒體基因都是以MT-開頭

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)

plot_zoom.png

小提琴的寬窄代表對應(yīng)縱坐標(biāo)的細胞數(shù)量多少睦刃。如果以后遇到點的數(shù)量太多，把小提琴都給蓋上的情況十酣，那就pt.size = 0涩拙，大小為0就是不畫點。
根據(jù)上面的小提琴圖里可以看到所有細胞的三個指標(biāo)的分布情況耸采，我們過濾就是把三個指標(biāo)比大部隊數(shù)值大的點給過濾掉兴泥。例如這個數(shù)據(jù)過濾標(biāo)準(zhǔn)可以是：

seu.obj = subset(seu.obj,
                nFeature_RNA < 6000 &
                nCount_RNA < 30000 &
                percent.mt < 18)

標(biāo)準(zhǔn)不是唯一的，大點兒小點兒問題不大虾宇。如果沒有尖尖的線搓彻，不過濾也可以。有的公共數(shù)據(jù)拿到時就已經(jīng)是過濾好的嘱朽。尖尖部分對應(yīng)的縱坐標(biāo)數(shù)值太大旭贬，屬于離群值，就不要了燥翅。 一般過濾掉尖尖部分就可以了骑篙，不要過濾的太狠，標(biāo)準(zhǔn)見下：

2024-06-18-151153.png

過濾后

plot_zoom2.png

4.降維聚類分群

4.1 理解降維這件事

對于表達矩陣來說森书，一個基因就是一個維度靶端，有幾萬個維度，太多了凛膏，需要減少杨名，稱之為降維。
PCA是主成分分析猖毫，完成線性的降維台谍，會把幾萬個維度轉(zhuǎn)換為幾十個新的維度，而UMAP和tSNE是非線性的降維吁断，線性降維不夠徹底趁蕊，整點高級的，進一步降維仔役。
UMAP 和tSNE二選一就行掷伙，沒有必須用哪個的說法，只能說UMAP用的多一些又兵，圖一般會緊湊一些任柜。

4.2 file.exists——存在即跳過

file.exists("filename")
#存在的文件名會返回TRUE卒废，而不存在與工作目錄下的文件名則會返回FALSE

因此下面的代碼是，結(jié)合if語句實現(xiàn)了分情況討論：

如果工作目錄下不存在”obj.Rdata”這個文件宙地，則運行大括號里的代碼(最后一句是save摔认，所以運行完也就保存了”obj.Rdata”)； 如果工作目錄下存在”obj.Rdata”這個文件(說明大括號里的代碼已經(jīng)運行過了)宅粥，則不運行大括號里的代碼参袱，也就跳過了這個耗時比較長的限速步驟
總之這個技巧可以避免多次重復(fù)運行限速步驟，非常實用粹胯。但需要注意?蓖柔，一旦輸入數(shù)據(jù)有改動，比如說前面的過濾標(biāo)準(zhǔn)改了风纠，那么就要刪除工作目錄下的obj.Rdata文件况鸣，才能重新運行這段代碼。

f = "obj.Rdata"
if(!file.exists(f)){
  seu.obj = seu.obj %>% 
  NormalizeData() %>%  
  FindVariableFeatures() %>%  
  ScaleData(features = rownames(.)) %>%  
  RunPCA(features = VariableFeatures(.))  %>%
  FindNeighbors(dims = 1:15) %>%   #選擇前15個維度竹观，15是個比較折中的值镐捧，選擇的數(shù)量越多計算量越大，而且信息越冗余臭增；選的太少又不具有代表性懂酱。一般是10-20
  FindClusters(resolution = 0.5) %>%   #分群的分辨率是0.5，這個參數(shù)的選擇范圍是0.1-1.5之間誊抛，數(shù)值越大分出來的群越多列牺，數(shù)值越小分出來的點越少
  RunUMAP(dims = 1:15) %>% 
  RunTSNE(dims = 1:15)
  save(seu.obj,file = f)
}
load(f)
ElbowPlot(seu.obj)

plot_zoom3.png

p1 <- DimPlot(seu.obj, reduction = "umap",label = T)+NoLegend();p1

plot_zoom4.png

標(biāo)記#的兩個地方可能會需要改動

一個是dims = 1:15,代表選擇前15個維度，15是個比較折中的值拗窃，選擇的數(shù)量越多計算量越大瞎领，而且信息越冗余；選的太少又不具有代表性随夸。一般是10-20九默，根據(jù)ElbowPlot來選擇拐點的橫坐標(biāo)(拐點就是在這個點之前縱坐標(biāo)下降比較快，在這個點之后縱坐標(biāo)下降比較慢)宾毒。不是很重要驼修，只要拐點在15之前，選15就一點問題都沒有诈铛。

一個是resolution = 0.5乙各，代表分群的分辨率是0.5，這個參數(shù)的選擇范圍是0.1-1.5之間幢竹，數(shù)值越大分出來的群越多耳峦，數(shù)值越小分出來的點越少。0.5也是一個比較折中的值妨退，如果后面注釋不困難妇萄，就不用改動。

如果決定要改動那么就不能只改代碼咬荷，要把obj.Rdata從工作目錄下刪掉冠句，再重新運行修改后的代碼。