10X genomics的基本原理
在這個教程中，主要將分析 10X Genomics 免費(fèi)提供的外周血單核細(xì)胞 (PBMC) 數(shù)據(jù)集配阵。在 Illumina NextSeq 500 上對 2,700 個單細(xì)胞進(jìn)行了測序】ㄇ可以在https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz此處找到原始數(shù)據(jù)桂敛。我們從讀取數(shù)據(jù)開始赡矢。 Read10X() 函數(shù)從 10X 讀取 cellranger 管道的輸出，返回一個唯一的分子識別 (UMI) 計數(shù)矩陣饱须。此矩陣中的值表示在每個單元格（列）中檢測到的每個特征（即基因域醇；行）的分子數(shù)。
Read10X() 函數(shù)是針對于整理好的10X Genomics 數(shù)據(jù)蓉媳，如果手頭的不是類似文件譬挚，可以將其進(jìn)行轉(zhuǎn)換，成為格式一致的文件酪呻。
接下來使用計數(shù)矩陣創(chuàng)建一個 Seurat 對象减宣。該對象用作包含單細(xì)胞數(shù)據(jù)集的數(shù)據(jù)（如計數(shù)矩陣）和分析（如 PCA 或聚類結(jié)果）的容器。例如玩荠，count matrix存儲在 pbmc[[“RNA”]]@counts 中漆腌。
library(dplyr)
library(Seurat)
library(patchwork)
創(chuàng)建對象
加載數(shù)據(jù)
# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
創(chuàng)建 Seurat 對象
### 2.創(chuàng)建Seurat對象
### counts 輸入的是數(shù)據(jù)，行是基因阶冈，列是細(xì)胞
### project 參數(shù)輸入的是項目名稱,出現(xiàn)在metadata的orig.ident這一列
### min.cells 限定的是基因：每個基因在至少多少個細(xì)胞中出現(xiàn)
### min.features 限定的是細(xì)胞: 每個細(xì)胞中最少有多少個基因
scobj <- CreateSeuratObject(counts = scdata,
project = "pbmc3k",
min.cells = 3,
min.features = 200)
count matrix是什么樣子?
count矩陣是稀松矩陣闷尿，可以減少占用空間
pbmc.data[c("IGF2BP2", "TCL1A", "MS4A1"), 1:30]
dense.size <- object.size(as.matrix(pbmc.data))
dense.size
sparse.size <- object.size(pbmc.data)
sparse.size
dense.size/sparse.size
預(yù)處理流程
計算線粒體含量
這是質(zhì)控的重要步驟，使用PercentageFeatureSet函數(shù)
### 主要PercentageFeatureSet函數(shù)計算線粒體含量
### 人類使用pattern = "^MT-"女坑，小鼠使用pattern = "^mt-"
scobj[["percent.mt"]] <- PercentageFeatureSet(scobj, pattern = "^MT-")
### 該操作會在metadata數(shù)據(jù)里面增加一列叫做percent.mt
metadata <- scobj@meta.data
一般情況下填具，可以認(rèn)為線粒體含量多，意味著細(xì)胞可能趨于死亡匆骗，這樣的細(xì)胞就應(yīng)該剔除灌旧。但是如果本身研究的就是和線粒體相關(guān)的內(nèi)容绑咱，例如藥物干預(yù)會引起線粒體的變化，那么就要根據(jù)具體情況來分析是否需要剔除枢泰。
質(zhì)控展示
### 質(zhì)控數(shù)據(jù)可視化描融，使用VlnPlot函數(shù)
### nFeature_RNA, number of Feature, 每個細(xì)胞中有多少個基因
### nCount_RNA, number of counts, 每個細(xì)胞中有多少個counts
### percent.mt, 我們自己增加的列, 每個細(xì)胞中線粒體基因的比例
VlnPlot(scobj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
大部分情況下，線粒體的含量應(yīng)該去除較大的衡蚂，這個可以根據(jù)具體的情況進(jìn)行分析窿克，一般10%以內(nèi)可以接受。但是依然是要根據(jù)相應(yīng)的干預(yù)和特點(diǎn)進(jìn)行取舍毛甲。
也可以展示特征之間的關(guān)系
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
正式篩選
這一步依然是根據(jù)metadata的數(shù)據(jù)進(jìn)行篩選年叮，這一部分依然是根據(jù)相關(guān)的具體情況進(jìn)行篩選。
### 正式篩選玻募，篩選的是細(xì)胞只损，最終細(xì)胞減少
### nFeature_RNA > 200
### nFeature_RNA < 2500
### percent.mt < 5
scobj <- subset(scobj, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
篩選之后，再次用小提琴圖展示七咧，是上面小提琴圖的縮減版本跃惫。
至于到底應(yīng)該確定到多少的篩選指標(biāo)，可以先按照常規(guī)指標(biāo)進(jìn)行設(shè)定艾栋，在后續(xù)的細(xì)胞聚類中爆存，看看是否有與線粒體相關(guān)的群聚類出來，然后再更改篩選指標(biāo)蝗砾，聚類指標(biāo)確實消失先较，那么就可以對篩選指標(biāo)確定。
數(shù)據(jù)標(biāo)準(zhǔn)化
每個細(xì)胞分別進(jìn)行檢測悼粮，所以如果要比較各個細(xì)胞闲勺，還是要進(jìn)行標(biāo)準(zhǔn)化
### 先除以總數(shù)，再乘以系數(shù)扣猫，然后取log
scobj <- NormalizeData(scobj, normalization.method = "LogNormalize", scale.factor = 10000)
### 默認(rèn)參數(shù)可以省略
scobj <- NormalizeData(scobj)
特征篩選
如果某些基因的表達(dá)再各個細(xì)胞之間表達(dá)是恒定的菜循，所以要區(qū)分各個細(xì)胞，就要使用變化差異大的基因苞笨。
尋找高變基因
scobj <- FindVariableFeatures(scobj, selection.method = "vst", nfeatures = 2000)
### 默認(rèn)參數(shù)可以省略
scobj <- FindVariableFeatures(scobj)
高變基因繪圖
### 使用VariableFeatures 函數(shù)提取高變基因
### 等同于 scobj@assays$RNA@var.features
top10 <- head(VariableFeatures(scobj), 10)
### 使用VariableFeaturePlot 畫圖
plot1 <- VariableFeaturePlot(scobj)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
縮放數(shù)據(jù)
縮放的數(shù)據(jù)達(dá)到的效果是：基因的平均值為0债朵，方差為1后专。便于不同類型之間的比較吏颖。
### 降維之前的必備操作
scobj <- ScaleData(scobj, features = rownames(scobj))
### 如果不限定參數(shù)，只會縮放高變基因
### scobj <- ScaleData(scobj)
### 縮放后的數(shù)據(jù)存放在scobj@assays$RNA@scale.data,會很大
scale.data <- scobj@assays$RNA@scale.data
這一步縮放數(shù)據(jù)會花費(fèi)比較長的時間蔬咬，原因是這一步的操作是在所有基因上展開粤咪。但是谚中，其實我們最需要的是變異大的基因，也就是可能是標(biāo)志基因的數(shù)據(jù)。所以宪塔，這個函數(shù)默認(rèn)是可以使用前2000的基因進(jìn)行縮放的磁奖，可以省略features參數(shù)的內(nèi)容。
scobj <- ScaleData(scobj)
需要說明的是某筐，雖然縮放會花費(fèi)一些時間比搭，數(shù)據(jù)也會變大，但是南誊，因為后續(xù)如果要繪制熱圖身诺，使用的依然是縮放數(shù)據(jù)，所以依然在這一步將所有基因的數(shù)據(jù)進(jìn)行縮放抄囚。
如果霉赡，保存數(shù)據(jù)的時候，還是可以把縮放數(shù)據(jù)去除后幔托，再進(jìn)行保存穴亏。下次需要的時候再次縮放即可
可以這樣操作
scobj@assays$RNA@scale.data <-matrix()
PCA線性降維
PCA操作的對象為縮放過的數(shù)據(jù)，因此重挑，需要保證前面的縮放數(shù)據(jù)嗓化。結(jié)果會保存在對象的reduction中，可以從中提取數(shù)據(jù)攒驰。
為什么要進(jìn)行主成分分析蟆湖，主要還是用少量的數(shù)據(jù)和維度故爵，盡可能的保留原來數(shù)據(jù)特點(diǎn)玻粪。
降維
scobj <- RunPCA(scobj, features = VariableFeatures(object = scobj))
DimPlot(scobj, reduction = "pca")
### PCA降維數(shù)據(jù)存放在scobj@reductions$pca中
pcadata = as.data.frame(scobj@reductions$pca@cell.embeddings)
ggplot(pcadata,aes(PC_1,PC_2,color="red"))+
geom_point()
### 選擇合適的PCA維度
ElbowPlot(scobj)
降維展示
# Examine and visualize PCA results a few different ways
print(scobj[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(scobj, dims = 1:2, reduction = "pca")
[圖片上傳失敗...(image-b45585-1685178814476)]
使用熱圖進(jìn)行展示
DimHeatmap(scobj, dims = 1, cells = 500, balanced = TRUE)
dims 是確定需要展示維度數(shù)，可以嘗試
[圖片上傳失敗...(image-f3b1c3-1685178814476)]
UMAP非線性降維
如何將高維數(shù)據(jù)壓縮到二維平面中诬垂，使用UMAP進(jìn)行處理劲室。其原理可以如下展示：
Understanding UMAP
### 依賴PCA的結(jié)果
### dims = 1:10 由上一步確定
scobj <- RunUMAP(scobj, dims = 1:10)
DimPlot(scobj, reduction = "umap")
### UMAP降維數(shù)據(jù)存放在scobj@reductions$umap中
umapdata = as.data.frame(scobj@reductions$umap@cell.embeddings)
ggplot(umapdata,aes(UMAP_1,UMAP_2,color="red"))+
geom_point()
聚類分群
二維空間中看到了細(xì)胞分群，那么在這個基礎(chǔ)上结窘，怎么進(jìn)行亞群分類很洋，這就是這一步的意義。
### 找緊鄰,dims = 1:10 跟UMAP相同
scobj <- FindNeighbors(scobj, dims = 1:10)
### 分群
scobj <- FindClusters(scobj, resolution = 0.5)
### 會在metadata中增加兩列數(shù)據(jù)"RNA_snn_res.0.5" "seurat_clusters"
metadata <- scobj@meta.data
resolution是這個主要參數(shù)隧枫，也就是亞群切割的分辨率喉磁，這個具體怎么挑選，可以按照以下的方法可視化查看官脓。
### 設(shè)置多個resolution選擇合適的resolution
scobj <- FindClusters(scobj, resolution = seq(0.2,1.2,0.1))
metadata <- scobj@meta.data
library(clustree)
clustree(scobj)
根據(jù)上面的可視化結(jié)果协怒，可以大致確定一個分辨率
### 選擇特定分辨率得到的分群此處為RNA_snn_res.0.5
scobj@meta.data$seurat_clusters <- scobj@meta.data$RNA_snn_res.0.5
Idents(scobj) <- "seurat_clusters"
DimPlot(scobj, reduction = "umap", label = T)
[圖片上傳失敗...(image-e8e7a6-1685178814476)]
大致得到了相應(yīng)的分群，后續(xù)就需要對這些分群進(jìn)行注釋卑笨。注釋是一個與個人知識背景很相關(guān)的分析過程孕暇，后續(xù)補(bǔ)充。
參考文章
Seurat - Guided Clustering Tutorial
Visualization of gene expression with Nebulosa (in Seurat)
Use regularized negative binomial regression to normalize UMI count data
