生物信息學習的正確姿勢

NGS系列文章包括NGS基礎、在線繪圖、轉錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫學設計實驗GEO數據分析 (step-by-step)）、批次效應處理等內容。

哈佛大學單細胞課程|筆記匯總 （五）

哈佛大學單細胞課程|筆記匯總 （四）

（六）Single-cell RNA-seq clustering analysis: aligning cells across conditions

我們的數據集包含來自兩個不同條件（對照和刺激）的兩個樣本，因此將這些樣本整合在一起以更好地進行比較。

聚類工作流程

聚類分析的目的是在待定義細胞類型的數據集中保留主要的變異來源，同時盡量屏蔽由于無用的變異來源（測序深度，細胞周期差異，線粒體表達，批次效應等）而產生的變異。

在我們的工作流程中需要應用到以下兩個資源：

• Satija Lab: Seurat v3 Guided Integration Tutorial(https://satijalab.org/seurat/v3.0/immune_alignment.html)

• Paul Hoffman: Cell-Cycle Scoring and Regression(http://satijalab.org/seurat/cell_cycle_vignette.html)

為了識別細胞亞群，我們將進行以下步驟：

（1）Normalization, variance stabilization, and regression of unwanted variation (e.g. mitochondrial transcript abundance, cell cycle phase, etc.) for each sample

（2）Integrationof the samples using shared highly variable genes (optional, but recommended to align cells from different samples/conditions if cell types are separating by sample/condition)

（3）Clustering cells based on top PCs (metagenes)

（4）Exploration of quality control metrics: determine whether clusters are unbalanced wrt UMIs, genes, cell cycle, mitochondrial content, samples, etc.

（5）Searching for expected cell types using known cell type-specific gene markers

以下流程為前兩步。

加載R包

# Single-cell RNA-seq analysis - clustering analysis# Load libraries library(Seurat) library(tidyverse) library(RCurl) library(cowplot)

Normalization, variance stabilization, and regression of unwanted variation for each sample

分析的第一步是將count矩陣標準化，以解決每個樣品每個細胞的測序深度差異。Seurat最近推出了一種叫sctransform的新方法，用于對scRNA-seq數據進行歸一化和方差穩定變化。

sctransform方法使用正則化負二項式模型對UMI counts進行建模，以消除由于測序深度（每個細胞的總nUMI）所引起的變異，同時基于具有相似豐度的基因之間的合并信息來調整方差。

Image credit: Hafemeister C and Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression, bioRxiv 2019 (https://doi.org/10.1101/576827)

模型結果殘差（residuals）是每個轉錄本標準化后的表達水平。

sctransform能通過回歸分析校對測序深度（nUMI）。但對于某些數據集，細胞周期可能是基因表達變化的顯著來源，對于其他數據集則不一定，我們需要檢查細胞周期是否是數據變化的主要來源并在需要時校對細胞周期的影響。

Cell cycle scoring

建議在執行sctransform方法之前檢查細胞周期階段。檢查之前先對數據做個標準化使得不同測序深度的細胞可比。

# Normalize the counts seurat_phase <- NormalizeData(filtered_seurat)

作者提供了人體細胞周期相關基因可以下載 (https://www.dropbox.com/s/hus4mrkueh1tfpr/cycle.rda?dl=1)和其他生物的細胞周期相關基因（https://github.com/hbctraining/scRNA-seq/blob/master/lessons/cell_cycle_scoring.md）。

把細胞周期相關基因讀入data目錄

# Load cell cycle markers load("data/cycle.rda")# Score cells for cell cycle seurat_phase <- CellCycleScoring(seurat_phase,g2m.features = g2m_genes,s.features = s_genes)# View cell cycle scores and phases assigned to cells View(seurat_phase@meta.data)

在對細胞進行細胞周期評分后，我們想使用PCA確定細胞周期是否是我們數據集中變異的主要來源。在此之前，我們首先需要對數據進行標準化（scale），并使用Seurat ScaleData（）函數。

• 調整每個基因的表達以使整個細胞的平均表達為0

• 縮放每個基因的表達以使細胞之間的表達值轉換標準差的倍數

# Identify the most variable genes seurat_phase <- FindVariableFeatures(seurat_phase,selection.method = "vst",nfeatures = 2000,verbose = FALSE)# Scale the counts seurat_phase <- ScaleData(seurat_phase)

NOTE: For the selection.method and nfeatures arguments the values specified are the default settings. Therefore, you do not necessarily need to include these in your code. We have included it here for transparency and inform you what you are using.

# Perform PCA seurat_phase <- RunPCA(seurat_phase)# Plot the PCA colored by cell cycle phase DimPlot(seurat_phase,reduction = "pca",group.by= "Phase",split.by = "Phase")

我們沒有看到在細胞周期方面具有很大的差異，也基于此圖，我們不會考慮細胞周期引起的變化。

SCTransform

我們此時可以放心的進行SCTransform，我們打算使用for循環對每個樣本進行 NormalizeData(), CellCycleScoring()和SCTransform() ，并使用 vars.to.regress 去除線粒體相關基因。

# Split seurat object by condition to perform cell cycle scoring and SCT on all samples split_seurat <- SplitObject(filtered_seurat, split.by = "sample")split_seurat <- split_seurat[c("ctrl", "stim")]for (i in 1:length(split_seurat)) {split_seurat[[i]] <- NormalizeData(split_seurat[[i]], verbose = TRUE)split_seurat[[i]] <- CellCycleScoring(split_seurat[[i]], g2m.features=g2m_genes, s.features=s_genes)split_seurat[[i]] <- SCTransform(split_seurat[[i]], vars.to.regress = c("mitoRatio"))}

NOTE: By default, after normalizing, adjusting the variance, and regressing out uninteresting sources of variation, SCTransform will rank the genes by residual variance and output the 3000 most variant genes. If the dataset has larger cell numbers, then it may be beneficial to adjust this parameter higher using the variable.features.n argument.

除了原始RNA計數外，assays slot中現在還有一個SCT組件。細胞之間表達變化最大的特征基因存儲于SCT assay中。

通常，在決定是否需要進行數據整合之前，我們總是先查看細胞在空間的分布。如果我們在Seurat對象中同時對兩個條件進行了歸一化并可視化，我們將看到特定于不同條件的聚類：

不同條件下細胞的聚類表明我們需要跨條件整合細胞。

NOTE: Seurat has a vignette for how to run through the workflow without integration（https://satijalab.org/seurat/v3.1/sctransform_vignette.html）. The workflow is fairly similar to this workflow, but the samples would not necessarily be split in the beginning and integration would not be performed.

Integrate samples using shared highly variable genes

為了進行整合（“harmony”整合不同平臺的單細胞數據之旅），我們將使用不同條件的細胞共有的高可變基因，然后，我們將“整合”或“協調”條件，以鑒定組之間相似或具有“共同的生物學特征”的細胞。

如果不確定要期望的簇或期望條件之間的某些不同細胞類型（例如腫瘤和對照樣品），可以先單獨處理某個條件，然后一起運行以查看兩種條件下是否存在針對特定細胞類型的特定簇。通常在不同條件下進行聚類分析時會出現基于條件的聚類，而數據整合可以幫助相同類型的細胞進行聚類。為了整合，我們首先需要先通過SCTransform篩選高變基因，然行“整合”或者“協調”（"integrate" or "harmonize"）條件，以便在不同條件下同種細胞具有相似的表達。這些條件包括：

（1）不同條件（control VS stimuli）

（2）不同數據集

（3）不同類型 (e.g. scRNA-seq and scATAC-seq)

整合最重要的目的就是使得不同條件和數據集之間細胞類型相同的數據表達譜一致。這也意味著部分細胞亞群的生物學狀態是相同的，下圖為整合步驟：

Image credit: Stuart T and Butler A, et al. Comprehensive integration of single cell data, bioRxiv 2018 (https://doi.org/10.1101/460147)

具體細節是：

先執行典型相關性分析（canonical correlation analysis (CCA)）：

CCA是PCA的一種形式，能識別數據中最大的差異，不過前提是組/條件之間共有這個最大變異（使用每個樣本中變化最大的3000個基因）。

NOTE: The shared highly variable genes are used because they are the most likely to represent those genes distinguishing the different cell types present.

在數據集中識別錨點（anchors）或相互最近的鄰居（mutual nearest neighbors ，MNN）：MMN可以被認為是“最佳搭檔”（'best buddies'）。對于每種條件的細胞來說：

“The difference in expression values between cells in an MNN pair provides an estimate of the batch effect, which is made more precise by averaging across many such pairs. A correction vector is obtained and applied to the expression values to perform batch correction.” ——Stuart and Bulter et al. (2018)（https://www.biorxiv.org/content/early/2018/11/02/460147）.

• 根據基因表達值確定細胞在其他情況下最接近的鄰居-它是'best buddies'。

• 進行二次分析，如果兩個細胞在兩個方向上都是'best buddies'，則這些細胞將被標記為將兩個數據集“錨定”在一起的錨點。

• 篩選錨點：通過其本地鄰域中的重疊來評估錨點對之間的相似性（不正確的錨點得分會較低）——相鄰細胞是否具有彼此相鄰的'best buddies'？

• 整合數據/條件：使用錨點和相應分數來轉換細胞表達值，從而可以整合條件/數據集（不同的樣本、條件、數據集、模態）

• NOTE: Transformation of each cell uses a weighted average of the two cells of each anchor across anchors of the datasets. Weights determined by cell similarity score (distance between cell and k nearest anchors) and anchor scores, so cells in the same neighborhood should have similar correction values.

  現在，使用我們的SCTransform對象作為輸入，進行數據整合：

  首先，我們需要指定使用由SCTransform識別的3000個高變基因進行整合。默認情況下，僅選擇前2000個基因。

  # Select the most variable features to use for integration integ_features <- SelectIntegrationFeatures(object.list = split_seurat,nfeatures = 3000)

  然后，我們需要準備SCTransform對象以進行整合。

  # Prepare the SCT list object for integration split_seurat <- PrepSCTIntegration(object.list = split_seurat,anchor.features = integ_features)

  現在，我們將執行canonical correlation analysis（CCA），找到最佳錨點并過濾不正確的錨點（Cell 深度 一套普遍適用于各類單細胞測序數據集的錨定整合方案）。

  # Find best buddies - can take a while to run integ_anchors <- FindIntegrationAnchors(object.list = split_seurat,normalization.method = "SCT",anchor.features = integ_features)

  最后，整合不同條件并保存。

  # Integrate across conditions seurat_integrated <- IntegrateData(anchorset = integ_anchors,normalization.method = "SCT")# Save integrated seurat object saveRDS(seurat_integrated, "results/integrated_seurat.rds")

  UMAP visualization

  我們利用降維技術可視化整合后的數據，例如PCA和Uniform Manifold Approximation and Projection（UMAP）。

  盡管PCA將能獲得所有PC，但繪圖時一般只使用兩個或3個。相比之下，UMAP則進一步整合多個PCs獲取更多信息，并在二維中進行展示。這樣，細胞之間的距離代表表達上的相似性。

  為了生成這些可視化，我們需要先運行PCA和UMAP方法。先從PCA開始：

  # Run PCA seurat_integrated <- RunPCA(object = seurat_integrated)# Plot PCA PCAPlot(seurat_integrated,split.by = "sample")

  我們看到在PCA映射中具有很好的重合。然后進行UMAP：

  # Run UMAP seurat_integrated <- RunUMAP(seurat_integrated,dims = 1:40,reduction = "pca")# Plot UMAP DimPlot(seurat_integrated)

  整合數據后

  整合數據前

  從圖中我們可以看到很好的整合，比起未經過CCA的數據要好太多了。

  往期精品(點擊圖片直達文字對應教程)

  后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結

以上是生活随笔為你收集整理的哈佛大学单细胞课程|笔记汇总 （六）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。