---
# System prepended metadata

title: 22-05： ES由来セルトリ細胞scRNAseq解析　FTSLC_scRNAseq 第２弾 (aggrしない) Sertoli 発現クラスターを抽出
tags: [解析, ノート, TS, single-cell RNA-seq]

---

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

# 22-05： ES由来セルトリ細胞scRNAseq解析　FTSLC_scRNAseq 第２弾 (aggrしない) Sertoli 発現クラスターを抽出

<br>

#### 2022/05/30
#### 【 作業場所＆保存場所 】

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


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

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

#### Sertoli markerが発現しているクラスター(Cluster1と6)を抽出
1. 再クラスタリング
2. granulosa, Sertoli cell markerのviolin plot
3. 結果を見ていくつかの遺伝子はfeature plot
4. Mki67(増殖マーカー）の発現をチェック（feature plot)

<br>

### 1. Sertoli markerが発現しているCluster(1と6)を抽出し再クラスタリング

```r=
#クラスターの抽出・再解析
#関数の定義
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=c("data.integrated.gsc","subsetUmapClust"), file="10_dataintegrated_gsc.RData")
```
resolutionを0.1〜0.3でふってみる。
```r=
#resolution=0.1
data.integrated.gsc16 <- subsetUmapClust(
    data = data.integrated.gsc,
    subset.name = "seurat_clusters",
    subset.value = c("1", "6"),
    npcs = 14,
    dims = 1:10,
    resolution=0.1)
#20220527_data.integrated.gsc16_DimPlot01
DimPlot(object = data.integrated.gsc16, label = TRUE)

#resolution=0.2
data.integrated.gsc16_02 <- subsetUmapClust(
    data = data.integrated.gsc,
    subset.name = "seurat_clusters",
    subset.value = c("1", "6"),
    npcs = 14,
    dims = 1:10,
    resolution=0.2)
#20220527_data.integrated.gsc16_DimPlot02
DimPlot(object = data.integrated.gsc16_02, label = TRUE)

#resolution=0.3
data.integrated.gsc16_03 <- subsetUmapClust(
    data = data.integrated.gsc,
    subset.name = "seurat_clusters",
    subset.value = c("1", "6"),
    npcs = 14,
    dims = 1:10,
    resolution=0.3)
#20220527_data.integrated.gsc16_DimPlot03
DimPlot(object = data.integrated.gsc16_03, label = TRUE)
```
plotは後でまとめて


### 2. granulosa, sertoli cell markerのviolin plot

```r=
granulosa = c("Kitl", "Inha", "Foxl2", "Runx1","Fst", "Akr1cl", "Cdkn1b", "Aard", "Hmgcs2")
sertoli = c("Sox9", "Mro", "Aard", "Amh", "Dhh", "Ptgds", "Hsd17b3")

#20220527_multiVlnPlot_gsc1601_
multiVlnPlot(object=data.integrated.gsc16, features = granulosa, title="granulosa markers , resolution=0.1")
multiVlnPlot(object=data.integrated.gsc16, features = sertoli, title="sertoli markers , resolution=0.1")

#20220527_multiVlnPlot_gsc1602_
multiVlnPlot(object=data.integrated.gsc16_02, features = granulosa, title="granulosa markers , resolution=0.2")
multiVlnPlot(object=data.integrated.gsc16_02, features = sertoli, title="sertoli markers , resolution=0.2")

#20220527_multiVlnPlot_gsc1603_
multiVlnPlot(object=data.integrated.gsc16_03, features = granulosa, title="granulosa markers , resolution=0.3")
multiVlnPlot(object=data.integrated.gsc16_03, features = sertoli, title="sertoli markers , resolution=0.3")

#20220527_multiVlnPlot_gsc1604_
multiVlnPlot(object=data.integrated.gsc16_04, features = granulosa, title="granulosa markers , resolution=0.4")
multiVlnPlot(object=data.integrated.gsc16_04, features = sertoli, title="sertoli markers , resolution=0.4")
```

![](https://i.imgur.com/O2JAPzW.png)


<br>

### 3. granulesaとsertoliの各マーカーのFeaturePlot
```r=
# 20220527_FeaturePlot_Kitl (5 x 35)
FeaturePlot(object=data.integrated.gsc16, features= "Kitl", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Inha
FeaturePlot(object=data.integrated.gsc16, features= "Inha", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Foxl2
FeaturePlot(object=data.integrated.gsc16, features= "Foxl2", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Runx1
FeaturePlot(object=data.integrated.gsc16, features= "Runx1", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Fst
FeaturePlot(object=data.integrated.gsc16, features= "Fst", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Akr1cl
FeaturePlot(object=data.integrated.gsc16, features= "Akr1cl", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Cdkn1b
FeaturePlot(object=data.integrated.gsc16, features= "Cdkn1b", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Aard
FeaturePlot(object=data.integrated.gsc16, features= "Aard", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Hmgcs2
FeaturePlot(object=data.integrated.gsc16, features= "Hmgcs2", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
```
```
pdfjoin --outfile 20220527_FeaturePlot_01_granulosa.pdf 20220527_FeaturePlot_01_Kitl.pdf 20220527_FeaturePlot_01_Inha.pdf 20220527_FeaturePlot_01_Foxl2.pdf 20220527_FeaturePlot_01_Runx1.pdf 20220527_FeaturePlot_01_Fst.pdf 20220527_FeaturePlot_01_Akr1cl.pdf 20220527_FeaturePlot_01_Cdkn1b.pdf 20220527_FeaturePlot_01_Aard.pdf 20220527_FeaturePlot_01_Hmgcs2.pdf
```
![](https://i.imgur.com/zrCWAH8.jpg)

```r=
# 20220527_FeaturePlot_Sox9 (5 x 35)
FeaturePlot(object=data.integrated.gsc16, features= "Sox9", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Mro
FeaturePlot(object=data.integrated.gsc16, features= "Inha", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Aard
FeaturePlot(object=data.integrated.gsc16, features= "Aard", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Amh
FeaturePlot(object=data.integrated.gsc16, features= "Amh", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Dhh
FeaturePlot(object=data.integrated.gsc16, features= "Dhh", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Ptgds
FeaturePlot(object=data.integrated.gsc16, features= "Ptgds", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
# 20220527_FeaturePlot_Hsd17b3
FeaturePlot(object=data.integrated.gsc16, features= "Hsd17b3", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
```
```
pdfjoin --outfile 20220527_FeaturePlot_sertoli.pdf 20220527_FeaturePlot_Sox9.pdf 20220527_FeaturePlot_Mro.pdf 20220527_FeaturePlot_Aard.pdf 20220527_FeaturePlot_Amh.pdf 20220527_FeaturePlot_Dhh.pdf 20220527_FeaturePlot_Ptgds.pdf 20220527_FeaturePlot_Hsd17b3.pdf
```
![](https://i.imgur.com/Y4QDUzM.png)

<br>

### 4. Mki67(増殖マーカー）とSryのFeaturePlot
```r=
# FeaturePlot_Mki67
FeaturePlot(object=data.integrated.gsc16, features= "Mki67", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)

# FeaturePlot_Sry
FeaturePlot(object=data.integrated.gsc16, features= "Sry",
            col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
```
![](https://i.imgur.com/1TzbYL7.png)
![](https://i.imgur.com/ftHCEi8.png)

<br>

### 5. vivo vs vitro (FTSLCs)のDEG解析

gonadal somatic cellsを抽出データからAmhまたはDhhまたはHsd3b1が発現している細胞（多分0以上というcut-offでよいと思います）を抽出してきて、
XY_11.5 vs XY_FTSLCs, XY12.5 vs XY_FTSLCsのDEG解析。

```r=
# gonadal somatic cellsを抽出データからAmhまたはDhhまたはHsd3b1が発現している細胞を抽出
gsc_ADH <- subset(x = data.integrated.gsc, subset = Amh > 0 | Dhh > 0 | Hsd3b1 > 0)

# Idents設定
Idents(object =gsc_ADH) <- gsc_ADH@meta.data$sample

## find DEGs  XY115 vs FTSLCs
diff_genes_XY115_FTSLCs <- FindMarkers(
object = gsc_ADH,
ident.1 = "XY_11.5",
ident.2 = "XY_FTSLCs",
min.pct = 0.25,
logfc.threshold = 0.25,
test.use = "MAST")
# export degs
diff_genes_XY115_FTSLCs %>%
dplyr::filter(avg_log2FC >= 0.5) %>%
dplyr::filter(p_val_adj < 0.01) -> degs

## find DEGs  XY125 vs FTSLCs
diff_genes_XY125_FTSLCs <- FindMarkers(
object = gsc_ADH,
ident.1 = "XY_12.5",
ident.2 = "XY_FTSLCs",
min.pct = 0.25,
logfc.threshold = 0.25,
test.use = "MAST")
# export degs2
diff_genes_XY125_FTSLCs %>%
dplyr::filter(avg_log2FC >= 0.5) %>%
dplyr::filter(p_val_adj < 0.01) -> degs2
```
XY_11.5 vs XY_FTSLCs 結果
```
> dim(degs)
[1] 28  5
                  p_val avg_log2FC pct.1 pct.2     p_val_adj
Cpa2      4.498229e-165  0.5884531 0.750 0.977 1.396835e-160
Serpina3a  1.883720e-95  0.5910985 0.940 0.998  5.849517e-91
Bmyc       6.305971e-91  1.0380883 0.982 0.934  1.958193e-86
Kctd14     6.225445e-89  1.2114003 0.937 0.946  1.933187e-84
Serpinb6b  1.053462e-86  0.5822037 0.834 0.951  3.271315e-82
Ifitm1     1.515101e-82  1.1116480 0.999 0.999  4.704842e-78
Wnt6       1.470495e-77  0.7512319 0.939 0.976  4.566330e-73
Irx3       1.424943e-74  0.6015774 0.991 0.973  4.424875e-70
Rgs2       2.445254e-67  0.6317426 0.939 0.960  7.593248e-63
Lgals7     4.915839e-66  0.9254498 0.999 0.999  1.526516e-61
Igfbp2     9.303883e-66  0.8199994 0.990 0.995  2.889135e-61
Ccdc61     2.455491e-58  0.6067288 0.960 0.971  7.625037e-54
Lhx9       5.623103e-58  0.5807709 0.957 0.935  1.746142e-53
Cst8       6.526003e-57  0.5350147 0.869 0.911  2.026520e-52
Bst2       1.923476e-56  0.5968330 0.840 0.915  5.972969e-52
Timm8a1    3.298512e-56  0.5684563 0.983 0.996  1.024287e-51
Actr6      2.719391e-53  0.5330789 0.946 0.974  8.444523e-49
Fst        1.213634e-50  0.8695522 0.969 0.912  3.768697e-46
Amhr2      1.173178e-48  0.5567155 0.967 0.954  3.643070e-44
Dbi        1.481149e-46  0.5724643 0.983 0.992  4.599411e-42
Hmgcs2     4.630715e-44  0.5825996 0.774 0.553  1.437976e-39
Serpinb6a  1.063653e-42  0.5675884 0.927 0.968  3.302963e-38
Akr1cl     1.423428e-41  0.7311550 0.868 0.877  4.420172e-37
Tkt        1.064314e-39  0.5931712 0.969 0.977  3.305013e-35
Gadd45g    2.059316e-39  0.6166051 0.828 0.674  6.394793e-35
Laptm4b    1.989857e-38  0.5216030 0.970 0.995  6.179103e-34
Pglyrp1    4.313674e-31  0.5127387 0.819 0.867  1.339525e-26
Cdkn1c     9.838466e-20  0.5469654 0.997 0.963  3.055139e-15
```

XY12.5 vs XY_FTSLCs結果
```
> dim(degs2)
[1] 15  5
                p_val avg_log2FC pct.1 pct.2     p_val_adj
Serpina3a 1.566152e-171  0.6145488 0.837 0.998 4.863372e-167
Bst2      3.111447e-158  0.5352388 0.752 0.915 9.661975e-154
Cpa1      8.607218e-128  0.5289280 0.715 0.950 2.672799e-123
Wnt6      2.806593e-104  0.7318353 0.950 0.976 8.715314e-100
Akr1cl    2.789744e-102  0.7546402 0.803 0.877  8.662991e-98
Rgs2      3.224254e-102  0.5161843 0.894 0.960  1.001228e-97
Cpa2       2.979597e-98  0.5273586 0.885 0.977  9.252544e-94
Lgals7     1.082519e-96  0.9862038 0.999 0.999  3.361547e-92
Serpinb6a  2.130182e-87  0.5960677 0.910 0.968  6.614856e-83
Kctd14     3.703850e-84  1.0500462 0.980 0.946  1.150157e-79
Bmyc       1.424592e-81  0.8331254 0.978 0.934  4.423786e-77
Irx3       1.277675e-75  0.6233372 0.993 0.973  3.967564e-71
Ifitm1     2.965591e-70  0.9713347 1.000 0.999  9.209049e-66
Fst        5.470956e-52  0.8063158 0.969 0.912  1.698896e-47
Cdkn1c     2.474446e-29  0.5622118 0.999 0.963  7.683899e-25
```

発現が下がったのも。
dplyr::filter(avg_log2FC >= 0.5)を dplyr::filter(avg_log2FC <= -0.5) で。

```r=
diff_genes_XY115_FTSLCs %>%
    dplyr::filter(avg_log2FC <= -0.5) %>%
    dplyr::filter(p_val_adj < 0.01) -> Ddegs

diff_genes_XY125_FTSLCs %>%
    dplyr::filter(avg_log2FC <= -0.5) %>%
    dplyr::filter(p_val_adj < 0.01) -> Ddegs2
```
```
> dim(Ddegs)
[1] 25  5
> 
> dim(Ddegs2)
[1] 32  5
```

<br>

### 6. その他、追加

gonadalsomaticcells抽出-> 1,6クラスター抽出0.2 の サンプルごとのUMAP
```r=
# サンプルごとのUMAP
DimPlot(object = data.integrated.gsc16_02 , reduction = "umap" , label = TRUE , split.by = "sample")
```
![](https://i.imgur.com/EUW8xPT.png)

gonadalsomaticcells抽出-> 1,6クラスター抽出0.2 の HeatMap
```r=
# クラスター特異的遺伝の同定
integrated.markers <- FindAllMarkers(
    object = data.integrated.gsc16_02,
    only.pos = FALSE,
    min.pct = 0.25,
    logfc.threshold = 0.25,
    test.use = "MAST")

# 計算に少し時間がかかるのでこのオブジェクトはsaveしておく
save(data.integrated.gsc16_02, file="11_data.integrated.gsc16_02.RData")

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

#20220601_HeatMap_Top5.pdf
DoHeatmap(object = data.integrated.gsc16_02, features=top5.posi$gene, size=3)
```
![](https://i.imgur.com/KfpeRgk.png)

gonadalsomaticcells抽出-> 1,6クラスター抽出0.2 のPax8のFeaturePlot
```r=
# 20220527_FeaturePlot_Pax8
FeaturePlot(object=data.integrated.gsc16_02, features= "Pax8", col=c("gray", "red"), split.by="sample", min.cutoff=0.1, pt.size=1)
```
![](https://i.imgur.com/w76Oj2A.png)
