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

# 22-04： ES由来セルトリ細胞scRNAseq解析　FTSLC_scRNAseq 第２弾 (aggrしない) gonadal somatic cellsを抽出

<br>

#### 2026/05/16
#### 【 作業場所＆保存場所 】

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

---

#### ※ aggrしないでSeuratした方

【参考】
scRNA-seqデータの解析(応用編)
https://hackmd.io/@takahirosuzuki/ByoEBNVPL
<br>

#### 血管内皮細胞、血液細胞、生殖細胞などを取り除いてgonadal somatic cellsを抽出
1. 再クラスタリング
2. マーカー発現check(violin plot)(大まかな細胞種分け）
3. 全データでクラスター間のcorrelation matirx
4. XY_E12.5とXY_FTSLCsで各クラスターの割合(conposition)を解析
5. pseudotime analysis (DDRTree)

<br>

### 1. gonadal somatic cellsを抽出し再クラスタリング

関数の定義

#--------------------------------------------
# subsetUmapClust
#--------------------------------------------
#' Extract subset and do clustering
#' 
#' Extract the subsets from a seurat object and do scaling, and PCA, UMAP and tSNE dimensionality reductions
#' 
#' @param data A seurat object
#' @param subset.name Parameter to subset on. Eg, the name of a gene, PC_1, a column name in [e-mail], etc. Any argument that can be retreived using FetchData
#' @param subset.value Returns cells with the subset name equal to this value
#' @param resolution resolution for clustering
#' 
#' @importFrom Seurat Read10X CreateSeuratObject RenameCells NormalizeData FindVariableFeatures
#' 
#' @return a seurat object
#' 
#' @keywords seurat 10x
#' @export 
#' 
subsetUmapClust <- function(data,  subset.name = "seurat_clusters", subset.value, npcs = 30, dims = 1:30, resolution = 0.6){
    if(length(subset.value) == 1){
        data.subset <- eval(parse(text=paste0("subset(data, subset = ", subset.name, " == ", subset.value, ")" )))
    }else{
        multisub <- paste0("subset(data, ", 
                                paste(paste(subset.name, paste0("'",subset.value,"'"), sep="=="), collapse = "|"),
                                ")"
        )
        data.subset <- eval(parse(text=multisub))
    }
    data.subset <- RunPCA(object = data.subset, npcs = npcs, verbose = FALSE)
    data.subset <- RunUMAP(object = data.subset, reduction = "pca", dims = dims)
    data.subset <- FindNeighbors(object = data.subset, dims = dims)
    data.subset <- FindClusters(object = data.subset, reduction.type = "umap", resolution =  resolution)
    return(data.subset)
}


## ここで保存
save(list=ls(), file="8_20220516.RData")

gonadal somatic cellsを抽出 （resolution=0.3におちついた）

CL012346_reso03_DimPlot (8x10インチ)

load("8_20220516.RData")

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


data.integrated.gsc <- subsetUmapClust(
    data = data.integrated,
    subset.name = "seurat_clusters",
    subset.value = c("0", "1", "2", "3", "4", "6"),
    npcs = 14,
    dims = 1:10,
    resolution=0.3
)

DimPlot(object = data.integrated.gsc, label = TRUE)

![](https://i.imgur.com/JyKEdU7.jpg)

DimPlot(object = data.integrated.gsc, reduction = "umap" , label = TRUE , split.by = "sample") + scale_color_ucscgb()



<br>

### 2. マーカー発現check(violin plot)(大まかな細胞種分け）

progenitor = c("Sox11", "Ecm1", "Nr2f1")
granulosa = c("Kitl", "Inha", "Foxl2", "Runx1","Fst", "Akr1cl", "Cdkn1b", "Aard", "Hmgcs2")
stromal = c("Wnt5a", "Pdgfra", "Tcf21", "Acta2", "Arx", "Gng13", "Lgr5", "Apoc1", "Gpc3", "Hmcn1")
set1 = c("Nr5a1", "Arx", "Sfrp1", "Nr2f2")
set2 = c("Cdkn1b", "Fst", "Amhr2", "Mt1")
sertoli = c("Sox9", "Mro", "Aard", "Amh", "Dhh", "Ptgds", "Hsd17b3")
pre_sertoli = c("Runx1", "Sry", "Nr0b1")
fetal_leydig = c("Cyp11a1", "Cyp17a1", "Star", "Insl3", "Hsd3b1")

#MultiVlnPlot_gsc_res03_" (5x10)
multiVlnPlot(object=data.integrated.gsc, features = progenitor, title="progenitor markers , resolution=0.3")
multiVlnPlot(object=data.integrated.gsc, features = granulosa, title="granulosa markers , resolution=0.3")
multiVlnPlot(object=data.integrated.gsc, features = stromal, title="stromal, markers , resolution=0.3")
multiVlnPlot(object=data.integrated.gsc, features = set1, title="set1 markers , resolution=0.3")
multiVlnPlot(object=data.integrated.gsc, features = set2, title="set2 markers , resolution=0.3")
multiVlnPlot(object=data.integrated.gsc, features = sertoli, title="sertoli markers , resolution=0.3")
multiVlnPlot(object=data.integrated.gsc, features = pre_sertoli, title="pre_sertoli markers , resolution=0.3")
multiVlnPlot(object=data.integrated.gsc, features = fetal_leydig, title="fetal_leydig, markers , resolution=0.3")

save(list=ls(), file="9_gsc_res03_20220516.RData")


<br>

### 3. 全データでクラスター間のcorrelation matirx

load("5_20220516_cl02.RData")

library(Seurat)
library(dplyr)
library(patchwork)
library(MAST)
library(ggsci)
library("ComplexHeatmap")
library("circlize")
library("RColorBrewer")
library(scales)

av.exp <- AverageExpression(data.integrated.gsc)$integrated
cor.exp <- cor(av.exp)

#類似度に基づいたクラスタリング
hc <- hclust(dist(cor.exp))
breaks1 <- seq(0.95, 1, length=500)
mycols1 <- colorRamp2(breaks1,colorRampPalette(brewer.pal(9, "GnBu"))(500))

ht1 <- Heatmap(cor.exp,
                col = mycols1,
                name = "Correlation Coefficient", #title of legend, 
                row_names_gp = gpar(fontsize = 12),   
                column_names_gp = gpar(fontsize = 12), 
                column_names_rot = 0,
                row_dend_width = unit(2, "cm"),   
                column_dend_height = unit(2, "cm"), 
                row_dend_side = "right",
                column_dend_side = "bottom",
                cluster_rows = hc,
                cluster_columns = hc,
                cell_fun = function(j, i, x, y, width, height, fill){
                        grid.rect(x = x, y = y, width = width, height = height, gp = gpar(col = "black", fill = NA))
                }
            )

ht1


<br>

### 4. XY_E12.5とXY_FTSLCsで各クラスターの割合(conposition)を解析

#割合の計算
prop_cluster <- prop.table(x=table([e-mail], [e-mail]$sample), margin =2)

#plot用にデータフレームに変換
Cluster_name=factor(rownames(prop_cluster), levels=rownames(prop_cluster))    
df <- data.frame(Cluster=Cluster_name,XY_12.5=prop_cluster[,7], XY_FTSLCs=prop_cluster[,8])
df2 <- tidyr::gather(df, key=Sample, value=Proportion, -Cluster, factor_key = TRUE)

#ggplotでvidualization
theme <- theme(panel.background = element_blank(),    # initialization
               panel.border = element_blank(),
               panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               strip.background = element_blank(),
               plot.title = element_text(size = 10,hjust = 0.5),
               axis.text.x = element_text(size=10, angle=90),
               axis.text.y = element_text(size=10),
               axis.title.x = element_blank(),
               axis.title.y = element_text(size=10),
               axis.ticks = element_line(size=0.5),
               axis.line = element_line(size=0.5),
               plot.margin = unit(c(1,1,1,1),"line")
)
 
g <- ggplot(df2, aes(x= Sample, y = Proportion, fill=Cluster))
g <- g + geom_bar(stat = "identity")
g <- g + scale_y_continuous(labels = percent)    # display as %
g <- g + theme

g


<br>

### 5. pseudotime analysis (DDRTree)

library(monocle3)
library(SeuratWrappers)

# Converting to monocle3
cds <- as.cell_data_set(data.integrated.gsc, assay = "RNA")

# monocle3  clustering
cds <- cluster_cells(cds = cds, verbose = TRUE)

# monocle3 learning graph 
cds <- learn_graph(cds, use_partition = FALSE , close_loop = FALSE)

#Plot
p <- plot_cells(cds,
	label_groups_by_cluster = FALSE,
	group_label_size = 5,
	cell_size = 0.5,
	show_trajectory_graph = TRUE,
	color_cells_by = "ident",
	label_leaves = TRUE,
	label_branch_points = TRUE,
	label_cell_groups=FALSE,
	trajectory_graph_segment_size = 1.5)

p