scRNA-seqやってみる（その3）

時系列や各種処理を与えたサンプル間で比較するなら、それらをマージして同じクラスターで重ね合わせて考えたい。
まずはumapチャートを出力するまで。4つの実験区の比較を行う例をやってみたい。Seuratのハンズオンでは2区の重ね合わせなんだけど、実際に10xの出力として保存されたデータからではなく、ライブラリからサンプルデータを読み出して処理しており、実際的ではない。なのでCellrangerの出力するfiltered_feature_bc_matrixフォルダに3つのファイル(barcode.tsv.gz, features.tsv.gz, matrix.mtx.gz)が保存されているところをスタートとした解析を試してみる。

if(!require(Seurat)){install.packages("Seurat")}
if(!require(patchwork)){install.packages("patchwork")}
if(!require(dplyr)){install.packages("dplyr")}

library(Seurat)
library(patchwork)

install.packages('BiocManager')
# Bioconductorのバージョンを設定
BiocManager::install(version = "3.19")
BiocManager::install('glmGamPoi')

#####Remove all objects before running#####
rm(list = ls())

input.dir <- "/path/to/sample/"

# サンプル1のデータセットを読み込む
sample1 <- Read10X(paste(input.dir, "sample1/filtered_feature_bc_matrix", sep=""))

# サンプル2のデータセットを読み込む
sample2 <- Read10X(paste(input.dir, "sample2/filtered_feature_bc_matrix", sep=""))

# サンプル3のデータセットを読み込む
sample3 <- Read10X(paste(input.dir, "sample3/filtered_feature_bc_matrix", sep=""))

# サンプル4のデータセットを読み込む
sample4 <- Read10X(paste(input.dir, "sample4/filtered_feature_bc_matrix", sep=""))

sample1 <- CreateSeuratObject(counts = sample1, project = "sample1")
sample2 <- CreateSeuratObject(counts = sample2, project = "sample2")
sample3 <- CreateSeuratObject(counts = sample3, project = "sample3")
sample4 <- CreateSeuratObject(counts = sample4, project = "sample4")

# 各データセットにユニークなセル識別子を追加
sample1 <- RenameCells(sample1, new.names = paste("sample1", Cells(sample1), sep = "_"))
sample2 <- RenameCells(sample2, new.names = paste("sample2", Cells(sample2), sep = "_"))
sample3 <- RenameCells(sample3, new.names = paste("sample3", Cells(sample3), sep = "_"))
sample4 <- RenameCells(sample4, new.names = paste("sample4", Cells(sample4), sep = "_"))

# サンプルデータセットのリストを作成
datasets <- list(sample1, sample2, sample3, sample4)

# データセットの前処理
datasets <- lapply(datasets, function(x) {
  x <- SCTransform(x, verbose = FALSE)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
  x
})

# 統合のためのアンカーフィーチャーの選択
features <- SelectIntegrationFeatures(object.list = datasets, nfeatures = 2000)

# SCT正規化を使用した統合アンカーの発見
datasets <- PrepSCTIntegration(object.list = datasets, anchor.features = features)
anchors <- FindIntegrationAnchors(object.list = datasets, normalization.method = "SCT", anchor.features = features)

# データセットの統合
combined <- IntegrateData(anchorset = anchors, normalization.method = "SCT")

# 統合データの次の解析ステップ
combined <- ScaleData(combined, verbose = FALSE)
combined <- RunPCA(combined, npcs = 30, verbose = FALSE)
combined <- RunUMAP(combined, reduction = "pca", dims = 1:30)
combined <- FindNeighbors(combined, reduction = "pca", dims = 1:30)
combined <- FindClusters(combined, resolution = 0.5)

# ここでデータを一旦保存しておく
# 後で読み込むときはcombined <- readRDS(file =paste(input.dir, "combined.rds", sep=""))
saveRDS(combined, file = paste(input.dir, "combined.rds", sep=""))

# 結果の可視化
p1 <- DimPlot(combined, reduction = "umap", group.by = "orig.ident")
p2 <- DimPlot(combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

DimPlot(combined, reduction = "umap", split.by = "orig.ident")

# 出来上がったクラスターに名前をつけることができる
combined <- RenameIdents(combined, `0` = "Cluster 0", `1` = "Cluster 1", `2` = "Cluster 2", `3` = "Cluster 3")
DimPlot(combined, label = TRUE)

# マーカー遺伝子の可視化
markers.to.plot <- c("GeneA", "GeneB", "GeneC")
DotPlot(combined, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8, split.by = "orig.ident") + RotatedAxis()

# Find Differentially Expressed Genes:
fold_change_threshold <- 0.25
pct_threshold <- 0.25
p_value_threshold <- 0.01

# クラスタごとのDEG
all_markers <- FindAllMarkers(
  object = combined,
  only.pos = FALSE,
  min.pct = pct_threshold,
  logfc.threshold = fold_change_threshold,
  test.use = "wilcox"
)

# p値でフィルタリング
significant_markers <- all_markers[all_markers$p_val_adj < p_value_threshold, ]

# 結果をCSVファイルとして保存
write.csv(significant_markers, file = "Differentially_Expressed_Genes.csv", row.names = FALSE)

# 上位10個のDEGを可視化
top10 <- significant_markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

pdf(file = "Top10_Differentially_Expressed_Genes.pdf", width = 8.27, height = 11.69)
DoHeatmap(combined, features = top10$gene) + NoLegend()
dev.off()