---
title: "Comparison of WT vs KO AOM_DSS mouse model"
author: Anne Hartebrodt
date: "09.01.2024"
output:
html_document:
number_sections: true
---ZEB1 differential abundance analysis of myCAF and iCAF like phenotypes.
The goal of this analysis is to elucidate if there is a change in cell type composition between the two conditions. We will use integrated clustering and differential abundance analysis.
knitr::opts_chunk$set(message = FALSE, warning = FALSE, echo=TRUE, include=TRUE)require(scater)
require(scran)
require(Seurat)
require(stringr)
require(patchwork)
require(DAseq)
require(harmony)
require(data.table)
require(dplyr)
require(decoupleR)
library(dorothea)
library(dplyr)
library(tibble)
library(tidyr)
library(patchwork)
library(ggplot2)
library(pheatmap)Reproduction setup
Input files
In order to use this notebook please note that you should create an
analysis folder with a subdirectory data. In the course of
the notebook execution a results folder with various subdirectories will
be created automatically and an html version of the notebook will be
created in the directory where you execute this notebook.
## Your fully qualified base directory.
analysis_dir<-'~/Documents/external_analyses/marc_stemmler/'
data_dir<-file.path(analysis_dir, 'data')
input_data<-file.path(data_dir, 'merged.rescaled2.RData')
marker_file<-file.path(data_dir, 'marker_genes.tsv.csv')Python env
DASeq relies on python. Make sure to update the python environment accordingy.
python_env<-'~/anaconda3/envs/daseq2/bin/python3'Output settings
results.dir<-file.path(analysis_dir, 'results')
dir.create(results.dir)Data preparation
Load the data and set seed
sce.data<-get(load(input_data))
set.seed(12)Make a seurat object and make the gene names usable
sce.seurat <-as.Seurat(sce.data, data = NULL)
genes<-rownames(sce.seurat@assays$originalexp)
genes<-as.character(sapply(genes, function(x) str_split(x, '--', simplify = T)[[1]]))
samples <-sce.seurat@assays$originalexp@data@Dimnames[[2]]
# Create a new seurat object using only the original counts.
x <-sce.seurat@assays$originalexp@counts
rownames(x)<-genes
sce.seurat<-CreateSeuratObject(counts = x)
sce.seurat$batch<-sce.data$batch
sce.seurat$condition<-sce.data$conditionQC
RUN minimal QC (mitochondrial fraction, number of genes etc.)
qc.dir <- file.path(results.dir, 'qc')
dir.create(qc.dir)
sce.seurat[["percent.mt"]] <- PercentageFeatureSet(sce.seurat, pattern = "^mt-")
p<-VlnPlot(sce.seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
ggsave(p, file = file.path(qc.dir, 'before_qc.pdf'), height = 10)
fwrite(as.data.table(table(sce.seurat@meta.data['batch'])), file = file.path(qc.dir, 'cells_before_qc.tsv'), sep = '\t')
# Subset and plot after qc
sce.seurat <- subset(sce.seurat, subset = nFeature_RNA > 100 & nFeature_RNA < 6000 & percent.mt < 10)
p<-VlnPlot(sce.seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
ggsave(p, file = file.path(qc.dir, 'after_qc.pdf'), height = 10)
fwrite(as.data.table(table(sce.seurat@meta.data['batch'])), file = file.path(qc.dir, 'cells_after_qc.tsv'), sep = '\t')
p
Standard Seurat pipeline
Initial clustering
Here, a standard Seurat preprocessing pipeline is run with Data normalization, Variable Feature Extraction, PCA, Neighborhood embedding and UMAP representation.
sce.seurat <- NormalizeData(object = sce.seurat)
sce.seurat <- FindVariableFeatures(object = sce.seurat)
sce.seurat <- ScaleData(object = sce.seurat)
sce.seurat <- RunPCA(object = sce.seurat)
sce.seurat <- FindNeighbors(object = sce.seurat, k.param = 20)
sce.seurat <- FindClusters(object = sce.seurat)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1473
## Number of edges: 41812
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8728
## Number of communities: 13
## Elapsed time: 0 secondssce.seurat <- RunUMAP(object = sce.seurat, dims = 1:50)
DimPlot(object = sce.seurat, reduction = "umap")
Visualizing the clustering
standard.plots<-DimPlot(object = sce.seurat, group.by = c('seurat_clusters', 'batch', 'condition'))
ggsave(standard.plots, file = file.path(qc.dir, 'embeddings_before_harmony.pdf'), width = 25, units = 'cm', height = 10)
standard.plots+ plot_layout(ncol = 1)
There is no massive batch effect in the data, if any, there is a slight problem width SCE.005/ cluster 3. Since the other analysis used batch integration it is used here as well to avoid to strong differences.
To conform with prior analysis integrate the batches
We use harmony to integrate the batches.
sce.seurat <- RunHarmony(sce.seurat, group.by.vars = "batch")
sce.seurat <- RunUMAP(sce.seurat, reduction = "harmony", dims = 1:30)
sce.seurat <- FindNeighbors(sce.seurat, reduction = "harmony", dims = 1:30) %>% FindClusters()## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1473
## Number of edges: 49120
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8656
## Number of communities: 12
## Elapsed time: 0 secondsharmony.plot<-DimPlot(object = sce.seurat, group.by = c('seurat_clusters', 'batch', 'condition'))
ggsave(harmony.plot, file = file.path(qc.dir, 'embeddings_after_harmony.pdf'), width = 25, units = 'cm', height = 10)
harmony.plot+ plot_layout(ncol = 1)
Marker gene set base embedding
In order to annotate the clusters to phenotypes from the literature, we use known gene sets and score the data using the ModuleScore function of Seurat. There are several gene sets in the literature, and we use several to see which one fits the best.
markers<-fread(marker_file, header=FALSE)
iCAF <- markers[V2 == 'iCAF']$V1
myCAF<- markers[V2 == 'myCAF']$V1
apCAF<- markers[V2 == 'apCAF']$V1
mCAF<- markers[V2 == 'mCAF']$V1
dCAF<- markers[V2 == 'dCAF']$V1
vCAF<- markers[V2 == 'vCAF']$V1
score<-list(iCAF, myCAF, apCAF, mCAF, dCAF, vCAF)
names(score)<-c('iCAF', 'myCAF', 'apCAF', 'mCAF', 'dCAF', 'vCAF')
sce.seurat <- AddModuleScore(
object = sce.seurat,
features = score,
name =c('iCAF', 'myCAF', 'apCAF', 'mCAF', 'dCAF', 'vCAF'),
ctrl = 100)
p1<- FeaturePlot(sce.seurat, c('iCAF1'), pt.size= 3)+scale_color_gradient2(low = 'blue', mid = 'lightgrey',high = 'red', midpoint = 0) +ggtitle("")+xlab("")+ylab("")+theme(axis.ticks = element_blank(), axis.text = element_blank(), aspect.ratio = 1)+ guides(color="none")
p3<- FeaturePlot(sce.seurat, c('myCAF2'), pt.size= 3)+scale_color_gradient2(low = 'blue', mid = 'lightgrey',high = 'red', midpoint = 0)+ggtitle("")+xlab("")+ylab("")+theme(axis.ticks = element_blank(), axis.text = element_blank(), aspect.ratio = 1)+ guides(color="none")+ coord_fixed(ratio = 1)
p4<- FeaturePlot(sce.seurat, c('apCAF3'), pt.size= 3)+scale_color_gradient2(low = 'blue', mid = 'lightgrey',high = 'red', midpoint = 0)+ggtitle("")+xlab("")+ylab("")+theme(axis.ticks = element_blank(), axis.text = element_blank(), aspect.ratio = 1)+ guides(color="none")+ coord_fixed(ratio = 1)
p2<- FeaturePlot(sce.seurat, c('mCAF4'), pt.size= 3)+scale_color_gradient2(low = 'blue', mid = 'lightgrey',high = 'red', midpoint = 0)+ggtitle("")+xlab("")+ylab("")+theme(axis.ticks = element_blank(), axis.text = element_blank(), aspect.ratio = 1)+ guides(color="none")+ coord_fixed(ratio = 1)
p5<- FeaturePlot(sce.seurat, c('dCAF5'), pt.size= 3)+scale_color_gradient2(low = 'blue', mid = 'lightgrey',high = 'red', midpoint = 0)+ggtitle("")+xlab("")+ylab("")+theme(axis.ticks = element_blank(), axis.text = element_blank(), aspect.ratio = 1)+ guides(color="none")+ coord_fixed(ratio = 1)
p6<- FeaturePlot(sce.seurat, c('vCAF6'), pt.size= 3)+scale_color_gradient2(low = 'blue', mid = 'lightgrey',high = 'red', midpoint = 0)+ggtitle("")+xlab("")+ylab("")+ coord_fixed(ratio = 1)+theme(axis.ticks = element_blank(), axis.text = element_blank(), aspect.ratio = 1)+ guides(color="none")
cellstates.dir<-file.path(results.dir, 'cellstates')
dir.create(cellstates.dir)
pp<-p1+p3+p4+p2+p5+p6+ plot_layout(ncol = 2)
ggsave(pp, file = file.path(cellstates.dir, 'cellstates.pdf'), height = 30, units = 'cm')
pp
pnames<-c('iCAF', 'myCAF', 'apCAF', 'mCAF', 'dCAF', 'vCAF')
i<-1
for (p in list(p1, p2, p3, p4, p5, p6)){
ggsave(p, file = file.path(cellstates.dir, paste0(pnames[i], '.pdf')), height = 8, width = 8)
ggsave(p, file = file.path(cellstates.dir, paste0(pnames[i], '.png')), height = 8, width = 8)
i<-i+1
}Differential abundance using DASeq
In order to identify whether the proportion of different cell populations shift between WT and KO, we run differential abundance analysis using DASEQ. DASeq is able to identify groups of cells which contain systematically more cells of one of the conditions. Here, DAseq is set up such that low scores in the DASeq result correspond to higher number of WT cells and high scores indicate higher numbers of KO cells. The analysis here follows the standard procedure as indicated in the tutorial: https://klugerlab.github.io/DAseq/articles/tutorial.html
Set up data
daseq.dir<-file.path(results.dir, 'daseq')
dir.create(daseq.dir)
sce.pca<-sce.seurat@reductions$pca@cell.embeddings
sce.labels<-rownames(sce.pca)
lables.WT<- rownames(sce.seurat@meta.data[sce.seurat$condition=='WT' ,])
lables.KO<- rownames(sce.seurat@meta.data[sce.seurat$condition=='KO' ,])
da_cells <- getDAcells(
X = sce.pca,
cell.labels = sce.labels,
labels.1 = lables.WT,
labels.2 = lables.KO,
k.vector = seq(15, 500, 50),
plot.embedding = sce.seurat@reductions$umap@cell.embeddings[,1:2]
)## Calculating DA score vector.
## Running GLM.
## Test on random labels.
## Setting thresholds based on permutationda_cells <- updateDAcells(
X = da_cells, pred.thres = c(-0.8,0.8),
plot.embedding = sce.seurat@reductions$umap@cell.embeddings[,1:2]
)
p<-da_cells$da.cells.plot
ggsave(p, file = file.path(daseq.dir, 'daseq_cells.pdf'))Show the statistically enriched regions
da_regions <- getDAregion(
X = sce.pca,
da.cells = da_cells,
cell.labels = sce.labels,
labels.1 = lables.WT,
labels.2 = lables.KO,
resolution = 0.01,
min.cell = 15,
plot.embedding = sce.seurat@reductions$umap@cell.embeddings[,1:2],
size= 3
)## Removing 4 DA regions with cells < 15.tab<-da_regions$DA.stat
fwrite(tab, file.path(daseq.dir, 'daseq_stats.tsv'))
tab## DA.score pval.wilcoxon pval.ttest
## [1,] 0.9351875 5.369764e-14 7.701717e-15
## [2,] 1.0000000 1.347801e-12 2.272813e-13
## [3,] -1.0000000 5.623935e-06 1.166155e-05There are 3 regions which show a statistical enrichment of either WT or KO cell types.
Make plots
okabe <- c("grey", "red", "black", "blue")
pp<-da_regions$da.region.plot+scale_color_manual (values = okabe)+xlab("")+ylab("")+theme(axis.ticks = element_blank(), axis.text = element_blank(), aspect.ratio = 1)
pp
ggsave(pp, file = file.path(daseq.dir, 'daseq_regions_wiht_number.pdf'), height = 8, width = 9)
ggsave(pp, file = file.path(daseq.dir, 'daseq_regions_wiht_number.png'), height = 8, width = 9, bg = "white")
pp<-pp+scale_color_manual(labels = c("neutral", "KO", "KO", 'WT'), values = okabe)+xlab("")+ylab("")+theme(axis.ticks = element_blank(), axis.text = element_blank())+ coord_fixed(ratio = 1)
ggsave(pp, file = file.path(daseq.dir, 'daseq_regions_wiht_label.pdf'), height = 8, width = 9)
pp
We see that the black KO region corresponds to a cluster which contains mainly M005 cells. It is hard to evaluate whether this group of cells represents an artifact or a genuine discovery. Here, the choice was made to disregard this cluster from the analysis, because it may not represent a generalizable discovery for all KO mice.
Find markers for DA regions
Find representative genes for the regions.
STG_markers <- STGmarkerFinder(
X = sce.seurat@assays$RNA@data,
da.regions = da_regions,
lambda = 1.5, n.runs = 5, return.model = T,
python.use = python_env,
GPU = ''
)## num_samples : 114
## Epoch: 1000 train loss=1.046695828 valid loss= 1.043797851 valid acc= 1.000000000
## train reg_fs: 0.9917962551116943
## Epoch: 2000 train loss=0.974549890 valid loss= 0.975865960 valid acc= 1.000000000
## train reg_fs: 0.9371622800827026
## Epoch: 3000 train loss=0.908763587 valid loss= 0.919636726 valid acc= 1.000000000
## train reg_fs: 0.8771899938583374
## Optimization Finished!
## test loss: 0.9196367263793945, test acc: 1.0
## 716
## num_samples : 114
## Epoch: 1000 train loss=0.844590425 valid loss= 0.885044634 valid acc= 1.000000000
## train reg_fs: 0.8142205476760864
## Epoch: 2000 train loss=0.765868723 valid loss= 0.827748299 valid acc= 0.944444418
## train reg_fs: 0.7495274543762207
## Epoch: 3000 train loss=0.701694310 valid loss= 0.760835052 valid acc= 1.000000000
## train reg_fs: 0.6849058866500854
## Optimization Finished!
## test loss: 0.7608350515365601, test acc: 1.0
## 97
## num_samples : 114
## Epoch: 1000 train loss=0.651630700 valid loss= 0.665259421 valid acc= 1.000000000
## train reg_fs: 0.6218627691268921
## Epoch: 2000 train loss=0.582408190 valid loss= 0.606567800 valid acc= 1.000000000
## train reg_fs: 0.561877965927124
## Epoch: 3000 train loss=0.519773245 valid loss= 0.553079665 valid acc= 1.000000000
## train reg_fs: 0.5060113668441772
## Optimization Finished!
## test loss: 0.5530796647071838, test acc: 1.0
## 65
## num_samples : 114
## Epoch: 1000 train loss=0.501071990 valid loss= 0.539566278 valid acc= 0.944444418
## train reg_fs: 0.4571804702281952
## Epoch: 2000 train loss=0.467387408 valid loss= 0.486466348 valid acc= 0.944444418
## train reg_fs: 0.41316932439804077
## Epoch: 3000 train loss=0.397382945 valid loss= 0.437929630 valid acc= 0.944444418
## train reg_fs: 0.3737944960594177
## Optimization Finished!
## test loss: 0.437929630279541, test acc: 0.9444444179534912
## 54
## num_samples : 114
## Epoch: 1000 train loss=0.351103395 valid loss= 0.373511761 valid acc= 1.000000000
## train reg_fs: 0.33759021759033203
## Epoch: 2000 train loss=0.318421364 valid loss= 0.340284646 valid acc= 1.000000000
## train reg_fs: 0.30549705028533936
## Epoch: 3000 train loss=0.287361622 valid loss= 0.312105000 valid acc= 1.000000000
## train reg_fs: 0.2773233950138092
## Optimization Finished!
## test loss: 0.3121050000190735, test acc: 1.0
## 48
## num_samples : 84
## Epoch: 1000 train loss=1.000762701 valid loss= 1.001584053 valid acc= 1.000000000
## train reg_fs: 0.987221360206604
## Epoch: 2000 train loss=0.936720610 valid loss= 0.943250418 valid acc= 1.000000000
## train reg_fs: 0.9279508590698242
## Epoch: 3000 train loss=0.872979283 valid loss= 0.901977658 valid acc= 1.000000000
## train reg_fs: 0.863950788974762
## Optimization Finished!
## test loss: 0.9019776582717896, test acc: 1.0
## 536
## num_samples : 84
## Epoch: 1000 train loss=0.813439727 valid loss= 0.989051878 valid acc= 0.949999988
## train reg_fs: 0.7984195947647095
## Epoch: 2000 train loss=0.739849627 valid loss= 1.161535978 valid acc= 0.750000000
## train reg_fs: 0.7314081788063049
## Epoch: 3000 train loss=0.672549307 valid loss= 1.184564590 valid acc= 0.600000024
## train reg_fs: 0.6649122834205627
## Optimization Finished!
## test loss: 1.1845645904541016, test acc: 0.6000000238418579
## 49
## num_samples : 84
## Epoch: 1000 train loss=0.616423845 valid loss= 1.061165094 valid acc= 0.649999976
## train reg_fs: 0.6004747152328491
## Epoch: 2000 train loss=0.549860060 valid loss= 0.878835440 valid acc= 0.800000012
## train reg_fs: 0.5398499965667725
## Epoch: 3000 train loss=0.494013548 valid loss= 0.721590817 valid acc= 0.949999988
## train reg_fs: 0.4839315712451935
## Optimization Finished!
## test loss: 0.7215908169746399, test acc: 0.949999988079071
## 26
## num_samples : 84
## Epoch: 1000 train loss=0.442976654 valid loss= 0.619562447 valid acc= 0.949999988
## train reg_fs: 0.4334219992160797
## Epoch: 2000 train loss=0.394765377 valid loss= 0.523002207 valid acc= 0.949999988
## train reg_fs: 0.3881967067718506
## Epoch: 3000 train loss=0.354995072 valid loss= 0.440699160 valid acc= 0.949999988
## train reg_fs: 0.3481259047985077
## Optimization Finished!
## test loss: 0.44069916009902954, test acc: 0.949999988079071
## 22
## num_samples : 84
## Epoch: 1000 train loss=0.322823763 valid loss= 0.374619067 valid acc= 1.000000000
## train reg_fs: 0.31274858117103577
## Epoch: 2000 train loss=0.289415061 valid loss= 0.327162355 valid acc= 1.000000000
## train reg_fs: 0.2817081809043884
## Epoch: 3000 train loss=0.261621624 valid loss= 0.288426757 valid acc= 1.000000000
## train reg_fs: 0.2545800507068634
## Optimization Finished!
## test loss: 0.2884267568588257, test acc: 1.0
## 20
## num_samples : 32
## Epoch: 1000 train loss=1.034658551 valid loss= 1.159772754 valid acc= 1.000000000
## train reg_fs: 0.9919085502624512
## Epoch: 2000 train loss=0.953719437 valid loss= 1.069233060 valid acc= 1.000000000
## train reg_fs: 0.9354907274246216
## Epoch: 3000 train loss=0.892450809 valid loss= 1.001452565 valid acc= 1.000000000
## train reg_fs: 0.8736256957054138
## Optimization Finished!
## test loss: 1.0014525651931763, test acc: 1.0
## 737
## num_samples : 32
## Epoch: 1000 train loss=0.823999465 valid loss= 1.149300575 valid acc= 0.833333313
## train reg_fs: 0.8098279237747192
## Epoch: 2000 train loss=0.754330635 valid loss= 1.132464886 valid acc= 0.833333313
## train reg_fs: 0.7429871559143066
## Epoch: 3000 train loss=0.683510005 valid loss= 1.102321982 valid acc= 0.833333313
## train reg_fs: 0.6761536598205566
## Optimization Finished!
## test loss: 1.102321982383728, test acc: 0.8333333134651184
## 74
## num_samples : 32
## Epoch: 1000 train loss=0.642935932 valid loss= 0.885854781 valid acc= 0.833333313
## train reg_fs: 0.6176947951316833
## Epoch: 2000 train loss=0.595218062 valid loss= 0.835131049 valid acc= 0.833333313
## train reg_fs: 0.5611858367919922
## Epoch: 3000 train loss=0.531115532 valid loss= 0.771598816 valid acc= 0.833333313
## train reg_fs: 0.5077548027038574
## Optimization Finished!
## test loss: 0.7715988159179688, test acc: 0.8333333134651184
## 52
## num_samples : 32
## Epoch: 1000 train loss=0.469875276 valid loss= 0.703926504 valid acc= 0.833333313
## train reg_fs: 0.4567704200744629
## Epoch: 2000 train loss=0.422494829 valid loss= 0.654970527 valid acc= 0.833333313
## train reg_fs: 0.4104986786842346
## Epoch: 3000 train loss=0.379862428 valid loss= 0.603472531 valid acc= 0.833333313
## train reg_fs: 0.36914631724357605
## Optimization Finished!
## test loss: 0.6034725308418274, test acc: 0.8333333134651184
## 41
## num_samples : 32
## Epoch: 1000 train loss=0.339740425 valid loss= 0.525905669 valid acc= 0.833333313
## train reg_fs: 0.3319685161113739
## Epoch: 2000 train loss=0.304611951 valid loss= 0.483825505 valid acc= 0.833333313
## train reg_fs: 0.29925212264060974
## Epoch: 3000 train loss=0.276807666 valid loss= 0.446785033 valid acc= 0.833333313
## train reg_fs: 0.27059632539749146
## Optimization Finished!
## test loss: 0.446785032749176, test acc: 0.8333333134651184
## 38dp<-DimPlot(sce.seurat, split.by = 'batch' , ncol = 3)
ggsave(dp, file = file.path(daseq.dir, 'batches.split.pdf'))
label<-function(x, threshold){
if (x>threshold){
return ('up')
}
else if(x< -threshold) {
return ('down')
}
else{
return('neutral')
}
}
sce.seurat$da.pred<-da_cells$da.pred
sce.seurat$da.pred.bin<-round(da_cells$da.pred)
sce.seurat$da.pred.bin.label<-sapply(sce.seurat$da.pred, function(x) label(x, 0.8))
da_clusters<-list(c(0), c(8), c(7,9))
for (i in 1:length(STG_markers$da.markers)){
region_2<-as.data.table(STG_markers$da.markers[[as.character(i)]])
region_2<-region_2[p_value<0.01]
print(region_2)
fwrite(region_2, file = file.path(daseq.dir, paste0('region_', i,'.tsv')))
Idents(sce.seurat)<-'seurat_clusters'
sce.sub<-subset(sce.seurat, subset = seurat_clusters %in% da_clusters[[i]] )
pv<-VlnPlot(sce.sub, features = region_2[1:10]$gene , group.by ='da.pred.bin.label' )
ggsave(pv, file = file.path(daseq.dir, paste0('top_region_', i,'.pdf')), width = 20, height = 30, units = 'cm')
}## gene avg_logFC p_value
## 1: Apod 3.6087416 4.400616e-61
## 2: Osr2 3.3450601 3.869057e-56
## 3: Gas1 2.8932365 2.630517e-54
## 4: Fst 3.2248433 6.781890e-54
## 5: Ctla2a 4.5678976 9.968275e-50
## ---
## 306: C1d 0.9995466 7.453457e-03
## 307: Gstm5 1.0829472 8.089937e-03
## 308: Blvrb 1.2124816 8.273648e-03
## 309: Glud1 0.4116755 8.820157e-03
## 310: Strn3 0.4098708 9.151682e-03## gene avg_logFC p_value
## 1: Trim29 6.44231783 9.611054e-227
## 2: Pkp1 5.49426401 6.880149e-190
## 3: Celsr2 5.80215778 5.733870e-186
## 4: Dsc3 5.18794218 8.509447e-176
## 5: 1190003J15Rik 5.61830914 4.003599e-172
## ---
## 466: Ppp3ca 0.34634591 6.926656e-03
## 467: Tpm2 0.08361507 7.070530e-03
## 468: Junb -0.22190389 7.584926e-03
## 469: Anxa2 -0.19999891 8.417456e-03
## 470: Rn45s -0.10410975 9.372423e-03## gene avg_logFC p_value
## 1: Mmp10 3.2565049 1.300282e-25
## 2: Col6a4 3.2625322 1.086675e-20
## 3: Fbn2 3.7727421 1.072979e-18
## 4: Foxl1 3.5903656 2.891393e-18
## 5: Inhba 3.0842734 7.246118e-18
## ---
## 187: Pim1 -1.4133427 8.787636e-03
## 188: Ndufs3 1.2077530 9.362957e-03
## 189: Acta2 0.7538293 9.619177e-03
## 190: Eif3a 0.6812234 9.756517e-03
## 191: Cox6a1 0.5890282 9.906662e-03mdf<-as.data.table(sce.seurat@meta.data)
summary.da<- mdf[, .N, by = c( 'seurat_clusters','condition', 'da.pred.bin.label')] %>% arrange(seurat_clusters)
summary.da<-as.data.table(summary.da)
fwrite(summary.da, file.path(daseq.dir, 'summary_daseq.tsv'), sep = '\t')Differential Gene expression global WT vs. KO
To get some general understanding we run a differential analysis of all WT vs all KO cells.
de.dir<-file.path(results.dir, 'deg')
dir.create(de.dir)
Idents(sce.seurat)<-'condition'
m<-FindMarkers(sce.seurat, ident.1 = 'WT', ident.2 = 'KO')
genes<-rownames(m)
m<-as.data.table(m)
m$gene<-genes
fwrite(m[p_val_adj<0.01], file = file.path(de.dir, 'DEG_WT_KO_p0.01.tsv'), sep = '\t')
vln<-VlnPlot(sce.seurat, features = m$gene[1:10])
ggsave(vln, file = file.path(de.dir, 'vln_plots_top10.pdf'))Make plots and assign tentative phenotypes
In order to characterize cell states further, it is useful to assign tentative labels to each group, even if a categorical cell type classification/definition most likely does not reflect reality. Based on the cell scores above, we used the mcaf, icaf, apcaf definition also used by Elyada et. al to characterize our cells. Additionally a cluster of vCAF cells has been identified.
##############################################################################
## 6. grouped UMAP
numbered.umap<-DimPlot(sce.seurat, group.by = 'seurat_clusters', label = TRUE)
numbered.umap
ggsave(numbered.umap, file = file.path(de.dir, 'umap_with_number.pdf'))
numbered.umap<-DimPlot(sce.seurat, group.by = 'condition', label = TRUE)
numbered.umap
ggsave(numbered.umap, file = file.path(de.dir, 'umap_with_condition.pdf'))
## Assign labels.
sce.seurat$phenotype<-recode(sce.seurat$seurat_clusters, '0' = 'mCAF-iCAF', '1' = 'apCAF', '2'= 'apCAF', '3' = 'mCAF-iCAF', '4' = 'vCAF', '5' = 'mCAF-iCAF', '6' = 'apCAF' ,'7' = 'mCAF-iCAF', '8'='??-0', '9'='mCAF-iCAF', '10' = '??-1', '11' = 'apCAF')
sce.seurat$phenotype_new<-recode(sce.seurat$seurat_clusters, '0' = 'iCAF1', '1' = 'apCAF', '2'= 'apCAF', '3' = 'iCAF2', '4' = 'vCAF', '5' = 'i/myCAF', '6' = 'apCAF' ,'7' = 'myCAF', '8'='ND.3', '9'='myCAF', '10' = 'ND.1', '11' = 'ND.2')
##Plot
labelle.umap<-DimPlot(sce.seurat, group.by = 'phenotype', label = TRUE)
ggsave(labelle.umap, file = file.path(de.dir, 'labelled_umap.pdf'))
labelle.umap<-DimPlot(sce.seurat, group.by = 'phenotype_new', label = TRUE, label.size = 3, label.box = T, repel = T, pt.size = 3)+xlab("")+ylab("")+theme(axis.text = element_blank(), axis.ticks = element_blank())+ggtitle("")+xlab("")+ylab("")+theme(axis.ticks = element_blank(), axis.text = element_blank(), aspect.ratio = 1)+ guides(color="none")
labelle.umap
ggsave(labelle.umap, file = file.path(de.dir, 'labelled_umap.pdf'), height = 8, width = 8)
ggsave(labelle.umap, file = file.path(de.dir, 'labelled_umap.png'), height = 8, width = 8)Phenotype scores.
The ModuleScores calculated above shown for each of the cell phenotype groups.
Idents(sce.seurat)<-'phenotype_new'
icaf1.score<-VlnPlot(sce.seurat, features = 'iCAF1', idents = c('iCAF1', 'iCAF2', 'i/myCAF','myCAF'), split.by = 'condition')+ggtitle('iCAF scores in phenotype clusters')+xlab('')
icaf1.score
ggsave(icaf1.score, file = file.path(cellstates.dir, 'iCAF_scores_in_pops.pdf'), width = 30, height = 10, units = 'cm')
Idents(sce.seurat)<-'phenotype_new'
mycaf.score<-VlnPlot(sce.seurat, features = 'myCAF2', idents = c('iCAF1', 'iCAF2', 'i/myCAF','myCAF'), split.by = 'condition') +ggtitle('myCAF scores in phenotype clusters')+xlab('')
mycaf.score
ggsave(mycaf.score, file = file.path(cellstates.dir, 'myCAF_scores_in_pops.pdf'), width = 30, height = 10, units = 'cm')Correlation between scores by Seurat cluster
meta<-as.data.table(sce.seurat@meta.data)
for (g in list(c(7,9), c(3,5), c(0))){
for(phe in c('WT', 'KO')){
co<-cor.test(meta[condition == phe & seurat_clusters %in% g]$iCAF1, meta[condition == phe & seurat_clusters %in% g]$myCAF2)
chars<-capture.output(co)
fileConn<-file(file.path(cellstates.dir, paste0('corr.', paste(g, collapse = ''), '.', phe, '.txt')))
writeLines(chars, con =fileConn )
close(fileConn)
}
}
scatter<-ggplot(meta, aes(iCAF1, myCAF2))+geom_point()+facet_wrap(c("phenotype_new", "condition"), scales = 'free')
ggsave(scatter, file = file.path(cellstates.dir, 'scatter_correlations.pdf'), height = 30, width = 30, units = 'cm')
Idents(sce.seurat)<-'seurat_clusters'
icaf.score<-VlnPlot(sce.seurat, features = 'iCAF1', split.by = 'condition')+ggtitle('iCAF score by Leiden clusters')
ggsave(icaf.score, file = file.path(cellstates.dir, 'iCAF_scores_by_seurat.pdf'))
Idents(sce.seurat)<-'phenotype_new'
icaf.score<-VlnPlot(sce.seurat, features = 'iCAF1', split.by = 'condition')+ggtitle('iCAF score by Phenotype clusters')
icaf.score
ggsave(icaf.score, file = file.path(cellstates.dir, 'iCAF_scores_by_phenotype.pdf'))
summary_table<-table(sce.seurat@meta.data[, c("condition", "phenotype_new")])
fwrite(summary_table, file = file.path(de.dir, 'cells_per_pop.tsv'), sep= '\t')What are these clusters? ==> DEG between clusters.
Cluster iCAF1 vs all other iCAF/myCAF
Idents(sce.seurat)<-'phenotype_new'
mm<-FindMarkers(sce.seurat, ident.1 = 'iCAF1', ident.2 = c('myCAF', 'i/myCAF', 'iCAF2'))
g<-rownames(mm)
mm<-as.data.table(mm)
mm$genes<-g
fwrite(mm, file = file.path(de.dir, 'cluster0_vs_restofmCAFiCAF.tsv'), sep = '\t')Cluster mycaf vs i/myCAF
Idents(sce.seurat)<-'phenotype_new'
mm<-FindMarkers(sce.seurat, ident.1 = 'i/myCAF' , ident.2 = c('myCAF'))
g<-rownames(mm)
mm<-as.data.table(mm)
mm$genes<-g
fwrite(mm, file = file.path(de.dir, 'clusterimyCAF_vs_clustermyCAF.tsv'), sep = '\t')9. Add tests between WT and KO within clusters.
all_genes<-c()
Idents(sce.seurat)<-'phenotype_new'
phenotypes<-c('myCAF', 'i/myCAF', 'iCAF2', 'vCAF', 'iCAF1', 'apCAF')
for (pop in phenotypes){
sub.seurat<-subset(sce.seurat, subset = phenotype_new == pop)
Idents(sub.seurat)<-'condition'
m<-FindMarkers(sub.seurat, ident.1 = 'WT', ident.2 = 'KO')
genes<-rownames(m)
m<-as.data.table(m)
m$genes<-genes
m<-m[p_val_adj < 0.05]
pop <- gsub('/', '_', pop)
fwrite(m, file = file.path(de.dir, paste0(pop, 'DEG_WT_KO_p0.05.tsv')), sep = '\t')
sub.dir<-file.path(de.dir, pop)
dir.create(sub.dir)
for(g in m$genes){
pl<-VlnPlot(sub.seurat, features = g)
ggsave(pl, file = file.path(sub.dir, paste0(g, '.pdf')))
}
all_genes<-c(all_genes, m$genes)
} Save data object
saveRDS(sce.seurat, file = file.path(results.dir, 'seurat.Rds'), compress = F)
sce.seurat <- readRDS(file.path(results.dir, 'seurat.Rds'))Check transcription factor activity
tf.dir<-file.path(results.dir, 'tf-target')
dir.create(tf.dir)
net <- get_collectri(organism='mouse', split_complexes=FALSE)
fwrite(as.data.table(net), file = file.path(tf.dir, 'collectri.tsv'), sep = '\t')
net1<- net %>% filter(source =='Zeb1')
net2<- net %>% filter(source %in% net1$target)
apcaf.de<-fread(file.path(de.dir, 'apCAFDEG_WT_KO_p0.05.tsv'))
vcaf.de<-fread(file.path(de.dir, 'vCAFDEG_WT_KO_p0.05.tsv'))
icafmcaf.de.0<-fread(file.path(de.dir, 'iCAF1DEG_WT_KO_p0.05.tsv'))
icafmcaf.de.35<-fread(file.path(de.dir, 'i_myCAFDEG_WT_KO_p0.05.tsv'))
icafmcaf.de.79<-fread(file.path(de.dir, 'myCAFDEG_WT_KO_p0.05.tsv'))
icafmcaf.de.1<-fread(file.path(de.dir, 'iCAF2DEG_WT_KO_p0.05.tsv'))
intersect(net1$target, apcaf.de$genes)## character(0)intersect(net1$target, icafmcaf.de.0$genes)## [1] "Cdh13"p.prim<-VlnPlot(sce.seurat, features =intersect(net1$target, icafmcaf.de.0$genes), split.by = 'condition')
ggsave(p.prim, file = file.path(tf.dir, 'primary_targets_in_iCAFmCAF0.pdf'))
intersect(net1$target, vcaf.de$genes)## character(0)intersect(net2$target, apcaf.de$genes)## character(0)intersect(net2$target, icafmcaf.de.0$genes)## [1] "C3" "Egr1" "Cd44" "Fos" "Krt14" "Zeb1"
## [7] "Dhcr24" "Des" "Fosb" "Lcn2" "Serpinh1" "Nt5e"
## [13] "Hk2"intersect(net2$target, icafmcaf.de.1$genes)## character(0)intersect(net2$target, icafmcaf.de.35$genes)## character(0)intersect(net2$target, icafmcaf.de.79$genes)## character(0)intersect(net2$target, vcaf.de$genes)## [1] "S100a6"targets.targets.zeb1.icaf<-VlnPlot(sce.seurat, features = intersect(net2$target, icafmcaf.de.0$genes), split.by = 'condition')
ggsave(targets.targets.zeb1.icaf, file = file.path(tf.dir, 'targets_of_targets_of_zeb1_icaf_mcaf_0.pdf'), height = 60, width = 40, units = 'cm', limitsize = F)
targets.targets.zeb1.vacf<-VlnPlot(sce.seurat, features = intersect(net2$target, vcaf.de$genes), split.by = 'condition')
ggsave(targets.targets.zeb1.vacf, file = file.path(tf.dir, 'targets_of_targets_of_zeb1_vCAF.pdf'), height = 60, width = 40, units = 'cm', limitsize = F)sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=de_DE.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=de_DE.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] pheatmap_1.0.12 tidyr_1.3.0
## [3] tibble_3.2.1 dorothea_1.6.0
## [5] decoupleR_2.7.1 dplyr_1.1.3
## [7] data.table_1.14.8 harmony_1.1.0
## [9] Rcpp_1.0.11 DAseq_1.0.0
## [11] patchwork_1.1.3 stringr_1.5.0
## [13] SeuratObject_4.1.4 Seurat_4.4.0
## [15] scran_1.22.1 scater_1.22.0
## [17] ggplot2_3.4.4 scuttle_1.4.0
## [19] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0
## [21] Biobase_2.54.0 GenomicRanges_1.46.1
## [23] GenomeInfoDb_1.30.1 IRanges_2.28.0
## [25] S4Vectors_0.32.4 BiocGenerics_0.40.0
## [27] MatrixGenerics_1.6.0 matrixStats_1.0.0
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 scattermore_1.2
## [3] ModelMetrics_1.2.2.2 ragg_1.2.6
## [5] bit64_4.0.5 knitr_1.44
## [7] irlba_2.3.5.1 DelayedArray_0.20.0
## [9] rpart_4.1.16 hardhat_1.3.0
## [11] RCurl_1.98-1.12 generics_0.1.3
## [13] ScaledMatrix_1.2.0 RhpcBLASctl_0.23-42
## [15] cowplot_1.1.1 RANN_2.6.1
## [17] proxy_0.4-27 future_1.33.0
## [19] bit_4.0.5 tzdb_0.4.0
## [21] spatstat.data_3.0-1 xml2_1.3.5
## [23] lubridate_1.9.3 httpuv_1.6.11
## [25] viridis_0.6.4 gower_1.0.1
## [27] xfun_0.40 hms_1.1.3
## [29] jquerylib_0.1.4 evaluate_0.22
## [31] promises_1.2.1 fansi_1.0.5
## [33] progress_1.2.2 readxl_1.4.3
## [35] igraph_1.5.1 htmlwidgets_1.6.2
## [37] spatstat.geom_3.2-5 purrr_1.0.2
## [39] ellipsis_0.3.2 selectr_0.4-2
## [41] backports_1.4.1 deldir_1.0-9
## [43] sparseMatrixStats_1.6.0 vctrs_0.6.4
## [45] ROCR_1.0-11 abind_1.4-5
## [47] caret_6.0-94 cachem_1.0.8
## [49] withr_2.5.1 progressr_0.14.0
## [51] checkmate_2.3.0 vroom_1.6.4
## [53] sctransform_0.4.0 prettyunits_1.2.0
## [55] goftest_1.2-3 cluster_2.1.2
## [57] lazyeval_0.2.2 crayon_1.5.2
## [59] spatstat.explore_3.2-3 glmnet_4.1-8
## [61] edgeR_3.36.0 recipes_1.0.8
## [63] pkgconfig_2.0.3 labeling_0.4.3
## [65] nlme_3.1-155 vipor_0.4.5
## [67] nnet_7.3-17 rlang_1.1.1
## [69] globals_0.16.2 lifecycle_1.0.3
## [71] miniUI_0.1.1.1 rsvd_1.0.5
## [73] cellranger_1.1.0 polyclip_1.10-6
## [75] lmtest_0.9-40 Matrix_1.6-1.1
## [77] zoo_1.8-12 beeswarm_0.4.0
## [79] ggridges_0.5.4 png_0.1-8
## [81] viridisLite_0.4.2 bitops_1.0-7
## [83] KernSmooth_2.23-20 pROC_1.18.4
## [85] DelayedMatrixStats_1.16.0 shape_1.4.6
## [87] parallelly_1.36.0 spatstat.random_3.1-6
## [89] readr_2.1.4 beachmat_2.10.0
## [91] scales_1.2.1 magrittr_2.0.3
## [93] plyr_1.8.9 ica_1.0-3
## [95] zlibbioc_1.40.0 compiler_4.1.2
## [97] dqrng_0.3.1 RColorBrewer_1.1-3
## [99] fitdistrplus_1.1-11 cli_3.6.1
## [101] XVector_0.34.0 listenv_0.9.0
## [103] pbapply_1.7-2 MASS_7.3-55
## [105] tidyselect_1.2.0 stringi_1.7.12
## [107] textshaping_0.3.7 yaml_2.3.7
## [109] BiocSingular_1.10.0 locfit_1.5-9.8
## [111] ggrepel_0.9.4 grid_4.1.2
## [113] sass_0.4.7 bcellViper_1.30.0
## [115] tools_4.1.2 timechange_0.2.0
## [117] future.apply_1.11.0 parallel_4.1.2
## [119] rstudioapi_0.15.0 bluster_1.4.0
## [121] foreach_1.5.2 metapod_1.2.0
## [123] gridExtra_2.3 prodlim_2023.08.28
## [125] farver_2.1.1 Rtsne_0.16
## [127] digest_0.6.33 shiny_1.7.5.1
## [129] lava_1.7.2.1 later_1.3.1
## [131] RcppAnnoy_0.0.21 httr_1.4.7
## [133] colorspace_2.1-0 rvest_1.0.3
## [135] tensor_1.5 reticulate_1.34.0
## [137] splines_4.1.2 uwot_0.1.16
## [139] statmod_1.5.0 OmnipathR_3.11.1
## [141] spatstat.utils_3.0-3 sp_2.1-1
## [143] plotly_4.10.2 systemfonts_1.0.5
## [145] xtable_1.8-4 jsonlite_1.8.7
## [147] timeDate_4022.108 ipred_0.9-14
## [149] R6_2.5.1 pillar_1.9.0
## [151] htmltools_0.5.6.1 mime_0.12
## [153] glue_1.6.2 fastmap_1.1.1
## [155] BiocParallel_1.28.3 BiocNeighbors_1.12.0
## [157] class_7.3-20 codetools_0.2-18
## [159] utf8_1.2.3 lattice_0.20-45
## [161] bslib_0.5.1 spatstat.sparse_3.0-2
## [163] logger_0.2.2 curl_5.1.0
## [165] ggbeeswarm_0.7.2 leiden_0.4.3
## [167] survival_3.2-13 limma_3.50.3
## [169] rmarkdown_2.25 munsell_0.5.0
## [171] e1071_1.7-13 GenomeInfoDbData_1.2.7
## [173] iterators_1.0.14 reshape2_1.4.4
## [175] gtable_0.3.4