Figure S8

#Script used to make the figures used in supplementary figure S8

#Libraries needed
library(SCpubr)
library(Seurat)
library(Signac)
library(ggplot2)
library(dplyr)
library(ComplexHeatmap)
library(UCell)
library(circlize)

#Set the directory to where you have your data or want to store your output
setwd("/home/ipaassen/ATRT_sc_atlas/")

#Load seurat object and other datasets needed
ATRT <- readRDS("/data/projects/p868_CRISPR-RNA_TP53-RB1/irene/231207_ATRT_multiome_EnrichmentUMAP_newclusters.RDS")
motifs <- readRDS("motif_to_genesymbol.RDS") #To convert motif IDs to TF binding to it
gene_markers <- readRDS("10X_v3_frozen_samples_normalized_TME_and_TB_annotated_integrated_with_metacell_mapping_reannotated_top100_markers.rds") #Marker genes as identified in analysis pipeline of the snRNA dataset

#Define colour coding
subtype_colors <- c("ATRT_SHH" = "#253A79", "ATRT_MYC" = "#056A37", "ATRT_TYR" = "#931A1D", "ecMRT_BrainMet" = "#10C663")
idents_colors <- c("SHH.unspecified" = "#A4C2D6","CP.like" = "#be660e","Cilia.like" = "#be0e0e","OPC.like" = "#0435c8", 
                   "Rest" = "#84D3E5", 'NPC.like' = "#0466c8",'Hypoxic' = "#2A7072",
                   'Mesenchymal.like' = "#0ebe66","IPC.like" = "#be920e", 
                   "MYC.TYR.unspecified" = "#D6A4A9", "RG.like" = "#0497c8")
TF_colours <- c("MYC.Mes" = "#024431", "SHH.OPC" =  "#211D5C", "SHH.NPC" = "#20428A", "TYR.cilia" = "#8C1730", "SHH.shared" = "#0078BB")



#
#
############### Script for figures #######################
#
#


#Supplementary Figure 6, Panel A
#
#Stacked barplot of identities represented ber patient sample of the snMultiome data
do_BarPlot(ATRT, split.by = 'sample1', group.by = 'Enrich.scores',colors.use = idents_colors, 
           position = 'fill', order = F) 

#
#
############## Calculation of Differential motif activity per mature-like cells vs. IPC-like cells #########################
#
#


#Settings for differenital motif activityt calling
ATRT <- SetIdent(ATRT, value = "Enrich.scores")
DefaultAssay(ATRT) <- 'chromvar'

#
#
######### Supplementary Figure S8, Panel B #############
#
#

#Mesenchymal-like
Mesenchymal.TFs <- FindMarkers(
  ATRT,
  ident.1 = "IPC.like",
  ident.2 = "Mesenchymal.like",
  logfc.threshold = 0,
  min.pct = 0.1,
  only.pos = F)

#Make vulcanoplot of p-value and log2FC
SCpubr::do_VolcanoPlot(sample = ATRT,
                       de_genes = Mesenchymal.TFs,
                       pval_cutoff = 1e-40,
                       FC_cutoff = 5,
                       plot.title = "IPC.vs.Mesenchymal",add_gene_tags = F) + scale_color_brewer(palette = "OrRd", direction = -1)


#Filter for highest MYC marker genes
Mesenchymal.specific <-Mesenchymal.TFs[Mesenchymal.TFs$avg_log2FC < -5 & Mesenchymal.TFs$p_val_adj < 1e-40,]
Mesenchymal.specific$gene <- motifs$symbol[motifs$motif %in% rownames(Mesenchymal.specific)]

#
#
######### Supplementary Figure S8, Panel C #############
#
#

#NPC-like
NPC.TFs <- FindMarkers(
  ATRT,
  ident.1 = "IPC.like",
  ident.2 = "NPC.like",
  logfc.threshold = 0,
  min.pct = 0.1,
  only.pos = F,subset.ident = SHH)

#Make vulcanoplot of p-value and log2FC
SCpubr::do_VolcanoPlot(sample = ATRT,
                       de_genes = NPC.TFs,
                       pval_cutoff = 1e-70,
                       FC_cutoff = 1,
                       plot.title = "ATRT-SHH: IPC.vs.NPC",add_gene_tags = F) + scale_color_brewer(palette = "OrRd", direction = -1)


#Filter for highest SHH marker genes
NPC.specific <-NPC.TFs[NPC.TFs$avg_log2FC < -1 & NPC.TFs$p_val_adj < 1e-70,]
NPC.specific$gene <- motifs$symbol[motifs$motif %in% rownames(NPC.specific)]


#
#
######### Supplementary Figure S8, Panel D #############
#
#

# OPC-like vs IPC. like
OPC.TFs <- FindMarkers(
  ATRT,
  ident.1 = "IPC.like",
  ident.2 = c("OPC.like"),
  logfc.threshold = 0,
  min.pct = 0.1,
  only.pos = F)

# Make vulcanoplot of p-value and log2FC
SCpubr::do_VolcanoPlot(sample = ATRT,
                       de_genes = OPC.TFs,
                       pval_cutoff = 1e-10,
                       FC_cutoff =0.5,
                       plot.title = "ATRT-SHH: IPC.vs.NPC+OPC",add_gene_tags = F) + scale_color_brewer(palette = "OrRd", direction = -1)


#Filter for highest SHH marker genes
OPC.specific <-OPC.TFs[OPC.TFs$avg_log2FC < -0.5 & OPC.TFs$p_val_adj < 1e-10,]
OPC.specific$gene <- motifs$symbol[motifs$motif %in% rownames(OPC.specific)]


#
#
######### Supplementary Figure S8, Panel E #############
#
#
#

#Motifs against NPC and OPC both together
OPC.and.NPC <- list(NPC.like = rownames(NPC.specific), OPC.like = rownames(OPC.specific))
OPC.and.NPC.motifs <- intersect(rownames(NPC.specific),rownames(OPC.specific)) #Subset for motifs found in both

overlap_genes <- motifs$symbol[motifs$motif %in% OPC.and.NPC.motifs] #Get TFs behind motif ID
overlap_genes <- gsub("\\s*\\([^\\)]+\\)","",overlap_genes)
Group_order <- list("Enrich.scores" = c("IPC.like","RG.like" ,"NPC.like","OPC.like","SHH.unspecified","Mesenchymal.like","CP.like","Cilia.like","MYC.TYR.unspecified"))


#
#
######### Supplementary Figure S8, Panel F #############
#
#
#

# Expression Heatmap of the chromvar scores
do_ExpressionHeatmap(ATRT, features = OPC.and.NPC.motifs, groups.order = Group_order,axis.text.x.angle = 90,
                     group.by = "Enrich.scores", assay = "chromvar", enforce_symmetry = T, features.order = OPC.and.NPC.motifs, 
                     diverging.palette = "PuOr",legend.title = "Chromvar score")

# Expression heatmap of the TFs underlying the motif IDs
do_ExpressionHeatmap(ATRT, features = overlap_genes, features.order = overlap_genes, groups.order = Group_order,axis.text.x.angle = 90,
                     group.by = "Enrich.scores", assay = "RNA",slot = "scale.data", enforce_symmetry = T)


#
#
## Supplementary Figure S8 Panel G
#
#

#Make enrichUMAP of RFX3 gene expression
do_FeaturePlot(ATRT, features = c("RFX4"), reduction = "umap", 
               assay = 'RNA', order = T, legend.title = "LMX1A expression",
               legend.position = "bottom")

#
#
## Supplementary Figure S8 Panel H
#
#

#RFX4 is motif ID MA0799.1
do_FeaturePlot(ATRT, features = c("MA0799.1"), reduction = "umap", assay = 'chromvar', order = T, enforce_symmetry = F,
               legend.title = "LMX1A chromvar score",
               label.size = 4,
               use_viridis = T,
               viridis.palette = "inferno",
               viridis.direction = -1,
               min.cutoff = 0,
               legend.position = "bottom")



#
#
######### Supplementary Figure 6 Panel I #############
#
#


#Cilia-like
Cilia.TFs <- FindMarkers(
  ATRT,
  ident.1 = "IPC.like",
  ident.2 = "Cilia.like",
  logfc.threshold = 0,
  min.pct = 0.1,
  only.pos = F)


#Filter for highest TYR marker genes
Cilia.specific <-Cilia.TFs[Cilia.TFs$avg_log2FC < -1 & Cilia.TFs$p_val_adj < 1e-40,]
Cilia.genes <- motifs$symbol[motifs$motif %in% rownames(Cilia.specific)]
Cilia.specific$gene <- motifs$symbol[motifs$motif %in% rownames(Cilia.specific)]

order_genes <- unique(c(rownames(Cilia.specific)))
TYR_motifs <- motifs[motifs$motif %in% order_genes,]
TYR_motifs$symbol <- gsub("\\s*\\([^\\)]+\\)","",TYR_motifs$symbol)
Group_order <- list("Enrich.scores" = c("IPC.like","RG.like" ,"NPC.like","OPC.like","SHH.unspecified","Mesenchymal.like","CP.like","Cilia.like","MYC.TYR.unspecified"))

#heatmap on chromvar scores
do_ExpressionHeatmap(ATRT, features = TYR_motifs$motif, groups.order = Group_order,
                     group.by = "Enrich.scores", assay = "chromvar", enforce_symmetry = T, features.order = TYR_motifs$motif, 
                     diverging.palette = "PuOr",legend.title = "Chromvar score",axis.text.x.angle = 90)

#heatmap on TF expression 
do_ExpressionHeatmap(ATRT, features = TYR_motifs$symbol, features.order = TYR_motifs$symbol, groups.order = Group_order,
                     group.by = "Enrich.scores", assay = "RNA",slot = "scale.data", enforce_symmetry = T, axis.text.x.angle = 90)