10X Multiome Cell Lines

library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
library(GenomeInfoDb)
set.seed(1234)
library(stringr)
library(ggplot2)
library(infercnv)

setwd("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/")

# get gene annotations for hg38
annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
ucsc.levels <- str_replace(string=paste("chr",seqlevels(annotation),sep=""), pattern="chrMT", replacement="chrM")
seqlevels(annotation) <- ucsc.levels #standard seq level change threw error, using a string replace instead

# set up sample loop to load the RNA and ATAC data, save to seurat object
setupseurat<-function(i){
  setwd(paste0("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/",i,"/outs"))
  counts <- Read10X_h5("filtered_feature_bc_matrix.h5") #count data
  fragpath <- "atac_fragments.tsv.gz" #atac fragments
  metadata_cellranger<-read.csv("per_barcode_metrics.csv") #metadata
  row.names(metadata_cellranger)<-metadata_cellranger$barcode

  # create a Seurat object containing the RNA adata
  dat <- CreateSeuratObject(
    counts = counts$`Gene Expression`,
    assay = "RNA"
  )

  # create ATAC assay and add it to the object
  dat[["ATAC"]] <- CreateChromatinAssay(
    counts = counts$Peaks,
    sep = c(":", "-"),
    fragments = fragpath,
    annotation = annotation
  )

  #QC cells
  DefaultAssay(dat) <- "ATAC"

  dat <- NucleosomeSignal(dat)
  dat <- TSSEnrichment(dat)
  dat<-AddMetaData(dat,metadata=metadata_cellranger)

  plt<-VlnPlot(
    object = dat,
    features = c("nCount_RNA", "nCount_ATAC", "TSS.enrichment", "nucleosome_signal"),
    ncol = 4,
    pt.size = 0
  )
  ggsave(plt,file=paste0(i,".qc.pdf"))
  system(paste0("slack -F ",i,".qc.pdf ryan_todo"))
  saveRDS(dat,file=paste0(i,".SeuratObject.rds"))
}

#generate all seurat objects
lapply(c("YW_si_Control_plus_E2","YW_si_Control_plus_Veh","YW_si_KLF4_plus_E2","YW_si_KLF4_plus_Veh"),setupseurat)

#combine seurat objects and filter
yw1<-readRDS("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/YW_si_Control_plus_E2/outs/YW_si_Control_plus_E2.SeuratObject.rds") ; yw1$sample<-"YW_si_Control_plus_E2"
yw2<-readRDS("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/YW_si_Control_plus_Veh/outs/YW_si_Control_plus_Veh.SeuratObject.rds") ; yw2$sample<-"YW_si_Control_plus_Veh"
yw3<-readRDS("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/YW_si_KLF4_plus_E2/outs/YW_si_KLF4_plus_E2.SeuratObject.rds") ; yw3$sample<-"YW_si_KLF4_plus_E2"
yw4<-readRDS("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/YW_si_KLF4_plus_Veh/outs/YW_si_KLF4_plus_Veh.SeuratObject.rds") ; yw4$sample<-"YW_si_KLF4_plus_Veh"
dat <- merge(yw1, y = c(yw2,yw3,yw4), add.cell.ids = c("YW_si_Control_plus_E2","YW_si_Control_plus_Veh","YW_si_KLF4_plus_E2","YW_si_KLF4_plus_Veh"), project = "YW")
saveRDS(dat,file="/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/yw_merged.SeuratObject.rds")

dat<-readRDS("yw_merged.SeuratObject.rds") #ensure it reads in properly

####FILTERING AND QC####
  dat
  table(dat$sample)
  # YW_si_Control_plus_E2 YW_si_Control_plus_Veh     YW_si_KLF4_plus_E2
  #                  7357                   6435                   7727
  #   YW_si_KLF4_plus_Veh
  #                  6785
  # filter out low quality cells
  dat <- subset(
    x = dat,
    subset = nCount_ATAC < 100000 &
      nCount_RNA < 25000 &
      nCount_ATAC > 500 &
      nCount_RNA > 500 &
      nucleosome_signal < 2 &
      TSS.enrichment > 1
  )
  dat
  table(dat$sample)
  # YW_si_Control_plus_E2 YW_si_Control_plus_Veh     YW_si_KLF4_plus_E2
  #                  7030                   6023                   7552
  #   YW_si_KLF4_plus_Veh
  #                  6350
  #Not much filtered out, all passing QC

####INITIAL PROCESSING####
  # call peaks using MACS2
  peaks <- CallPeaks(dat, macs2.path = "/home/groups/CEDAR/mulqueen/src/miniconda3/bin/macs2")

  # remove peaks on nonstandard chromosomes and in genomic blacklist regions
  peaks <- keepStandardChromosomes(peaks, pruning.mode = "coarse")
  peaks <- subsetByOverlaps(x = peaks, ranges = blacklist_hg38_unified, invert = TRUE)

  # quantify counts in each peak
  macs2_counts <- FeatureMatrix(
    fragments = Fragments(dat),
    features = peaks,
    cells = colnames(dat)
  )

  # create a new assay using the MACS2 peak set and add it to the Seurat object
  dat[["peaks"]] <- CreateChromatinAssay(
    counts = macs2_counts,
    fragments = dat@assays$ATAC@fragments,
    annotation = annotation
  )

#RNA Processing
  DefaultAssay(dat) <- "RNA"
  dat <- SCTransform(dat)
  dat <- RunPCA(dat)
  dat<- RunUMAP(
    object = dat,
    reduction.name="rna_umap",
    reduction="pca",
    assay = "SCT",
    verbose = TRUE,
    dims=1:50
  )
  p1<-DimPlot(dat,reduction="rna_umap",group.by="sample")+ggtitle("RNA UMAP")

#DNA Accessibility processing
  DefaultAssay(dat) <- "peaks"
  dat <- FindTopFeatures(dat, min.cutoff = 5)
  dat <- RunTFIDF(dat)
  dat <- RunSVD(dat)
  dat<- RunUMAP(
    object = dat,
    reduction.name="atac_umap",
    reduction="lsi",
    assay = "peaks",
    verbose = TRUE,
    dims=2:40
  )
  p2<-DimPlot(dat,reduction="atac_umap",group.by="sample")+ggtitle("ATAC UMAP")


# build a joint neighbor graph using both assays
  dat <- FindMultiModalNeighbors(
    object = dat,
    reduction.list = list("pca", "lsi"), 
    dims.list = list(1:50, 2:40),
    modality.weight.name = "RNA.weight",
    verbose = TRUE
  )

# build a joint UMAP visualization
  dat <- RunUMAP(
    object = dat,
    nn.name = "weighted.nn",
    reduction.name="multimodal_umap",
    assay = "RNA",
    verbose = TRUE
  )
  p3<-DimPlot(dat,reduction="multimodal_umap",group.by="sample")+ggtitle("Multimodal UMAP")

#Cluster on multimodal graph
dat <- FindClusters(dat, resolution = 0.8, verbose = FALSE,graph="wknn")
p4<-DimPlot(dat,reduction="multimodal_umap",group.by="seurat_clusters")+ggtitle("Multimodal UMAP Clusters")

#Finally Plot results
plt<-(p1 | p2)/(p3 | p4)
  ggsave(plt,file="YW_multimodal.umap.pdf",width=15)
  system("slack -F YW_multimodal.umap.pdf ryan_todo")
  saveRDS(dat,file="yw_merged.SeuratObject.rds")



plt<-VlnPlot(
    object = dat,
    features = c("nCount_RNA", "nCount_ATAC", "TSS.enrichment", "nucleosome_signal"),
    ncol = 4,
    pt.size = 0,
    group.by="sample"
  )
  ggsave(plt,file="yw_merged.qc.pdf")
  system("slack -F yw_merged.qc.pdf ryan_todo")

####Run InferCNV
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
library(GenomeInfoDb)
set.seed(1234)
library(stringr)
library(ggplot2)
library(infercnv)
library(ComplexHeatmap)
library(ggdendro)
library(dendextend)
library(circlize)
library(dendsort)
library(patchwork)


setwd("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/")

####RUNNING INFERCNV#####
#https://bioconductor.org/packages/devel/bioc/manuals/infercnv/man/infercnv.pdf
dat<-readRDS("yw_merged.SeuratObject.rds") #ensure it reads in properly

DefaultAssay(dat)<-"RNA"
#write out raw counts matrix
  counts=as.matrix(dat@assays$RNA@counts[,colnames(dat)])
  write.table(counts,file="YW_inferCNV.counts.txt",sep="\t",col.names=T,row.names=T,quote=F)
#write out cell annotation
  cell_annotation=as.data.frame(cbind(row.names(dat@meta.data),dat@meta.data["sample"]))
  write.table(cell_annotation,file="YW_inferCNV.annotation.txt",sep="\t",col.names=F,row.names=F,quote=F)

# get gene annotations for hg38
  annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
  ucsc.levels <- str_replace(string=paste("chr",seqlevels(annotation),sep=""), pattern="chrMT", replacement="chrM")
  seqlevels(annotation) <- ucsc.levels #standard seq level change threw error, using a string replace instead

#write out gene order list
  gene_order<-annotation[!duplicated(annotation$gene_name),]
  gene_order<-as.data.frame(gene_order[gene_order$gene_name %in% row.names(dat),])
  gene_order<-gene_order[c("gene_name","seqnames","start","end")]
  chrorder<-paste0("chr",c(1:22,"X","Y","M"))
  gene_order$seqnames<-factor(gene_order$seqnames,levels=chrorder) # set chr order
  gene_order<-with(gene_order, gene_order[order(seqnames, start),]) #order by chr and start position
  write.table(gene_order,file="YW_inferCNV.gene_order.txt",sep="\t",col.names=F,row.names=F,quote=F)

#run infercnv
  infercnv_obj = CreateInfercnvObject(raw_counts_matrix="YW_inferCNV.counts.txt",
                                      annotations_file="YW_inferCNV.annotation.txt",
                                      delim="\t",
                                      gene_order_file="YW_inferCNV.gene_order.txt",
                                      ref_group_names=NULL)

saveRDS(infercnv_obj,file="YW_inferCNV.Rds")

#read in where inferCNV keeps hanging (step 14)
infercnv_obj<-readRDS("YW_inferCNV.Rds")

  infercnv_obj = infercnv::run(infercnv_obj,
                               cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics
                               out_dir="./YW_inferCNV", 
                               denoise=TRUE,
                               HMM=FALSE,
                               num_threads=10,
                               cluster_references=FALSE,
                               resume_mode=FALSE
                               )
                              #other options (removed for now)
                                #k_obs_groups=2,
                               #cluster_by_groups=FALSE, 
                               #HMM_type="i6",
                               #
                               #output_format="pdf")


saveRDS(infercnv_obj,file="YW_inferCNV.Rds")

#Run function keeps hanging on step 14, just going to cluster it myself.

infercnv_obj<-readRDS("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/YW_inferCNV/14_invert_log_transform.infercnv_obj")
mat<-as.data.frame(t(infercnv_obj@expr.data))
mat<-mat[,row.names(infercnv_obj@gene_order[infercnv_obj@gene_order$chr=="chr17",])] #limit to chr 17 
mat<-mat[sample(1:nrow(mat),5000),] #limit to 5000 cells, which will be enough to see cluster bias
plt<-Heatmap(mat,
  column_order=1:ncol(mat),
  show_row_names=FALSE,
  show_column_names=FALSE,
  km=2)

pdf("YW.infercnv.heatmap.pdf")
plt<-draw(plt)
dev.off()
system("slack -F YW.infercnv.heatmap.pdf ryan_todo")

celltypes<-row_order(plt) #get row order from heatmap, this includes the slice
metadata_cluster<-as.data.frame(cbind(cellID=rownames(mat),celltype=0))
metadata_cluster[unlist(celltypes[[1]]),]$celltype<-"1"
metadata_cluster[unlist(celltypes[[2]]),]$celltype<-"2"
table(metadata_cluster$celltype)
row.names(metadata_cluster)<-metadata_cluster$cellID
metadata_cluster$cnv_profile<-ifelse(metadata_cluster$celltype==1,"MCF7","T47D") #this is due to the mcf7 amplification we have seen in other data

cell_line<-setNames(metadata_cluster$cnv_profile,row.names(metadata_cluster))
dat<-AddMetaData(dat,metadata=cell_line,col.name="cnv_profile")#assign slices to cell names

#cluster with low resolution on peaks dataset
dat$peaks_cluster<-ifelse(as.numeric(dat@reductions$atac_umap@cell.embeddings[,1])<0,"T47D","MCF7") #it is clear from the atac data that one cluster is T47D and one is MCF7 based on chr17 profiles
table(dat$peaks_cluster)

plt1<-DimPlot(dat,split.by="cnv_profile",group.by="cnv_profile",na.value=NA)+ggtitle("CNV Profile Split")
plt2<-DimPlot(dat,group.by="sample")+ggtitle("ATAC Samples")
plt3<-DimPlot(dat,group.by="peaks_cluster")+ggtitle("Final Cell Type Assignment")

ggsave(plt1/(plt2|plt3),file="YW_celltype.pdf",width=10)
system("slack -F YW_celltype.pdf ryan_todo")

#finally save the results, and separate data files for MCF7 and T47D
saveRDS(dat,file="yw_merged.SeuratObject.rds")
dat_mcf7<-subset(dat,peaks_cluster=="MCF7")
dat_t47d<-subset(dat,peaks_cluster=="T47D")
saveRDS(dat_mcf7,file="yw_mcf7.SeuratObject.rds")
saveRDS(dat_t47d,file="yw_t47d.SeuratObject.rds")

####Run InferCNV
library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
library(GenomeInfoDb)
set.seed(1234)
library(stringr)
library(ggplot2)
library(infercnv)
library(ComplexHeatmap)
library(ggdendro)
library(dendextend)
library(circlize)
library(dendsort)
library(patchwork)


setwd("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/")


####RUNNING INFERCNV#####
#https://bioconductor.org/packages/devel/bioc/manuals/infercnv/man/infercnv.pdf
dat<-readRDS("yw_merged.SeuratObject.rds") #ensure it reads in properly

DefaultAssay(dat)<-"RNA"
#write out raw counts matrix
  counts=as.matrix(dat@assays$RNA@counts[,colnames(dat)])
  #limit to just control samples
  counts<-counts[,startsWith(colnames(counts),"YW_si_Control")] 
  write.table(counts,file="YW_inferCNV.control.counts.txt",sep="\t",col.names=T,row.names=T,quote=F)
#write out cell annotation
  cell_annotation=as.data.frame(cbind(row.names(dat@meta.data),dat@meta.data["sample"]))
  cell_annotation<-cell_annotation[startsWith(row.names(cell_annotation),"YW_si_Control"),] 
  write.table(cell_annotation,file="YW_inferCNV.annotation.control.txt",sep="\t",col.names=F,row.names=F,quote=F)

# get gene annotations for hg38
  annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
  ucsc.levels <- str_replace(string=paste("chr",seqlevels(annotation),sep=""), pattern="chrMT", replacement="chrM")
  seqlevels(annotation) <- ucsc.levels #standard seq level change threw error, using a string replace instead

#write out gene order list
  gene_order<-annotation[!duplicated(annotation$gene_name),]
  gene_order<-as.data.frame(gene_order[gene_order$gene_name %in% row.names(dat),])
  gene_order<-gene_order[c("gene_name","seqnames","start","end")]
  chrorder<-paste0("chr",c(1:22,"X","Y","M"))
  gene_order$seqnames<-factor(gene_order$seqnames,levels=chrorder) # set chr order
  gene_order<-with(gene_order, gene_order[order(seqnames, start),]) #order by chr and start position
  write.table(gene_order,file="YW_inferCNV.gene_order.txt",sep="\t",col.names=F,row.names=F,quote=F)

#run infercnv
  infercnv_obj = CreateInfercnvObject(raw_counts_matrix="YW_inferCNV.control.counts.txt",
                                      annotations_file="YW_inferCNV.annotation.control.txt",
                                      delim="\t",
                                      gene_order_file="YW_inferCNV.gene_order.txt",
                                      ref_group_names=NULL)

saveRDS(infercnv_obj,file="YW_inferCNV.control.Rds")

#read in where inferCNV keeps hanging (step 14)
infercnv_obj<-readRDS("YW_inferCNV.control.Rds")

  infercnv_obj = infercnv::run(infercnv_obj,
                               cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics
                               out_dir="./YW_inferCNV.control", 
                               denoise=TRUE,
                               HMM=TRUE,
                               num_threads=10,
                               cluster_references=FALSE,
                               resume_mode=FALSE,
                               k_obs_groups=2,
                               output_format="pdf"
                               )



saveRDS(infercnv_obj,file="YW_inferCNV.control.Rds")
infercnv_obj<-readRDS("YW_inferCNV.control.Rds")

dat_control<-subset(dat,sample %in% c("YW_si_Control_plus_E2","YW_si_Control_plus_Veh"))
cluster_assignment<-read.table("./YW_inferCNV.control/infercnv.observation_groupings.txt")

cluster_assignment$cnv_profile<-ifelse(cluster_assignment$Dendrogram.Group==1,"MCF7","T47D") #this is due to the mcf7 amplification on chr17 and 20 we have seen in other data

cell_line<-setNames(cluster_assignment$cnv_profile,row.names(cluster_assignment))
dat_control<-AddMetaData(dat_control,metadata=cell_line,col.name="cnv_profile")#assign slices to cell names

#cluster with low resolution on peaks dataset
dat_control$peaks_cluster<-ifelse(as.numeric(dat_control@reductions$atac_umap@cell.embeddings[,1])<0,"T47D","MCF7") #it is clear from the atac data that one cluster is T47D and one is MCF7 based on chr17 profiles
table(dat_control$peaks_cluster)
table(paste(dat_control$peaks_cluster, dat_control$cnv_profile))

#MCF7 MCF7 MCF7 T47D T47D MCF7 T47D T47D
#     6262       176        26      6589
plt1<-DimPlot(dat_control,split.by="cnv_profile",group.by="cnv_profile",na.value=NA)+ggtitle("CNV Profile Split")
plt2<-DimPlot(dat_control,group.by="sample")+ggtitle("ATAC Samples")
plt3<-DimPlot(dat_control,group.by="peaks_cluster")+ggtitle("Final Cell Type Assignment")

ggsave(plt1/(plt2|plt3),file="YW_control.celltype.pdf",width=10)
system("slack -F YW_control.celltype.pdf ryan_todo")

#######CISTOPIC PROCESSING PER CELL LINE#################
setwd("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/")
library(Signac)
library(Seurat)
library(SeuratWrappers)
library(cisTopic)
library(patchwork)
set.seed(1234)
library(org.Hs.eg.db)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(AUCell)
library(rtracklayer)
library(parallel)
library(ComplexHeatmap)

cistopic_generation<-function(x,name_out,outdir="/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi"){
  wd<-outdir
  outname<-name_out
  atac_sub<-x
  cistopic_counts_frmt<-atac_sub$peaks@counts
  row.names(cistopic_counts_frmt)<-sub("-", ":", row.names(cistopic_counts_frmt))
  sub_cistopic<-cisTopic::createcisTopicObject(cistopic_counts_frmt)
  print("made cistopic object")
  sub_cistopic_models<-cisTopic::runWarpLDAModels(sub_cistopic,topic=c(10:30),nCores=5,addModels=FALSE)
  sub_cistopic_models<-addCellMetadata(sub_cistopic_models, cell.data =x@meta.data)
  pdf(paste0(wd,"/",outname,"_model_selection.pdf"))
  par(mfrow=c(3,3))
  sub_cistopic_models<- selectModel(sub_cistopic_models, type='derivative')
  dev.off()
  system(paste0("slack -F ",paste0(wd,"/",outname,"_model_selection.pdf")," ryan_todo"))
  
  saveRDS(sub_cistopic_models,file=paste0(wd,"/",outname,".CisTopicObject.Rds"))
  sub_cistopic_models<-readRDS(file=paste0(wd,"/",outname,".CisTopicObject.Rds"))
  print("finshed running cistopic")

  #Add cell embeddings into seurat
  cell_embeddings<-as.data.frame(sub_cistopic_models@selected.model$document_expects)
  colnames(cell_embeddings)<-sub_cistopic_models@cell.names
  n_topics<-nrow(cell_embeddings)
  row.names(cell_embeddings)<-paste0("topic_",1:n_topics)
  cell_embeddings<-as.data.frame(t(cell_embeddings))

  #Add feature loadings into seurat
  feature_loadings<-as.data.frame(sub_cistopic_models@selected.model$topics)
  row.names(feature_loadings)<-paste0("topic_",1:n_topics)
  feature_loadings<-as.data.frame(t(feature_loadings))

  #combined cistopic results (cistopic loadings and umap with seurat object)
  cistopic_obj<-CreateDimReducObject(embeddings=as.matrix(cell_embeddings),loadings=as.matrix(feature_loadings),assay="peaks",key="topic_")
  atac_sub@reductions$cistopic<-cistopic_obj
  n_topics<-ncol(Embeddings(atac_sub,reduction="cistopic")) #add scaling for ncount peaks somewhere in here
  atac_sub<-RunUMAP(atac_sub,reduction="cistopic",dims=1:n_topics)
  atac_sub <- FindNeighbors(object = atac_sub, reduction = 'cistopic', dims = 1:n_topics ) 
  atac_sub <- FindClusters(object = atac_sub, verbose = TRUE, graph.name="peaks_snn", resolution=0.2 ) 
  plt1<-DimPlot(atac_sub,reduction="umap",group.by=c("sample","seurat_clusters"))
  plt2<-FeaturePlot(atac_sub,reduction="umap",features=c("nucleosome_signal","TSS.enrichment","nCount_peaks","nFeature_peaks"))
  pdf(paste0(wd,"/",outname,".umap.pdf"),width=10)
  print(plt1)
  print(plt2)
  dev.off()
  system(paste0("slack -F ",paste0(wd,"/",outname,".umap.pdf")," ryan_todo"))
  saveRDS(atac_sub,paste0(wd,"/",outname,".SeuratObject.rds"))
  }


#Running cistopic on just control experiment (FOR TITAN PAPER)
dat_mcf7_full<-readRDS(file="yw_mcf7.SeuratObject.rds")
dat_mcf7<-subset(dat_mcf7_full, sample %in% c("YW_si_Control_plus_E2","YW_si_Control_plus_Veh"))
cistopic_generation(x=dat_mcf7,name_out="yw_mcf7.control")

dat_t47d_full<-readRDS(file="yw_t47d.SeuratObject.rds")
dat_t47d<-subset(dat_t47d_full, sample %in% c("YW_si_Control_plus_E2","YW_si_Control_plus_Veh"))
cistopic_generation(x=dat_t47d,name_out="yw_t47d.control")

#read in cistrome data for topic analysis
cistrome_db<-read.csv("/home/groups/CEDAR/mulqueen/ref/cistrome/human_factor_full_QC.txt",sep="\t") #has information on each download peaks files
cistrome_db<-cistrome_db[cistrome_db$Cell_line %in% c("MCF-7","T47D"),] #limit to our cell types
cistrome_db<-cistrome_db[cistrome_db$PeaksFoldChangeAbove10>1000,] #set lower limit for peaks
cistrome_dir="/home/groups/CEDAR/mulqueen/ref/cistrome/human_factor" #has individual peak files
ChIP_Seq_signatures <- paste(cistrome_dir, list.files(cistrome_dir), sep='/') #grab all files in directory
signature_dcids<-unlist(lapply(strsplit(unlist(lapply(strsplit(ChIP_Seq_signatures,"human_factor/"),"[",2)),"_sort"),"[",1)) #grab dcID from file name
ChIP_Seq_signatures <- ChIP_Seq_signatures[signature_dcids %in% cistrome_db$DCid] #limit to those in cistrome_db after filter
ChIP_Seq_signatures<-ChIP_Seq_signatures[order(as.numeric(unlist(lapply(strsplit(unlist(lapply(strsplit(ChIP_Seq_signatures,"human_factor/"),"[",2)),"_sort"),"[",1))))] #sort by numeric, like db
cistrome_db$label<-paste(cistrome_db$Cell_line,cistrome_db$Factor,cistrome_db$DCid,sep="_")
cistrome_db$file<-ChIP_Seq_signatures

hg38_chip<-paste(cistrome_dir, list.files(cistrome_dir,pattern=".bed"), sep='/') #grab all files in directory
hg38_chip<-hg38_chip[which(hg38_chip %in% cistrome_db$file)] #filter to our files of interest
labels<-cistrome_db[cistrome_db$file %in% hg38_chip,]$label #grab meaningful file labels



#Take cistopic data per cell line and perform bed file overlap using cistrome cell line specific ChIP-seq.

#run through further processing
#this function will perform more analysis on topics in the selected model, looking at 
#1. cell topic weights
#2. cistrome chipseq peak overlaps per topic
#3. topic annotations (if peaks fall primarily in promoter regions etc)

cistopic_processing<-function(x,y,outname){
  obj<-x
  cisTopicObject<-y
  cisTopicObject<-addCellMetadata(cisTopicObject, cell.data =obj@meta.data)
  #heatmap by topic
  cisTopicObject <- runUmap(cisTopicObject, target='cell')
  pdf(paste(outname,"celltopic_heatmap.pdf",sep="."))
  print(cellTopicHeatmap(cisTopicObject, method='Probability', colorBy=c("sample","seurat_clusters")))
  dev.off()
  system(paste("slack -F", paste(outname,"celltopic_heatmap.pdf",sep="."), "ryan_todo",sep=" "))

  #region score association
  pred.matrix <- predictiveDistribution(cisTopicObject)
  cisTopicObject <- getSignaturesRegions(cisTopicObject, hg38_chip, labels=labels, minOverlap = 0.01) #output new bed files to be read in and processed
  aucellRankings <- AUCell_buildRankings(pred.matrix, plot=FALSE, verbose=FALSE)   # Compute cell rankings
  cisTopicObject <- signatureCellEnrichment(cisTopicObject, aucellRankings, selected.signatures='all', aucMaxRank = 0.1*nrow(aucellRankings),plot=FALSE)   # Check signature enrichment in cells

  cisTopicObject <- getRegionsScores(cisTopicObject, method='NormTop', scale=TRUE)
  cisTopicObject <- binarizecisTopics(cisTopicObject, thrP=0.999, plot=TRUE)

  txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
  cisTopicObject <- annotateRegions(cisTopicObject, txdb=TxDb.Hsapiens.UCSC.hg38.knownGene, annoDb='org.Hs.eg.db')
  cisTopicObject <- runtSNE(cisTopicObject, target='region', perplexity=200, check_duplicates=FALSE) #cluster by regions
  saveRDS(cisTopicObject,paste0(outname,".CisTopicObject.Rds"))

  #save cistopic region data
  outdat<-lapply(cisTopicObject@binarized.cisTopics, function(x) {
    temp<-merge(x,cisTopicObject@region.data,by="row.names")
    temp$topic_assigned<-names(x)
    return(temp)})

  outdat<-data.table::rbindlist(outdat,use.names=FALSE)
  write.table(outdat,file=paste0(outname,".topic_region_assignment.txt"),col.names=T,row.names=T,quote=F,sep="\t")
  system(paste("slack -F ",paste0(outname,".topic_region_assignment.txt")," ryan_todo" ))
  pdf(paste0(outname,".topic_annotations.pdf"))
  par(mfrow=c(1,1))
  signaturesHeatmap(cisTopicObject, selected.signatures = 'annotation')
  plotFeatures(cisTopicObject, method='tSNE', target='region', topic_contr=NULL, colorBy=c('annotation'), cex.legend = 0.8, factor.max=.75, dim=2, legend=TRUE, intervals=20)
  dev.off()
  system(paste("slack -F ",paste0(outname,".topic_annotations.pdf")," ryan_todo" ))

  pdf(paste0(outname,".topic_scores.pdf"))
  par(mfrow=c(4,4))
  plotFeatures(cisTopicObject, method='Umap', target='cell', topic_contr='Probability', colorBy=NULL, cex.legend = 0.8, factor.max=.75, dim=2, legend=TRUE)
  dev.off()
  system(paste("slack -F ",paste0(outname,".topic_scores.pdf")," ryan_todo" ))

  obj[['cistrome_score']] <- CreateAssayObject(
    data = t(cisTopicObject@cell.data[(row.names(cisTopicObject@cell.data) %in% row.names(obj@meta.data)),
      (startsWith(colnames(cisTopicObject@cell.data),"MCF") | startsWith(colnames(cisTopicObject@cell.data),"T47D"))]
    ))
  saveRDS(obj,file=paste0(outname,".SeuratObject.rds"))

}



#add cistopic function to all cells
dat_mcf7<-readRDS(file="yw_mcf7.control.SeuratObject.rds")
dat_mcf7_cistopic<-readRDS(file="yw_mcf7.control.CisTopicObject.Rds")
cistopic_processing(x=dat_mcf7,y=dat_mcf7_cistopic,outname="yw_mcf7.control")
#x=dat_mcf7;y=dat_mcf7_cistopic;outname="yw_mcf7.control"

dat_t47d<-readRDS(file="yw_t47d.control.SeuratObject.rds") 
dat_t47d_cistopic<-readRDS(file="yw_t47d.control.CisTopicObject.Rds") 
cistopic_processing(x=dat_t47d,y=dat_t47d_cistopic,outname="yw_t47d.control")
#x=dat_t47d;y=dat_t47d_cistopic;outname="yw_t47d.control"

  library(Signac)
  library(Seurat)
  library(JASPAR2020)
  library(TFBSTools)
  library(BSgenome.Hsapiens.UCSC.hg38)
  library(patchwork)
  set.seed(1234)
library(BiocParallel)
register(SerialParam()) #using single core mode

setwd("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/")

# Get a list of motif position frequency matrices from the JASPAR database
pfm <- getMatrixSet(
  x = JASPAR2020,
  opts = list(species =9606, all_versions = FALSE))

# Scan the DNA sequence of each peak for the presence of each motif, using orgo_atac for all objects (shared peaks)
mcf7<-readRDS(file="yw_mcf7.control.SeuratObject.rds")
 peaks<-granges(mcf7[["peaks"]])
  peaks<-peaks[seqnames(peaks) %in% c(paste0("chr",c(1:22,"X","Y"))),]
    motif.matrix.hg38 <- CreateMotifMatrix(features = peaks, pwm = pfm, genome = BSgenome.Hsapiens.UCSC.hg38, use.counts = FALSE)
    motif.hg38 <- CreateMotifObject(data = motif.matrix.hg38, pwm = pfm)
  mcf7 <- SetAssayData(object = mcf7, assay = 'peaks', slot = 'motifs', new.data = motif.hg38)
  mcf7 <- RegionStats(object = mcf7, genome = BSgenome.Hsapiens.UCSC.hg38,assay="peaks")
  mcf7 <- RunChromVAR( object = mcf7,genome = BSgenome.Hsapiens.UCSC.hg38,assay="peaks")
  saveRDS(mcf7,file="yw_mcf7.control.SeuratObject.rds")


#now run t47d
t47d<-readRDS(file="yw_t47d.control.SeuratObject.rds")
 peaks<-granges(t47d[["peaks"]])
  peaks<-peaks[seqnames(peaks) %in% c(paste0("chr",c(1:22,"X","Y"))),]
    motif.matrix.hg38 <- CreateMotifMatrix(features = peaks, pwm = pfm, genome = BSgenome.Hsapiens.UCSC.hg38, use.counts = FALSE)
    motif.hg38 <- CreateMotifObject(data = motif.matrix.hg38, pwm = pfm)
  t47d <- SetAssayData(object = t47d, assay = 'peaks', slot = 'motifs', new.data = motif.hg38)
  t47d <- RegionStats(object = t47d, genome = BSgenome.Hsapiens.UCSC.hg38,assay="peaks")
  t47d <- RunChromVAR( object = t47d,genome = BSgenome.Hsapiens.UCSC.hg38,assay="peaks")
  saveRDS(t47d,file="yw_t47d.control.SeuratObject.rds")

  library(Signac)
  library(Seurat)
  library(SeuratWrappers)
  library(ggplot2)
  library(patchwork)
  #library(monocle3)
  library(cicero)
  library(SeuratObjects)
  library(EnsDb.Hsapiens.v86)
  setwd("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/")
  mcf7<-readRDS(file="yw_mcf7.control.SeuratObject.rds")
  t47d<-readRDS(file="yw_t47d.control.SeuratObject.rds")

  #Cicero processing function
  cicero_processing<-function(object_input,prefix){

      #Generate CDS format from Seurat object
      atac.cds <- as.CellDataSet(object_input,assay="peaks",reduction="umap")

      # convert to CellDataSet format and make the cicero object
      print("Making Cicero format CDS file")
      atac.cicero <- make_cicero_cds(atac.cds, reduced_coordinates = atac.cds@reducedDimS)
      saveRDS(atac.cicero,paste(prefix,"atac_cicero_cds.Rds",sep="_"))
      
      # extract gene annotations from EnsDb
      annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
      seqlevels(annotations)<-paste0("chr",seqlevels(annotations))

      # change to UCSC style since the data was mapped to hg38
      #seqlevelsStyle(annotations) <- 'UCSC'
      genome(annotations) <- "hg38"

      genome <- annotations@seqinfo # get the chromosome sizes from the Seurat object
      genome.df <- data.frame("chr" = genome@seqnames, "length" = genome@seqlengths) # convert chromosome sizes to a dataframe
      
      print("Running Cicero to generate connections.")
      conns <- run_cicero(atac.cicero, genomic_coords = genome.df) # run cicero
      saveRDS(conns,paste(prefix,"atac_cicero_conns.Rds",sep="_"))
      
      print("Generating CCANs")
      ccans <- generate_ccans(conns) # generate ccans
      saveRDS(ccans,paste(prefix,"atac_cicero_ccans.Rds",sep="_"))
      
      print("Adding CCAN links into Seurat Object and Returning.")
      links <- ConnectionsToLinks(conns = conns, ccans = ccans) #Add connections back to Seurat object as links
      DefaultAssay(object_input)<-"peaks"
      Links(object_input) <- links
      return(object_input)
  }

t47d<-cicero_processing(object_input=t47d,prefix="yw_t47d.control")
saveRDS(t47d,"yw_t47d.control.SeuratObject.unnormGA.Rds")

mcf7<-cicero_processing(object_input=mcf7,prefix="yw_mcf7.control")
saveRDS(mcf7,"yw_mcf7.control.SeuratObject.unnormGA.Rds")

# generate unnormalized gene activity matrix
  # gene annotation sample
  hg38_annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)

  pos <-as.data.frame(hg38_annotations,row.names=NULL)
  pos$chromosome<-paste0("chr",pos$seqnames)
  pos$gene<-pos$gene_id
  pos <- subset(pos, strand == "+")
  pos <- pos[order(pos$start),] 
  pos <- pos[!duplicated(pos$tx_id),] # remove all but the first exons per transcript
  pos$end <- pos$start + 1 # make a 1 base pair marker of the TSS

  neg <-as.data.frame(hg38_annotations,row.names=NULL)
  neg$chromosome<-paste0("chr",neg$seqnames)
  neg$gene<-neg$gene_id
  neg <- subset(neg, strand == "-")
  neg <- neg[order(neg$start,decreasing=TRUE),] 
  neg <- neg[!duplicated(neg$tx_id),] # remove all but the first exons per transcript
  neg$end <- neg$end + 1 # make a 1 base pair marker of the TSS

  gene_annotation<- rbind(pos, neg)
  gene_annotation <- gene_annotation[,c("chromosome","start","end","gene_name")] # Make a subset of the TSS annotation columns containing just the coordinates and the gene name
  names(gene_annotation)[4] <- "gene" # Rename the gene symbol column to "gene"

  geneactivity_processing<-function(cds_input,conns_input,prefix){
      atac.cds<- annotate_cds_by_site(cds_input, gene_annotation)
      unnorm_ga <- build_gene_activity_matrix(atac.cds, conns_input)
      saveRDS(unnorm_ga,paste(prefix,"unnorm_GA.Rds",sep="."))
  }

#t47d
   conns<-as.data.frame(readRDS("yw_t47d.control_atac_cicero_conns.Rds"))
   obj.cicero<-readRDS("yw_t47d.control_atac_cicero_cds.Rds")
   geneactivity_processing(cds_input=as.CellDataSet(t47d,assay="peaks",reduction="umap"),conns_input=conns,prefix="yw_t47d.control")

#mcf7
   conns<-as.data.frame(readRDS("yw_mcf7.control_atac_cicero_conns.Rds"))
   obj.cicero<-readRDS("yw_mcf7.control_atac_cicero_cds.Rds")
   geneactivity_processing(cds_input=as.CellDataSet(mcf7,assay="peaks",reduction="umap"),conns_input=conns,prefix="yw_mcf7.control")

#Read in unnormalized GA
cicero_gene_activities<-readRDS("yw_t47d.control.unnorm_GA.Rds")
t47d[['GeneActivity']]<- CreateAssayObject(counts = cicero_gene_activities) 
# normalize
   t47d <- NormalizeData(
     object = t47d,
     assay = 'GeneActivity',
     normalization.method = 'LogNormalize',
     scale.factor = median(t47d$nCount_GeneActivity)
   )
   saveRDS(t47d,"yw_t47d.control.SeuratObject.rds")


#Read in unnormalized GA
cicero_gene_activities<-readRDS("yw_mcf7.control.unnorm_GA.Rds")
mcf7[['GeneActivity']]<- CreateAssayObject(counts = cicero_gene_activities) 
# normalize
   mcf7 <- NormalizeData(
     object = mcf7,
     assay = 'GeneActivity',
     normalization.method = 'LogNormalize',
     scale.factor = median(mcf7$nCount_GeneActivity)
   )
   saveRDS(mcf7,"yw_mcf7.control.SeuratObject.rds")

library(Signac)
library(Seurat)
library(ComplexHeatmap)
library(JASPAR2020)
library(TFBSTools)
library(grid)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(parallel)
library(BSgenome.Hsapiens.UCSC.hg38)
library(patchwork)
library(circlize)
library(viridis)
library(reshape2)


setwd("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/")
mcf7_titan<-readRDS("/home/groups/CEDAR/doe/projects/Lydia_ATACRNA/justMCF7_wTopics.rds")
t47d_titan<-readRDS("/home/groups/CEDAR/doe/projects/Lydia_ATACRNA/justT47D_wTopics.rds")

mcf7<-readRDS(file="yw_mcf7.control.SeuratObject.rds")
t47d<-readRDS(file="yw_t47d.control.SeuratObject.rds")

#Add titan reduction to full seurat object
mcf7[["titan_lda"]]<-mcf7_titan[["lda"]]
t47d[["titan_lda"]]<-t47d_titan[["lda"]]

saveRDS(mcf7,file="yw_mcf7.control.SeuratObject.rds")
saveRDS(t47d,file="yw_t47d.control.SeuratObject.rds")

library(Signac)
library(Seurat)
set.seed(1234)
library(stringr)
library(ggplot2)

setwd("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/")
mcf7<-readRDS(file="yw_mcf7.control.SeuratObject.rds")
t47d<-readRDS(file="yw_t47d.control.SeuratObject.rds")

cluster_everything<-function(dat,outname){
  #RNA Processing
  DefaultAssay(dat) <- "RNA"
  dat <- SCTransform(dat)
  dat <- RunPCA(dat)
  dat<- RunUMAP(object = dat, reduction.name="rna_umap", reduction="pca", assay = "SCT", verbose = TRUE, dims=1:50 )
  #DNA Accessibility processing
  DefaultAssay(dat) <- "peaks"
  dat <- FindTopFeatures(dat, min.cutoff = 5)
  dat <- RunTFIDF(dat)
  dat <- RunSVD(dat)
  dat<- RunUMAP(object = dat, reduction.name="atac_umap", reduction="lsi", assay = "peaks", verbose = TRUE, dims=2:40 )
  # build a joint neighbor graph using both assays
  dat <- FindMultiModalNeighbors(object = dat, reduction.list = list("pca", "lsi"), dims.list = list(1:50, 2:40), modality.weight.name = c("SCT.weight","peaks.weight"), verbose = TRUE ) # build a joint neighbor graph using both assays
  dat <- RunUMAP(object = dat, nn.name = "weighted.nn", reduction.name="multimodal_umap", assay = "RNA", verbose = TRUE ) # build a joint UMAP visualization
  
  #run the topics UMAPs too
  dat<- RunUMAP(object = dat, reduction.name="titan_umap", reduction="titan_lda", assay = "peaks", verbose = TRUE, dims=1:ncol(dat@reductions$titan_lda@cell.embeddings) )
  dat<- RunUMAP(object = dat, reduction.name="cistopic_umap", reduction="cistopic", assay = "peaks", verbose = TRUE, dims=1:ncol(dat@reductions$cistopic@cell.embeddings) ) 
  # build a joint neighbor graph using both topic reductions
  dat <- FindMultiModalNeighbors(object = dat, reduction.list = list("titan_lda", "cistopic"), dims.list = list(1:ncol(dat@reductions$titan_lda@cell.embeddings), 1:ncol(dat@reductions$cistopic@cell.embeddings)), modality.weight.name = c("SCT.weight","peaks.weight"), verbose = TRUE )
  dat <- RunUMAP(object = dat, nn.name = "weighted.nn", reduction.name="multimodal_topic_umap", assay = "peaks", verbose = TRUE ) # build a joint UMAP visualization

#Finally Plot results
p1<-DimPlot(dat,reduction="rna_umap",group.by="sample")+ggtitle("RNA UMAP")
p2<-DimPlot(dat,reduction="atac_umap",group.by="sample")+ggtitle("ATAC UMAP")
p3<-DimPlot(dat,reduction="multimodal_umap",group.by="sample")+ggtitle("Multimodal UMAP")

p4<-DimPlot(dat,reduction="titan_umap",group.by="sample")+ggtitle("TITAN UMAP")
p5<-DimPlot(dat,reduction="cistopic_umap",group.by="sample")+ggtitle("CISTOPIC UMAP")
p6<-DimPlot(dat,reduction="multimodal_topic_umap",group.by="sample")+ggtitle("Multimodal Topic UMAP")
plt<-(p1 | p2)/(p3 | p4)/(p5 | p6)

  ggsave(plt,file=paste0(outname,".umap.pdf"),width=10,height=10)
  system(paste0("slack -F ",outname,".umap.pdf ryan_todo"))
return(dat)
}


mcf7<-cluster_everything(dat=mcf7,outname="yw_mcf7.control")
t47d<-cluster_everything(dat=t47d,outname="yw_t47d.control")
saveRDS(mcf7,file="yw_mcf7.control.SeuratObject.rds")
saveRDS(t47d,file="yw_t47d.control.SeuratObject.rds")

library(Signac)
library(Seurat)
library(ComplexHeatmap)
library(JASPAR2020)
library(TFBSTools)
library(grid)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(parallel)
library(BSgenome.Hsapiens.UCSC.hg38)
library(patchwork)
library(circlize)
library(viridis)
library(reshape2)
library(cisTopic)
library(ComplexHeatmap)

setwd("/home/groups/CEDAR/mulqueen/projects/multiome/220111_multi/")

mcf7<-readRDS(file="yw_mcf7.control.SeuratObject.rds")
t47d<-readRDS(file="yw_t47d.control.SeuratObject.rds")

#ESR1 and FOXM1 Topics identified by aysegul
#MCF7: FOXM1=topic4 , ESR1 = topic37
#T47D: FOXM1=topic11 , ESR1 = topic23
mcf7_foxm1_topic="lda_4";mcf7_esr1_topic="lda_37"
t47d_foxm1_topic="lda_11";t47d_esr1_topic="lda_23"

#Function to fit a linear regression to the relevant topics, and then bin data by quantiles
#maybe change binning to be based on ratio?

#Plot cistopic by cell matrix
dat_mcf7_cistopic<-readRDS(file="yw_mcf7.control.CisTopicObject.Rds")
dat_t47d_cistopic<-readRDS(file="yw_t47d.control.CisTopicObject.Rds")


make_heatmap<-function(seurat_object=mcf7,
  object=dat_mcf7_cistopic,
  method='Probability',
  colorBy=c('sample'),
  topic_order=mcf7_foxm1_topic,
  other_topic=mcf7_esr1_topic,
  outname){

    colorPal <- grDevices::colorRampPalette(c("white","red"))
    col_sample<-c("YW_si_Control_plus_Veh"="red","YW_si_Control_plus_E2"="blue")
    foxm1_list<-seurat_object@reductions$titan_lda@cell.embeddings[,topic_order]
    esr1_list<-seurat_object@reductions$titan_lda@cell.embeddings[,other_topic]
    col_foxm1<-colorRamp2(breaks=c(min(foxm1_list),max(foxm1_list)),colors=c("white","blue"))
    col_esr1<-colorRamp2(breaks=c(min(esr1_list),max(esr1_list)),colors=c("white","red"))
    #col_foxm1<-grDevices::colorRampPalette(c("white","black"))
    topic.mat <- modelMatSelection(object, "cell", method)
    rownames(topic.mat) <- paste("Topic", seq(1, nrow(topic.mat)))
    colnames(topic.mat) <- object@cell.names
    object.cell.data <- object@cell.data
    #set up ordering
    object.cell.data$sample_order<-ifelse(object.cell.data$sample=="YW_si_Control_plus_Veh",0,1)
    object.cell.data$foxm1_topic<-foxm1_list
    object.cell.data$esr1_topic<-esr1_list
    col_order<-with(object.cell.data,order(sample_order,foxm1_topic))
    object.cell.data<-object.cell.data[col_order,]
    topic.mat<-topic.mat[,col_order]


        heatmap <- ComplexHeatmap::Heatmap(topic.mat,
            col = colorPal(20), 
            top_annotation = HeatmapAnnotation(treatment = anno_simple(object.cell.data$sample,col=col_sample),
              foxm1_topic=anno_simple(object.cell.data$foxm1_topic,col=col_foxm1),
              esr1_topic=anno_simple(object.cell.data$esr1_topic,col=col_esr1)),
            column_split = object.cell.data[,"sample"],
            column_order=1:ncol(topic.mat),
            name = method,
            show_column_names = FALSE, 
            show_row_names = TRUE,
            column_title = "Topic contribution per cell", 
            column_title_gp = gpar(fontface = "bold"))
        pdf(outname)
        print(heatmap)
        dev.off()
        system(paste0("slack -F ",outname, " ryan_todo"))
}

make_heatmap(object=dat_mcf7_cistopic,
  topic_order=mcf7_foxm1_topic,
  other_topic=mcf7_esr1_topic,
  seurat_object=mcf7,
  outname="yw_mcf7_cellbyTopic_heatmap.pdf")

make_heatmap(object=dat_t47d_cistopic,
  topic_order=t47d_foxm1_topic,
  other_topic=t47d_esr1_topic,
  seurat_object=t47d,
  outname="yw_t47d_cellbyTopic_heatmap.pdf")




scatterplot_linearfit<-function(x,esr1_topic,foxm1_topic,outname){
  titan_dat<-as.data.frame(x@reductions$titan_lda@cell.embeddings) #get titan from y
  cor_mat<-titan_dat[,c(esr1_topic,foxm1_topic)]
  cor_mat$sample<-x$sample
  object<-x

  #generate a linear model
  x=as.numeric(cor_mat[,esr1_topic])
  y=as.numeric(cor_mat[,foxm1_topic])
  model_fit<-summary(lm(y~x,data=cor_mat))
  intercept<-model_fit$coefficients[1,1]
  slope<-model_fit$coefficients[2,1]
  std_err<-model_fit$coefficients[2,2]
  residuals<-model_fit$residuals
  print(model_fit)

  plt<-ggplot()+
    geom_point(aes(x=cor_mat[,esr1_topic],y=cor_mat[,foxm1_topic],color=cor_mat$sample),size=0.1)+
    theme_minimal()+
    geom_abline(intercept = intercept, slope = slope, col = "black")+
    ggtitle(paste(outname,"ESR1 and FOXM1 Topic Correlations"))+
    xlab(paste("ESR1",esr1_topic))+
    ylab(paste("FOXM1",foxm1_topic))+
    theme(axis.text = element_text(size = 6),axis.title=element_text(size=8),axis.ticks=element_blank())+
    scale_x_continuous(expand = expansion(mult = 0.5))+  scale_y_continuous(expand = expansion(mult = 0.5))

    ggsave(plt,file=paste0(outname,".ESR1_FOXM1.scatter.sample.pdf"),height=4,width=7,units="in")
    print(plt)
    dev.off()
    system(paste0("slack -F ",paste0(outname,".ESR1_FOXM1.scatter.sample.pdf")," ryan_todo"))

  #bins defined as top of topic distribution only in treated sample, to be used in later analysis
  cor_mat$topic_bin<-""
  #cor_mat[cor_mat$sample=="YW_si_Control_plus_E2" & cor_mat[,foxm1_topic]>quantile(cor_mat[,foxm1_topic],0.95) & cor_mat[,esr1_topic]<quantile(cor_mat[,esr1_topic],0.75),]$topic_bin<-"FOXM1" #was >0.95 and <0.75 for other, new is 0.75 and 0.5
  #cor_mat[cor_mat$sample=="YW_si_Control_plus_E2" &cor_mat[,esr1_topic]>quantile(cor_mat[,esr1_topic],0.95) & cor_mat[,foxm1_topic]<quantile(cor_mat[,foxm1_topic],0.75),]$topic_bin<-"ESR1" 
  cor_mat[cor_mat[,foxm1_topic]>quantile(cor_mat[,foxm1_topic],0.8) & cor_mat[,esr1_topic]<quantile(cor_mat[,esr1_topic],0.75),]$topic_bin<-"FOXM1" #was >0.95 and <0.75 for other, new is 0.75 and 0.5
  cor_mat[cor_mat[,esr1_topic]>quantile(cor_mat[,esr1_topic],0.8) & cor_mat[,foxm1_topic]<quantile(cor_mat[,foxm1_topic],0.75),]$topic_bin<-"ESR1" 
  print(table(cor_mat$topic_bin))


  plt<-ggplot()+
    geom_point(aes(x=cor_mat[,esr1_topic],y=cor_mat[,foxm1_topic],color=cor_mat$topic_bin),size=0.1)+
    theme_minimal()+
    scale_color_manual(values=c("NA"=NA,"FOXM1"="red","ESR1"="blue"))+
    geom_abline(intercept = intercept, slope = slope, col = "black")+
    ggtitle(paste(outname,"ESR1 and FOXM1 Topic Correlations"))+
    xlab(paste("ESR1",esr1_topic))+
    ylab(paste("FOXM1",foxm1_topic))+
    theme(axis.text = element_text(size = 6),axis.title=element_text(size=8),axis.ticks=element_blank())+
    scale_x_continuous(expand = expansion(mult = 0.5))+  scale_y_continuous(expand = expansion(mult = 0.5))

    ggsave(plt,file=paste0(outname,".ESR1_FOXM1.scatter.binned.pdf"),height=4,width=7,units="in")
    print(plt)
    dev.off()
    system(paste0("slack -F ",paste0(outname,".ESR1_FOXM1.scatter.binned.pdf")," ryan_todo"))

  bin_metadata<-setNames(as.character(cor_mat$topic_bin),rownames(cor_mat))
  object<-AddMetaData(object,metadata=bin_metadata,col.name="topic_bin")
  return(object)
}

mcf7<-scatterplot_linearfit(x=mcf7,outname="yw_mcf7.control.newbin",foxm1_topic=mcf7_foxm1_topic,esr1_topic=mcf7_esr1_topic) #p val 0.2161 #Adjusted R-squared:  8.235e-05
t47d<-scatterplot_linearfit(x=t47d,outname="yw_t47d.control.newbin",foxm1_topic=t47d_foxm1_topic,esr1_topic=t47d_esr1_topic) # pval p-value: < 2.2e-16 #Adjusted R-squared:  0.02529

#save now that there is titan data and its binned
saveRDS(mcf7,file="yw_mcf7.control.SeuratObject.rds")
saveRDS(t47d,file="yw_t47d.control.SeuratObject.rds")

table(mcf7$topic_bin)
table(t47d$topic_bin)
Idents(mcf7)<-mcf7$topic_bin
#      ESR1 FOXM1
# 4159  1152  1127
Idents(t47d)<-t47d$topic_bin
#       ESR1 FOXM1
# 4797   914   904

find_markers_topic_bin<-function(x,assay_name,outname,logfc.threshold_set=0,min.pct_set=0,cistrome_cellline="MCF",latent.vars="NA"){
  if(!(latent.vars=="NA")){
  x_da<-FindMarkers(object=x,ident.1 = "ESR1", ident.2="FOXM1", test.use = 'LR', logfc.threshold=logfc.threshold_set,latent.vars = latent.vars, only.pos=F, assay=assay_name,min.pct=min.pct_set)
  } else {
  x_da<-FindMarkers(object=x,ident.1 = "ESR1", ident.2="FOXM1", test.use = 'LR', logfc.threshold=logfc.threshold_set, only.pos=F, assay=assay_name,min.pct=min.pct_set)
  }
  #if statements to handle formating the different modalities
  if(assay_name=="cistrome_score"){ #if assay is cistrome limit to correct cell line
    x_da<-x_da[startsWith(row.names(x_da),cistrome_cellline),] #subset to same cell line cistrome data
    label_length<-length(strsplit(row.names(x_da),"-")[[1]])
    row.names(x_da)<-paste(unlist(lapply(strsplit(row.names(x_da),"-"),"[",label_length-1)),
                          unlist(lapply(strsplit(row.names(x_da),"-"),"[",label_length)),sep="-")
    }
  if(assay_name=="chromvar"){ #make chromvar names readable
      print("Translating Chromvar Motif Names for Plot")
      x_da$out_name <- unlist(lapply(unlist(lapply(row.names(x_da), function(x) getMatrixByID(JASPAR2020,ID=x))),function(y) name(y)))
  } else{
    x_da$out_name<-row.names(x_da)
  }
  if(assay_name=="peaks"){ #assign peak output to closest features
    closest_genes <- ClosestFeature(x,row.names(x_da))
    x_da <-cbind(x_da ,closest_genes)
  }

  x_da$sig<-ifelse(x_da$p_val_adj<=0.05,"sig","non_sig")
  x_da$label<-""

  if(assay_name!="peaks"){ #ignore adding peaks names (too many), and add top 10 significant names to either side of comparison
    ident_1_labels<- row.names(head(x_da[which((x_da$sig=="sig") & (x_da$avg_log2FC>0)),] ,n=10))#significant, is ident1 specific, is top 10
    ident_2_labels<- row.names(head(x_da[which((x_da$sig=="sig") & (x_da$avg_log2FC<0)),] ,n=10))#significant, is ident1 specific, is top 10
    x_da[c(ident_1_labels,ident_2_labels),]$label<-x_da[c(ident_1_labels,ident_2_labels),]$out_name
        #trim cistrome cell line from name (for readability)
  }
  x_scale<-max(abs(x_da$avg_log2FC))*1.1 #find proper scaling for x axis

  plt<-ggplot(x_da,aes(x=avg_log2FC,y=-log10(p_val_adj),color=sig,label=label,alpha=0.1))+
    geom_point(size=0.5)+
    theme_bw()+
    scale_fill_manual(values=c("non_sig"="#999999", "sig"="#FF0000"))+
    xlim(c(-x_scale,x_scale))+
    geom_vline(xintercept=0)+
    geom_hline(yintercept=-log10(0.05))+
    geom_text_repel(size=2,segment.size=0.1,max.overlaps=Inf,min.segment.length = 0, nudge_y = 1, segment.angle = 20,color="black") +
    theme(legend.position = "none",axis.title.x = element_blank(),axis.title.y = element_blank())
    ggsave(plt,file=paste0(outname,"_",assay_name,"_ESR1vFOXM1bins.pdf"),width=2,units="in",height=3)
  system(paste0("slack -F ",outname,"_",assay_name,"_ESR1vFOXM1bins.pdf"," ryan_todo"))
  write.table(x_da,file=paste0(outname,"_",assay_name,"_ESR1vFOXM1bins.markers.txt"),col.names=T,sep="\t")
  system(paste0("slack -F ",paste0(outname,"_",assay_name,"_ESR1vFOXM1bins.markers.txt")," ryan_todo"))
}

#Run DA/DE based on titan topic bins
find_markers_topic_bin(x=mcf7,outname="yw_mcf7.control",assay_name="chromvar")
find_markers_topic_bin(x=mcf7,outname="yw_mcf7.control",assay_name="cistrome_score",cistrome_cellline="MCF")
find_markers_topic_bin(x=mcf7,outname="yw_mcf7.control",assay_name="SCT")

find_markers_topic_bin(x=t47d,outname="yw_t47d.control",assay_name="chromvar")
find_markers_topic_bin(x=t47d,outname="yw_t47d.control",assay_name="cistrome_score",cistrome_cellline="T47")
find_markers_topic_bin(x=t47d,outname="yw_t47d.control",assay_name="SCT")

find_markers_topic_bin(x=mcf7,outname="yw_mcf7.control",assay_name="peaks",latent.vars="atac_peak_region_fragments")
find_markers_topic_bin(x=t47d,outname="yw_t47d.control",assay_name="peaks",latent.vars="atac_peak_region_fragments")


#Plot blend of foxm1 and esr1 topics over the umap
plt<-FeaturePlot(object = mcf7, features = c(mcf7_foxm1_topic, mcf7_esr1_topic),blend = T,col=c("white","red","blue"),reduction="multimodal_topic_umap",order=T)
ggsave(plt,file="mcf7_titantopic.umap.pdf",width=13)
system("slack -F mcf7_titantopic.umap.pdf ryan_todo")

plt<-FeaturePlot(object = t47d, features = c(t47d_foxm1_topic, t47d_esr1_topic),blend = T,col=c("white","red","blue"),reduction="multimodal_topic_umap",order=T)
ggsave(plt,file="t47d_titantopic.umap.pdf",width=13)
system("slack -F t47d_titantopic.umap.pdf ryan_todo")

#Correlate cistopic with titan topics
correlate_cistopic_titan<-function(x,outname,esr1_topic,foxm1_topic){
  cistopic_dat<-x@reductions$cistopic@cell.embeddings #get cistopic from x
  titan_dat<-as.data.frame(x@reductions$titan_lda@cell.embeddings) #get titan from y
  cells<-row.names(titan_dat)[which(row.names(titan_dat) %in% row.names(cistopic_dat))] # set list of shared cell names
  sum(cells %in% row.names(titan_dat))==sum(cells %in% row.names(titan_dat)) #check to make sure same number of cells
  titan_dat<-titan_dat[cells,]#reorder data frames
  cistopic_dat<-cistopic_dat[cells,]
  cor_mat<-data.frame(matrix(ncol = ncol(cistopic_dat), nrow = ncol(titan_dat))) #make an empty data frame to populate
  colnames(cor_mat)<-colnames(cistopic_dat)
  row.names(cor_mat)<-colnames(titan_dat)

  for (yi in 1:ncol(cistopic_dat)){
    for (xj in 1:ncol(titan_dat)){
      cor_mat[xj,yi]<-cor(x=as.numeric(cistopic_dat[,yi]),y=as.numeric(titan_dat[,xj]))
    }
  }

  #perform a weighted linear regression, TITAN topic is weight and cistopic topics are comparisons
  list_plots<-c()
  for (xi in colnames(cistopic_dat)){ 
    for (yi in colnames(cistopic_dat)){ 
      for (zj in c(foxm1_topic)){

      #generate a linear model
      x=as.numeric(cistopic_dat[,xi])
      y=as.numeric(cistopic_dat[,yi])
      x2_titan=as.numeric(titan_dat[,zj])
      model_fit<-summary(lm(y~x+x2_titan))
      intercept<-model_fit$coefficients[1,1]
      slope<-model_fit$coefficients[2,1]
      std_err<-model_fit$coefficients[2,2]
      residuals<-model_fit$residuals

    plt_dat<-data.frame(x=x,y=y,x2_titan=x2_titan)
    plt<-ggplot(plt_dat,aes(x=x,y=y,color=x2_titan))+
    geom_abline(intercept = intercept, slope = slope, col = "red")+
    geom_point(size=0.5,alpha=0.5)+
    scale_color_gradient2()+
    theme_void()+
    xlab(xi)+
    ylab(yi)
    list_plots[[paste(xi,yi)]]<-plt
    }
  }
}

  plt_list<-wrap_plots(list_plots,nrow=length(colnames(cistopic_dat)),guides="collect")+theme(legend.position = "none")
  ggsave(plt_list,file=paste0(outname,"_cistopic_scatterplot.png"),height=50,width=50,limitsize=FALSE)
  system(paste0("slack -F ",outname,"_cistopic_scatterplot.png ryan_todo"))

  out_plt<-Heatmap(cor_mat)
  columnorder<-column_order(out_plt)
  pdf(paste0(outname,".cistopic_titan_topiccorrelation.pdf"))
  print(out_plt)
  dev.off()
  system(paste0("slack -F ",paste0(outname,".cistopic_titan_topiccorrelation.pdf")," ryan_todo"))

  #add cistopic square correlation matrix
  cor_cistopic<-cor(cistopic_dat)
  out_plt<-Heatmap(cor_cistopic)
  pdf(paste0(outname,".cistopic_squarecorrelation.pdf"))
  print(out_plt)
  dev.off()
  system(paste0("slack -F ",paste0(outname,".cistopic_squarecorrelation.pdf")," ryan_todo"))

  #titan square correlation
  cor_titan<-cor(titan_dat)
  out_plt<-Heatmap(cor_titan)
  pdf(paste0(outname,".titan_squarecorrelation.pdf"))
  print(out_plt)
  dev.off()
  system(paste0("slack -F ",paste0(outname,".titan_squarecorrelation.pdf")," ryan_todo"))

  cor_mat_topics<-cor_mat[c(esr1_topic,foxm1_topic),]
  #cor_mat_topics<-rbind(cor_mat_topics,abs(cor_mat_topics[esr1_topic,]-cor_mat_topics[foxm1_topic,]))
  #row.names(cor_mat_topics)<-c(esr1_topic,foxm1_topic,"difference")
  cor_mat_topics<-rbind(abs(cor_mat_topics[esr1_topic,]-cor_mat_topics[foxm1_topic,]))
  row.names(cor_mat_topics)<-c("difference")
  col_fun = colorRamp2(c(0, 0.5), c("white", "red"))

  out_plt<-Heatmap(cor_mat_topics,
    col=col_fun,
    column_order=columnorder)
  pdf(paste0(outname,".cistopic_titan_topiccorrelation_specific.pdf"))
  print(out_plt)
  dev.off()
  system(paste0("slack -F ",paste0(outname,".cistopic_titan_topiccorrelation_specific.pdf")," ryan_todo"))

}

#x is seurat object with redcutions for cistopic and titan_lda
correlate_cistopic_titan(x=mcf7,outname="yw_mcf7.control",foxm1_topic=mcf7_foxm1_topic,esr1_topic=mcf7_esr1_topic)
correlate_cistopic_titan(x=t47d,outname="yw_t47d.control",foxm1_topic=t47d_foxm1_topic,esr1_topic=t47d_esr1_topic)

#now to correlate chromvar motifs with titan topics

#Correlate chromvar with titan topics
correlate_chromvar_titan<-function(x,outname,cellline="test",quantile_cutoff=0.95,esr1_topic,foxm1_topic){
  chromvar_dat<-as.data.frame(t(x@assays$chromvar@data)) #get chromvar from x
  tfList <- getMatrixByID(JASPAR2020, ID=colnames(chromvar_dat)) #set up readable chromvar names
  tfList <-unlist(lapply(names(tfList), function(x) name(tfList[[x]])))
  colnames(chromvar_dat)<-tfList
  titan_dat<-as.data.frame(x@reductions$titan_lda@cell.embeddings) #get titan from y  
  cells<-row.names(titan_dat)[which(row.names(titan_dat) %in% row.names(chromvar_dat))] # set list of shared cell names
  sum(cells %in% row.names(titan_dat))==sum(cells %in% row.names(titan_dat)) #check to make sure same number of cells
  titan_dat<-titan_dat[cells,]#reorder data frames
  chromvar_dat<-chromvar_dat[cells,]
  cor_mat<-data.frame(matrix(ncol = ncol(chromvar_dat), nrow = ncol(titan_dat))) #make an empty data frame to populate
  colnames(cor_mat)<-colnames(chromvar_dat)
  row.names(cor_mat)<-colnames(titan_dat)

  for (yi in 1:ncol(chromvar_dat)){
    for (xj in 1:ncol(titan_dat)){
      cor_mat[xj,yi]<-cor(x=chromvar_dat[,yi],y=titan_dat[,xj],use="complete.obs")
    }
  }


  out_plt<-Heatmap(t(cor_mat),
    row_names_gp = grid::gpar(fontsize = 4),
    row_km=5
    )
  pdf(paste0(outname,".chromvar_titan_topiccorrelation.pdf"),height=30)
  print(out_plt)
  dev.off()
  system(paste0("slack -F ",paste0(outname,".chromvar_titan_topiccorrelation.pdf")," ryan_todo"))

  cor_mat<-as.data.frame(t(cor_mat[c(esr1_topic,foxm1_topic),]))
 
  #generate a linear model
  x=as.numeric(cor_mat[,esr1_topic])
  y=as.numeric(cor_mat[,foxm1_topic])
  model_fit<-summary(lm(y~x,data=cor_mat))
  intercept<-model_fit$coefficients[1,1]
  slope<-model_fit$coefficients[2,1]
  std_err<-model_fit$coefficients[2,2]
  residuals<-model_fit$residuals

  #highlight points which were determined significant
  assay_name="chromvar"
  x_da<-read.table(paste0(outname,"_",assay_name,"_ESR1vFOXM1bins.markers.txt"))
  cor_mat$color_fill<-ifelse(abs(residuals)>std_err,"red","not_red")

  ident_1_labels<- row.names(head(x_da[which((x_da$sig=="sig") & (x_da$avg_log2FC>0)),] ,n=20))#significant, is ident1 specific, is top 20
  ident_2_labels<- row.names(head(x_da[which((x_da$sig=="sig") & (x_da$avg_log2FC<0)),] ,n=20))#significant, is ident1 specific, is top 20
  x_da[c(ident_1_labels,ident_2_labels),]$label<-x_da[c(ident_1_labels,ident_2_labels),]$out_name

  cor_mat$label<-unlist(lapply(1:nrow(cor_mat),function(y) {ifelse(row.names(cor_mat)[y] %in% x_da$label,row.names(cor_mat)[y],"")})) #add label to significant features

  plt<-ggplot()+
  geom_point(aes(x=cor_mat[,esr1_topic],y=cor_mat[,foxm1_topic],color=cor_mat$color_fill),size=0.1)+
  scale_color_manual(values=c("red"="red","not_red"="grey50"))+
  theme_minimal()+
  geom_text_repel(aes(x=cor_mat[,esr1_topic],y=cor_mat[,foxm1_topic],color=cor_mat$color_fill,label=cor_mat$label),force=5,max.overlaps=Inf,size=2,min.segment.length = 0,box.padding=0.3,segment.size=0.1)+
  geom_abline(intercept = intercept, slope = slope, col = "black")+
  ggtitle(paste(outname,"chromvar v. titan correlations"))+
  xlab(paste("ESR1",esr1_topic))+
  ylab(paste("FOXM1",foxm1_topic))+
  theme(axis.text = element_text(size = 6),axis.title=element_text(size=8),axis.ticks=element_blank(),legend.position="none")+
  scale_x_continuous(expand = expansion(mult = 0.5))+  scale_y_continuous(expand = expansion(mult = 0.5))

  ggsave(plt,file=paste0(outname,".chromvar_titan_topiccorrelation.comparison.pdf"),height=4,width=4,units="in")
  print(plt)
  dev.off()
  system(paste0("slack -F ",paste0(outname,".chromvar_titan_topiccorrelation.comparison.pdf")," ryan_todo"))
}

correlate_chromvar_titan(x=mcf7,outname="yw_mcf7.control",foxm1_topic=mcf7_foxm1_topic,esr1_topic=mcf7_esr1_topic)
correlate_chromvar_titan(x=t47d,outname="yw_t47d.control",foxm1_topic=t47d_foxm1_topic,esr1_topic=t47d_esr1_topic)

#Correlate cistrome data with titan topics
correlate_cistrome_titan<-function(x,outname,cellline="MCF",quantile_cutoff=0.95,esr1_topic,foxm1_topic){
  cistrome_dat<-as.data.frame(t(x@assays$cistrome_score@data)) #get cistrome from x
  titan_dat<-as.data.frame(x@reductions$titan_lda@cell.embeddings) #get titan from y  
  cells<-row.names(titan_dat)[which(row.names(titan_dat) %in% row.names(cistrome_dat))] # set list of shared cell names
  sum(cells %in% row.names(titan_dat))==sum(cells %in% row.names(titan_dat)) #check to make sure same number of cells
  titan_dat<-titan_dat[cells,]#reorder data frames
  cistrome_dat<-cistrome_dat[cells,]
  cor_mat<-data.frame(matrix(ncol = ncol(cistrome_dat), nrow = ncol(titan_dat))) #make an empty data frame to populate
  colnames(cor_mat)<-colnames(cistrome_dat)
  row.names(cor_mat)<-colnames(titan_dat)
  for (yi in 1:ncol(cistrome_dat)){
    for (xj in 1:ncol(titan_dat)){
      cor_mat[xj,yi]<-cor(x=cistrome_dat[,yi],y=titan_dat[,xj],use="complete.obs")
    }
  }

  cor_mat<-cor_mat[,startsWith(colnames(cor_mat),cellline)] #subset to same cell line cistrome data
  out_plt<-Heatmap(t(cor_mat),
    row_names_gp = grid::gpar(fontsize = 4),
    row_km=5
    )
  pdf(paste0(outname,".cistrome_titan_topiccorrelation.pdf"),height=30)
  print(out_plt)
  dev.off()
  system(paste0("slack -F ",paste0(outname,".cistrome_titan_topiccorrelation.pdf")," ryan_todo"))


  cor_mat<-as.data.frame(t(cor_mat[c(esr1_topic,foxm1_topic),]))
  
  #generate a linear model
  x=as.numeric(cor_mat[,esr1_topic])
  y=as.numeric(cor_mat[,foxm1_topic])
  model_fit<-summary(lm(y~x,data=cor_mat))
  intercept<-model_fit$coefficients[1,1]
  slope<-model_fit$coefficients[2,1]
  std_err<-model_fit$coefficients[2,2]
  residuals<-model_fit$residuals

  #highlight points which were determined significant
  assay_name="cistrome_score"
  x_da<-read.table(paste0(outname,"_",assay_name,"_ESR1vFOXM1bins.markers.txt"))
  ident_1_labels<- row.names(head(x_da[which((x_da$sig=="sig") & (x_da$avg_log2FC>0)),] ,n=10))#significant, is ident1 specific, is top 20
  ident_2_labels<- row.names(head(x_da[which((x_da$sig=="sig") & (x_da$avg_log2FC<0)),] ,n=10))#significant, is ident1 specific, is top 20
  sig_labels<-row.names(x_da[which(x_da$sig=="sig"),])#significant, is ident1 specific, is top 20
  x_da[c(ident_1_labels,ident_2_labels),]$label<-x_da[c(ident_1_labels,ident_2_labels),]$out_name
  cor_mat$label<-unlist(lapply(1:nrow(cor_mat),function(y) {ifelse(row.names(cor_mat)[y] %in% x_da$label,row.names(cor_mat)[y],"")})) #add label to significant features
  cor_mat$color_fill<-unlist(lapply(1:nrow(cor_mat),function(y) {ifelse(row.names(cor_mat)[y] %in% sig_labels,"red","not_red")})) #add label to significant features

  plt<-ggplot()+
  geom_point(aes(x=cor_mat[,esr1_topic],y=cor_mat[,foxm1_topic],color=cor_mat$color_fill),size=0.1)+
  scale_color_manual(values=c("red"="red","not_red"="grey50"))+
  theme_minimal()+
  geom_text_repel(aes(x=cor_mat[,esr1_topic],y=cor_mat[,foxm1_topic],color=cor_mat$color_fill,label=cor_mat$label),force=5,max.overlaps=Inf,size=2,min.segment.length = 0,box.padding=0.3,segment.size=0.1)+
  geom_abline(intercept = intercept, slope = slope, col = "black")+
  ggtitle(paste(outname,"cistrome v. titan correlations"))+
  xlab(paste("ESR1",esr1_topic))+
  ylab(paste("FOXM1",foxm1_topic))+
  theme(axis.text = element_text(size = 6),axis.title=element_text(size=8),axis.ticks=element_blank(),legend.position="none")+
  scale_x_continuous(expand = expansion(mult = 0.5))+  scale_y_continuous(expand = expansion(mult = 0.5))

  ggsave(plt,file=paste0(outname,".cistrome_titan_topiccorrelation.comparison.pdf"),height=4,width=4,units="in")
  print(plt)
  dev.off()
  system(paste0("slack -F ",paste0(outname,".cistrome_titan_topiccorrelation.comparison.pdf")," ryan_todo"))
}

correlate_cistrome_titan(x=mcf7,outname="yw_mcf7.control",foxm1_topic=mcf7_foxm1_topic,esr1_topic=mcf7_esr1_topic,cellline="MCF")
correlate_cistrome_titan(x=t47d,outname="yw_t47d.control",foxm1_topic=t47d_foxm1_topic,esr1_topic=t47d_esr1_topic,cellline="T47")



plt1<-FeaturePlot(mcf7,features=c("topic_20"),reduction="multimodal_topic_umap",cols=c("white","blue"))
plt2<-FeaturePlot(mcf7,features=c("topic_13"),reduction="multimodal_topic_umap",cols=c("white","red"))
ggsave(plt1|plt2,file="yw_mcf7.control.cistopic_differences.pdf",width=10)
system("slack -F yw_mcf7.control.cistopic_differences.pdf ryan_todo")


plt1<-FeaturePlot(t47d,features=c("topic_3"),reduction="multimodal_topic_umap",cols=c("white","blue"))
plt2<-FeaturePlot(t47d,features=c("topic_2"),reduction="multimodal_topic_umap",cols=c("white","red"))
ggsave(plt1|plt2,file="yw_t47d.control.cistopic_differences.pdf",width=10)
system("slack -F yw_t47d.control.cistopic_differences.pdf ryan_todo")