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

ちょっとごちゃごちゃしてきたので一旦整理すると

#Seuartパッケージのインストール
install.packages('Seurat')
library(Seurat)

#テストファイルの準備とdirectoryの移動
setwd("~/scrnaseqtest/immune_alignment_matrix/")

#data読み込み
ctrl.data <- read.table("immune_control_expression_matrix.txt", sep = '\t', row.names = 1)
stim.data <- read.table("immune_stimulated_expression_matrix.txt",  sep = "\t", row.names = 1)

#Seuratオブジェクトの作成
ctrl <- CreateSeuratObject(counts = ctrl.data, project = "IMMUNE_CTRL", min.cells = 5)
stim <- CreateSeuratObject(counts = stim.data, project = "IMMUNE_STIM", min.cells = 5)

#metadataの付与
ctrl@meta.data$ctrl <- "CTRL"
stim@meta.data$stim <- "STIM"

#フィルタリング
ctrl <- subset(ctrl, subset = nFeature_RNA > 500)
stim <- subset(stim, subset = nFeature_RNA > 500)

#正規化
ctrl <- NormalizeData(ctrl, scale.factor = 10000)
stim <- NormalizeData(stim, scale.factor = 10000)

#Z-Scoreの算出
ctrl <- ScaleData(ctrl)
stim <- ScaleData(stim)

#使用する遺伝子の選定
ctrl <- FindVariableFeatures(ctrl)
stim <- FindVariableFeatures(stim)

# 上位1000の変動遺伝子を取得
g.1 <- head(VariableFeatures(ctrl), 1000)
g.2 <- head(VariableFeatures(stim), 1000)

#g1、g2で共通して発現変動の大きい遺伝子をunique関数で抽出
genes.use <- unique(c(g.1, g.2))
scaled_genes <- rownames(GetAssayData(ctrl, slot = "scale.data"))
genes.use <- intersect(genes.use, scaled_genes)
scaled_genes <- rownames(GetAssayData(stim, slot = "scale.data"))
genes.use <- intersect(genes.use, scaled_genes)

とりあえずここまでは良さげ。

データを可視化していく。

CCA（Canonical Correlation Analysis）によるオブジェクトのマージ

# データセットをリストにまとめる
immune.list <- list(ctrl = ctrl, stim = stim)

# アンカーを見つける
anchors <- FindIntegrationAnchors(object.list = immune.list, anchor.features = genes.use)

# データセットの統合
immune.combined <- IntegrateData(anchorset = anchors, dims = 1:30)

# 統合データのスケーリングと主成分分析
immune.combined <- ScaleData(immune.combined)
immune.combined <- RunPCA(immune.combined)

PCAのプロットを出してみる

# 次元削減結果のプロット
p1 <- DimPlot(object = immune.combined, reduction = "pca", group.by = "stim", pt.size = 0.5)
p1

PC_1をy軸に取ったViolin plot

p2 <- VlnPlot(object = immune.combined, features = "PC_1", group.by = "stim")
p2

各次元の重要な遺伝子のヒートマップを表示

DimHeatmap(object = immune.combined, dims = 1:2, cells = 500, balanced = TRUE)

とりあえずできた気になっているが、自分のデータで意味のある解析ができるかというとかなり怪しい。