###### 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")
```
<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
```
```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
```
<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)
```
<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")
```
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)
```
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)
```
Sign in with Wallet
Connect another wallet