###### tags: `single-cell RNA-seq` `Seurat`
Seuratによる10Xデータの解析(基本編)
===
* Cell Rangerでの一次解析まで終わっている前提
* Spermatogenesisの４つのサンプル(vivo9, vitro9, vivo14, vitro14)を例に、一度に読み込んでQC,データの統合、UMAPでの次元圧縮まで行う。
* ここで定義する関数は`source("/osc-fs_home/t-suzuki/TenX_invitro_spermatogenesis/Functions.r")`でまとめて読み込める。
2020/07/18 Update
* QCプロットのドットを小さくしてバイオリンプロットを見やすくした。
* 3rd analysisに合わせてQCのcut-off値を変更
### パッケージの読み込み
```R
library(Seurat)
library(dplyr)
```
## データの読み込み
Cell Rangerでの解析結果フォルダ中のデータをSeurat objectとして読み込み
関数を作成して読み込みを自動化
### 関数の定義
```R
#--------------------------------------------
# scCountProc
#--------------------------------------------
#' Reading cellranger processed 10x data as a seurat object.
#'
#' Read the cell ranger processed 10x data, assign a sample name, normalize, and find veriable genes
#'
#' @param data_path path of ceranger preocessed result (one upper from "outs" directory)
#' @param sample_name charector a sample name
#' @param normalization logical, if TRUE, normalization is performed
#' @param nfetures integer, number of variable genes to be found
#'
#' @importFrom Seurat Read10X CreateSeuratObject RenameCells NormalizeData FindVariableFeatures
#'
#' @return a seurat object
#'
#' @keywords seurat 10x
#' @export
scCountProc <- function(data_path, sample_name = sample_name, normalization = FALSE, nfeatures =2000){
## Count
library(Seurat)
data.counts <- Read10X(data.dir = paste0(data_path, "/outs/filtered_feature_bc_matrix/"))
data <- CreateSeuratObject(counts = data.counts) # conversion to Seurat object
data$sample <- sample_name #サンプル名の付加
data <- RenameCells(object=data, add.cell.id = sample_name)
if (normalization == TRUE){
data <- NormalizeData(object = data) #正規化
}
data <- FindVariableFeatures(object = data, selection.method = "vst", nfeatures = nfeatures) #変動の大きい遺伝子の抽出
return(data)
}
```
### `scCountProc`を使用してデータの読み込み
```r
data.vivo09 <- scCountProc(
data_path = "/osc-fs_home/t-suzuki/TenX_invitro_spermatogenesis/vivo9_count", #Cell Rangerの結果のディレクトリ
sample_name = "vivo9", #サンプル名（任意）
normalization = FALSE Normalizeは後でやるのでここではFALSE
)
data.vitro09 <- scCountProc(
data_path = "/osc-fs_home/t-suzuki/TenX_invitro_spermatogenesis/vitro9_count",
sample_name = "vitro9",
normalization = FALSE
)
data.vivo14 <- scCountProc(
data_path = "/osc-fs_home/t-suzuki/TenX_invitro_spermatogenesis/vivo14_count",
sample_name = "vivo14",
normalization = FALSE
)
data.vitro14 <- scCountProc(
data_path = "/osc-fs_home/t-suzuki/TenX_invitro_spermatogenesis/vitro14_count",
sample_name = "vitro14",
normalization = FALSE
)
```
### データをリスト化
```r
data.list <- list(vivo09=data.vivo09, vitro09=data.vitro09, vivo14=data.vivo14, vitro14=data.vitro14)
#個別のオブジェクトは以降必要ないので消去（メモリ節約のため）
rm(list=c("data.vivo09", "data.vitro09", "data.vivo14", "data.vitro14"))
```
## QC
それぞれのデータのそれぞれの細胞について、1. 検出遺伝子数、2 ミトコンドリアゲノムRNAの割合をもとに、クオリティーの低い細胞を除去。
### それぞれの情報を可視化するためのplot objectを作成
```r
vp_list <- list() #violin plotのリストを定義
combp_list <- list() #scatter plotのリストを定義
for (i in 1:length(data.list)) {
#全体の分子数中のミトコンドリア遺伝子の割合を計算
data.list[[i]][["percent.mt"]] <- PercentageFeatureSet(object = data.list[[i]], pattern = "^mt-")
#violin plotのオブジェクト作成
vp_list <- c(vp_list,
list(VlnPlot(data.list[[i]],
features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
pt.size = 0.1,
ncol = 3)
)
)
#scatter plotのオブジェクト作成
plot1 <- FeatureScatter(data.list[[i]],
feature1 = "nCount_RNA",
feature2 = "percent.mt"
)
plot2 <- FeatureScatter(data.list[[i]],
feature1 = "nCount_RNA",
feature2 = "nFeature_RNA"
)
combp_list <- c(combp_list, list(CombinePlots(plots = list(plot1, plot2))))
}
```
### プロット
```r
vp_list
combp_list
```
#### 例）vivo9のプロット
このプロットをもとに検出遺伝子数が多すぎもなく少なすぎもなく、ミトコンドリアゲノムの割合が多すぎではないものだけを**マニュアル**で取ってくる。（結構適当）
### violin plotとscatter plotをみてcut-offを決めて範囲内の細胞を抽出
```r
data.list[["vivo09"]] <- subset(data.list[["vivo09"]], subset = nFeature_RNA > 300 & nFeature_RNA < 6500 & percent.mt < 75 )
data.list[["vitro09"]] <- subset(data.list[["vitro09"]], subset = nFeature_RNA > 300 & nFeature_RNA < 6000 & percent.mt < 75)
data.list[["vivo14"]] <- subset(data.list[["vivo14"]], subset = nFeature_RNA > 300 & nFeature_RNA < 6500 & percent.mt < 75)
data.list[["vitro14"]] <- subset(data.list[["vitro14"]], subset = nFeature_RNA > 300 & nFeature_RNA < 6500 & percent.mt < 75)
```
## Normalizationと変動遺伝子の検出
```r
for (i in 1:length(data.list)) {
data.list[[i]] <- NormalizeData(object = data.list[[i]], verbose = FALSE)
data.list[[i]] <- FindVariableFeatures(object = data.list[[i]], selection.method = "vst", nfeatures = 2000, verbose = FALSE)
}
```
## Seurat objectを統合
```r
data.anchors <- FindIntegrationAnchors(object.list = data.list, dims = 1:30)
all_genes <- data.anchors@object.list[[1]]@assays$RNA@counts@Dimnames[[1]]
data.integrated <- IntegrateData(anchorset = data.anchors, dims = 1:30, features.to.integrate=all_genes)
## デフォルトをintegrated assayに変更
DefaultAssay(object = data.integrated) <- "integrated"
## 一度ここで保存
save(data.integrated, file="data.integrated.QC.RData")
```
## Cell-Cycleの影響を調整
同じ細胞種でもcell cycleの影響で異なるクラスターに分かれてしまう可能性があるので、必要であれば、G1, G2/M, S期の影響を調整する。ただし、分化途中の細胞などを解析する際は必ずしも調整しない方がよい場合もあるのでケースバイケース。Spermatogenesisの場合はG2/MとSの影響のみを調整。
### cell cycle markerリスト(Tirosh et al, 2015)の読み込み
```r
s.genes <- sapply(cc.genes.updated.2019$s.genes, function(x){paste0(substr(x,1,1), tolower(substr(x,2,nchar(x))))})
g2m.genes <- sapply(cc.genes.updated.2019$g2m.genes, function(x){paste0(substr(x,1,1), tolower(substr(x,2,nchar(x))))})
```
### PCAでCell-Cycleの影響の確認
Cell cycleの補正前はG1, S, G2/Mに偏りがある
```r
# データ全体をscaling
data.integrated <- ScaleData(object =data.integrated, features = rownames(data.integrated), verbose = FALSE)
# PC
data.integrated <- RunPCA(data.integrated, features = c(s.genes, g2m.genes))
DimPlot(data.integrated)
```
G1, G2/M, Sで分かれている。
### cell-cycle スコアリング
G1, G2/M, Sをそれぞれ特異的遺伝子に基づいてスコアリングする
```r
data.integrated <- CellCycleScoring(data.integrated, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
```
### Cell-Cycleの影響の調整（半日から１日くらい時間がかかる)
G1, G2/M, Sすべてを調整する(分化途中の細胞などを見る場合はG2/M, Sのみの調整)
```r
data.integrated <- ScaleData(data.integrated, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(data.integrated))
```
(またはG2/M、Sのみの調整の場合。今回はこちらを採用)
```r
data.integrated$CC.Difference <- data.integrated$S.Score - data.integrated$G2M.Score
data.integrated <- ScaleData(data.integrated, vars.to.regress = "CC.Difference", features = rownames(data.integrated))
```
### Cell-Cycleの調整の確認
```r
data.integrated <- RunPCA(data.integrated, features = c(s.genes, g2m.genes))
DimPlot(data.integrated)
## ここで一度save
save(data.integrated, file="data.integrated.QC.ccreg.RData")
```
G2/M, Sの偏りが改善
## クラスタリングとUMAPでの次元圧縮・可視化
---
### デフォルトをintegrated assayに変更
```r
DefaultAssay(object = data.integrated) <- "integrated"
```
### Dimensionalityの決定
クラスタリングをする際に、何個目までの主成分を使用するかを決める。JackStraw解析をQQ-plotしたものから一様分布の累積分布関数（点線）と一致する一つ手前の次元数、またはElbow Plotで水平になる付近の次元数を**目視で**みて次元数を決める。
```r
data.integrated <- JackStraw(data.integrated, num.replicate = 100)
data.integrated <- ScoreJackStraw(data.integrated, dims = 1:20)
JackStrawPlot(data.integrated, dims = 1:20)
ElbowPlot(data.integrated, ndims = 50)
```
Jack strewで一様分布のp-value（点線）と一致する一つ前のPC、もしくはelbowplotのelbowのPCまでをクラスタリングのdimensionalityとして使用する。
ここではJackStrawPlotからPC14までとした。
### クラスタリング
```r
data.integrated <- RunPCA(object = data.integrated, npcs = 20, verbose = FALSE)
data.integrated <- RunUMAP(object = data.integrated, reduction = "pca", dims = 1:14) #ここのdimsをJackStrawPlotかElbowPlotの結果で調整
data.integrated <- FindNeighbors(object = data.integrated, dims = 1:14) #ここのdimsをJackStrawPlotかElbowPlotの結果で調整
data.integrated <- FindClusters(object = data.integrated, algorithm = 1, resolution = 0.1)
```
### umapで可視化
```r
DimPlot(object = data.integrated, reduction = "umap")
```
umapの結果を見ながら、FindClustersのresolutionを調整しte適度なクラスタリングになるまで最適化する。resolutionがちいさいほど細かくクラスタリングされる。
Sign in with Wallet
Connect another wallet