###### tags: `ノート` `TS` `解析` `single-cell RNA-seq`

# 22-02： ES由来セルトリ細胞scRNAseq解析　FTSLC_scRNAseq (aggrしない)

<br>

#### 2026/04/06
#### 【 作業場所＆保存場所 】

• 作業場所
　~/osc-fs/22_TS_FTSLC_scRNAseq_2022/Seurat_Non_aggr

---

### aggrしないでSeuratする

【参考】
Seuratによる10Xデータの解析(基本編)
https://hackmd.io/@takahirosuzuki/B1Tmsn5LI
<br>

#### 1. 関数の定義

library(Seurat)
library(dplyr)
library(patchwork)
library(MAST)

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)
}

<br>

#### 2. scCountProcを使用してデータの読み込み

data.E10_5_XX <- scCountProc(
	data_path = "//osc-fs_home/s-maeda/22_TS_FTSLC_scRNAseq_2022/E10-5_XX",　
	sample_name = "E10_5_XX",　　　
	normalization = FALSE　　　
)
data.E11_5_XX <- scCountProc(
	data_path = "//osc-fs_home/s-maeda/22_TS_FTSLC_scRNAseq_2022/E11-5_XX",
	sample_name = "E11_5_XX",
	normalization = FALSE
)
data.E12_5_XX <- scCountProc(
	data_path = "//osc-fs_home/s-maeda/22_TS_FTSLC_scRNAseq_2022/E12-5_XX",
	sample_name = "E12_5_XX",
	normalization = FALSE
)
data.FOSLCsXX <- scCountProc(
	data_path = "//osc-fs_home/s-maeda/22_TS_FTSLC_scRNAseq_2022/Induced_D6",　
	sample_name = "FOSLCs",　　　
	normalization = FALSE　　　
)

data.E10_5_XY <- scCountProc(
	data_path = "//osc-fs_home/s-maeda/22_TS_FTSLC_scRNAseq_2022/E10-5_XY",　
	sample_name = "E10_5_XY",　　　
	normalization = FALSE　　　
)
data.E11_5_XY <- scCountProc(
	data_path = "//osc-fs_home/s-maeda/22_TS_FTSLC_scRNAseq_2022/E11-5_XY",
	sample_name = "E11_5_XY",
	normalization = FALSE
)
data.E12_5_XY <- scCountProc(
	data_path = "//osc-fs_home/s-maeda/22_TS_FTSLC_scRNAseq_2022/E12-5_XY",
	sample_name = "E12_5_XY",
	normalization = FALSE
)
data.FTSLCsXY <- scCountProc(
	data_path = "//osc-fs_home/s-maeda/22_TS_FTSLC_scRNAseq_2022/10x-ig168",　
	sample_name = "FTSLCsXY",　　　
	normalization = FALSE　　　
)


#データをリスト化(アルファベットの順番注意!)
data.list <- list(XX_E10_5=data.E10_5_XX, XX_E11_5=data.E11_5_XX, XX_E12_5=data.E12_5_XX, XX_FOSLCs=data.FOSLCsXX,  XY_E10_5=data.E10_5_XY, XY_E11_5=data.E11_5_XY, XY_E12_5=data.E12_5_XY, XY_FTSLCs=data.FTSLCsXY)

#個別のオブジェクトは以降必要ないので消去（メモリ節約のため）
rm(list=c("data.E10_5_XY", "data.E11_5_XY", "data.E12_5_XY", "data.E10_5_XX", "data.E11_5_XX", "data.E12_5_XX",data.FOSLCsXX,data.FTSLCsXY))

<br>

#### 3. QC

vp_list <- list()　　　
combp_list <- list()　　　

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))))
}

vp_list
combp_list


20220415_vln.pdf と 20220415_cmb.pdf

<br>

#### 4. violin plotとscatter plotをみてcut-offを決めて範囲内の細胞を抽出

data.list[["XX_E10_5"]] <- subset(data.list[["XX_E10_5"]], subset = nFeature_RNA > 1000 & nFeature_RNA < 4000 & percent.mt < 12)
data.list[["XX_E11_5"]] <- subset(data.list[["XX_E11_5"]], subset = nFeature_RNA > 2000 & nFeature_RNA < 6000 & percent.mt < 10)
data.list[["XX_E12_5"]] <- subset(data.list[["XX_E12_5"]], subset = nFeature_RNA > 2000 & nFeature_RNA < 5000 & percent.mt < 10 )
data.list[["XX_E10_5"]] <- subset(data.list[["XX_E10_5"]], subset = nFeature_RNA > 1500 & nFeature_RNA < 4000 & percent.mt < 12)
data.list[["XY_E10_5"]] <- subset(data.list[["XY_E10_5"]], subset = nFeature_RNA > 1000 & nFeature_RNA < 4000 & percent.mt < 10 )
data.list[["XY_E11_5"]] <- subset(data.list[["XY_E11_5"]], subset = nFeature_RNA > 500 & nFeature_RNA < 5000 & percent.mt < 20)
data.list[["XY_E12_5"]] <- subset(data.list[["XY_E12_5"]], subset = nFeature_RNA > 2000 & nFeature_RNA < 5000 & percent.mt < 8 )
data.list[["XY_E10_5"]] <- subset(data.list[["XY_E10_5"]], subset = nFeature_RNA > 1000 & nFeature_RNA < 4000 & percent.mt < 12 )

<br>

#### 5. Normalizationと変動遺伝子の検出

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)
}

<br>

#### 6. Seurat objectを統合

data.anchors <- FindIntegrationAnchors(object.list = data.list, dims = 1:30)                                         #######ここ20220415、15：00
all_genes <- [e-mail][[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(list=ls(), file="1_20220421.RData")

<br>

#### 7. Cell-Cycleの影響を調整

cell cycle markerリスト(Tirosh et al, 2015)の読み込み

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に偏りがあるため

#データ全体をscaling
data.integrated <- ScaleData(object =data.integrated, features = rownames(data.integrated), verbose = FALSE)

#cell-cycle スコアリング
#G1, G2/M, Sをそれぞれ特異的遺伝子に基づいてスコアリングする
data.integrated <- CellCycleScoring(data.integrated, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# PC
data.integrated <- RunPCA(data.integrated, features = c(s.genes, g2m.genes))
DimPlot(data.integrated) 

save(list=ls(), file="2_20220422.RData")

#### Cell-Cycleの影響の調整（時間かかる）
G1, G2/M, Sすべてを調整　

data.integrated <- ScaleData(data.integrated, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(data.integrated))

#Cell-Cycleの調整の確認
data.integrated  <- RunPCA(data.integrated, features = VariableFeatures(data.integrated), nfeatures.print = 10)
data.integrated  <- RunPCA(data.integrated, features = c(s.genes, g2m.genes))

pdf("20220422_after_cellcycle.pdf")
DimPlot(data.integrated) 
dev.off()

save(list=ls(), file="3_20220423.RData")


<br>

#### 8.　 クラスタリングとUMAPでの次元圧縮・可視化

#デフォルトをintegrated assayに変更
DefaultAssay(object = data.integrated) <- "integrated"

#Dimensionalityの決定
data.integrated <- JackStraw(data.integrated, num.replicate = 100, dims = 50)
data.integrated <- ScoreJackStraw(data.integrated, dims = 1:50)

pdf("20220422_JackStrawPlot.pdf")
JackStrawPlot(data.integrated, dims = 1:50)
dev.off()

pdf("20220422_ElbowPlot.pdf")
ElbowPlot(data.integrated,  ndims = 50)
dev.off()

## 一度ここで保存
save(data.integrated, file="4_20220423.RData")


<br>

---
2022/4/26

クラスタリング(dims=1:10、resolution = 0.2 に決定)

#クラスタリング
data.integrated <- RunPCA(object = data.integrated, npcs = 20, verbose = FALSE)
data.integrated <- RunUMAP(object = data.integrated, reduction = "pca", dims = 1:10)   #ここのdimsをJackStrawPlotかElbowPlotの結果で調整
data.integrated <- FindNeighbors(object = data.integrated, dims =  1:10)    #ここのdimsをJackStrawPlotかElbowPlotの結果で調整
data.integrated <- FindClusters(object = data.integrated, algorithm = 1, resolution = 0.2) #resolutionを0.2にする

## 一度ここで保存
save(data.integrated, file="5_20220426_cl02.RData")

<br>

#### 9. meta.dataの追加

サンプル名の変更

data.integrated[["sample"]][data.integrated[["sample"]] == "E10_5_XX"] <- "XX_10.5"
data.integrated[["sample"]][data.integrated[["sample"]] == "E11_5_XX"] <- "XX_11.5"
data.integrated[["sample"]][data.integrated[["sample"]] == "E12_5_XX"] <- "XX_12.5"
data.integrated[["sample"]][data.integrated[["sample"]] == "FOSLCs"] <- "XX_FOSLCs"
data.integrated[["sample"]][data.integrated[["sample"]] == "E10_5_XY"] <- "XY_10.5"
data.integrated[["sample"]][data.integrated[["sample"]] == "E11_5_XY"] <- "XY_11.5"
data.integrated[["sample"]][data.integrated[["sample"]] == "E12_5_XY"] <- "XY_12.5"
data.integrated[["sample"]][data.integrated[["sample"]] == "FTSLCsXY"] <- "XY_FTSLCs"

#factorを設定することで可視化の際の並び順を指定できる。
levels(data.integrated[["sample"]]) <- c("XX_10.5", "XX_11.5", "XX_12.5", "XX_FOSLCs", "XY_10.5", "XY_11.5", "XY_12.5", "XY_FTSLCs")　

## 一度ここで保存
save(data.integrated, file="5_20220516_cl02.RData")

#umapで可視化
#20220426_UMAP.pdf
DimPlot(object = data.integrated , reduction = "umap" , label = TRUE , split.by = "sample")


<br>

#### 10. クラスター特異的遺伝の同定

integrated.markers <- FindAllMarkers(
    object = data.integrated,
    only.pos = FALSE,
    min.pct = 0.25,
    logfc.threshold = 0.25,
    test.use = "MAST"
)

## ここで保存
save(data.integrated, file="6_20220426_MAST.RData")

<br>

#### 11. 特異的発現遺伝子Top5のヒートマップ

integrated.markers %>%
    filter(avg_log2FC >= 0) %>%
    group_by(cluster) %>%
    top_n(5, avg_log2FC) -> top5.posi

#20220426_HeatMap_Top5.pdf
DoHeatmap(object = data.integrated, features=top5.posi$gene, size=3)

save(data.integrated, file="7_20220426_data.integrated.RData")