CCAN_analyses.Rmd

---
title: "R CCAN analysee"
output: html_notebook
---




Internal functions
```{r}
corrected_reads_CCAN<-function(object,assay,ccan,group_name,species=NULL){

list_of_matrices=list()
  
if(is.null(species)) {
species=c("mouse","human")
}

for (sp in species){
#do it per species
print(sp)
object_sub<-subset(object,subset = species == paste0(sp))
object_sub<-BinarizeCounts(object_sub)

average_count_per_group<-AverageExpression(object = object_sub,assay = assay,slot = "counts",group.by = group_name)
groups<-as.factor(object_sub[[]][,group_name])
#at this point we have the probability of accessibility for each cluster. To adjust for complexity we have to calculate the median number of sites accessible in each cluster we take the average of this and divide them by individual groups ()

median_nfeature_category<-aggregate(as.data.frame(object_sub$nFeature_cre),by=list(groups),FUN=median)

median_nfeature_category$norm_complexity<-mean(log10(median_nfeature_category$`object_sub$nFeature_cre`))/log10(median_nfeature_category$`object_sub$nFeature_cre`)
norm_factor<-as.data.frame(median_nfeature_category$norm_complexity)
row.names(norm_factor)<-median_nfeature_category$Group.1

return_matrix<-average_count_per_group$cre*t(replicate(dim(average_count_per_group$cre)[1],norm_factor$`median_nfeature_category$norm_complexity`))


ccan$CCAN<-paste0("CCAN_",ccan$CCAN)
ccan$Peak<-NULL
return_matrix<-merge(return_matrix,ccan,by="row.names")
row.names(return_matrix)<-return_matrix$Row.names
return_matrix$Row.names<-NULL

CCAN_groups<-as.factor(return_matrix$CCAN)
return_matrix$CCAN<-NULL
return_matrix<-aggregate(return_matrix,by=list(CCAN_groups),FUN = mean)


list_of_matrices[[paste0(sp)]]=return_matrix
}


return(list_of_matrices)


}

  # region can be a string, name of a gene, or GRanges object
FindRegion <- function(
  object,
  region,
  sep = c("-", "-"),
  assay = NULL,
  extend.upstream = 0,
  extend.downstream = 0
) {
  if (!is(object = region, class2 = "GRanges")) {
    # if separators are present in the string and we can convert the
    # start to a number, assume we're using genomic coordinates
    if (all(sapply(X = sep, FUN = grepl, x = region))) {
      region <- StringToGRanges(regions = region, sep = sep)
    } else {
      region <- LookupGeneCoords(object = object, assay = assay, gene = region)
      if (is.null(x = region)) {
        stop("Gene not found")
      }
    }
  }
  region <- suppressWarnings(expr = Extend(
    x = region,
    upstream = extend.upstream,
    downstream = extend.downstream
  )
  )
  return(region)
}


```
CCAN analyses
```{r}
library(Signac)
library(Seurat)
library(SeuratWrappers)
library(ggplot2)
library(patchwork)
library(GenomicRanges)

human<-readRDS("/brahms/torkenczyk/H_M_storage/analysis/objects/Select_human_tissues.rds")
mouse<-readRDS("/brahms/torkenczyk/H_M_storage/analysis/objects/Select_mouse_tissues.rds")



human
mouse
library(cicero)


#hg38 human all
human.cds <- as.cell_data_set(x = human)
human.cicero <- make_cicero_cds(human.cds, reduced_coordinates = reducedDims(human.cds)$UMAP)
genome <- seqlengths(human)
# convert chromosome sizes to a dataframe
genome.df <- data.frame("chr" = names(genome), "length" = genome)

# run cicero
conns_human <- run_cicero(human.cicero, genomic_coords = genome.df, sample_num = 100)
write.table(conns_human,file="Human_connections_cicero.txt",sep = "\t",quote = F,row.names = T,col.names = T)

ccans_human <- generate_ccans(conns_human,coaccess_cutoff_override = 0.15)
ccans_human <- generate_ccans(conns_human)



write.table(ccans_human,file="Human_ccans_cicero.txt",sep = "\t",quote = F,row.names = T,col.names = T)
```



```{r}
DefaultAssay(mouse)<-"mm9"
#do same for moust but in mm9
#mouse[["cre"]]<-NULL
mouse[["cre05"]]<-NULL
mouse[["cre99"]]<-NULL
mouse[["Deviation_DA_mouse"]]<-NULL
mouse[["Deviation_DA_human"]]<-NULL
mouse[["Deviation_DA_int"]]<-NULL
mouse[["chromvar"]]<-NULL
mouse[["mm10"]]<-NULL

#mm9 mouse but did it with hg38 as well
mouse.cds <- as.cell_data_set(x = mouse)
mouse.cicero <- make_cicero_cds(mouse.cds, reduced_coordinates = reducedDims(mouse.cds)$UMAP)
genome <- seqlengths(mouse)
# convert chromosome sizes to a dataframe
genome.df <- data.frame("chr" = names(genome), "length" = genome)

# run cicero
conns_mouse <- run_cicero(mouse.cicero, genomic_coords = genome.df, sample_num = 100)
write.table(conns_mouse,file="Mouse_connections_cicero.txt",sep = "\t",quote = F,row.names = T,col.names = T)

#ccans_mouse <- generate_ccans(conns_mouse,coaccess_cutoff_override = 0.15)
ccans_mouse <- generate_ccans(conns_mouse)


write.table(ccans_mouse,file="Mouse_ccans_cicero.txt",sep = "\t",quote = F,row.names = T,col.names = T)


#now mm10
#mm10 mouse but did it with hg38 as well
DefaultAssay(mouse)<-"mm10"
mouse[["mm9"]]<-NULL
mouse[["cre"]]<-NULL

mouse.cds <- as.cell_data_set(x = mouse)
mouse.cicero <- make_cicero_cds(mouse.cds, reduced_coordinates = reducedDims(mouse.cds)$UMAP)
genome <- seqlengths(mouse)
# convert chromosome sizes to a dataframe
genome.df <- data.frame("chr" = names(genome), "length" = genome)
# run cicero
conns_mouse_mm10 <- run_cicero(mouse.cicero, genomic_coords = genome.df, sample_num = 100)
write.table(conns_mouse_mm10,file="Mouse_connections_cicero_mm10.txt",sep = "\t",quote = F,row.names = T,col.names = T)
#ccans_mouse <- generate_ccans(conns_mouse,coaccess_cutoff_override = 0.15)
ccans_mouse_mm10 <- generate_ccans(conns_mouse_mm10)
write.table(ccans_mouse_mm10,file="Mouse_ccans_cicero_mm10.txt",sep = "\t",quote = F,row.names = T,col.names = T)

mouse<-readRDS("/brahms/torkenczyk/H_M_storage/analysis//active_use_objects/Select_mouse_tissues.rds")

DefaultAssay(mouse)<-"cre"
mouse[["mm9"]]<-NULL
mouse[["mm10"]]<-NULL

mouse.cds <- as.cell_data_set(x = mouse)
mouse.cicero <- make_cicero_cds(mouse.cds, reduced_coordinates = reducedDims(mouse.cds)$UMAP)
genome <- seqlengths(mouse)
# convert chromosome sizes to a dataframe
genome.df <- data.frame("chr" = names(genome), "length" = genome)
# run cicero
conns_mouse_hg38 <- run_cicero(mouse.cicero, genomic_coords = genome.df, sample_num = 100)
write.table(conns_mouse_hg38,file="Mouse_connections_cicero_mm10.txt",sep = "\t",quote = F,row.names = T,col.names = T)
#ccans_mouse <- generate_ccans(conns_mouse,coaccess_cutoff_override = 0.15)
ccans_mouse_hg38 <- generate_ccans(conns_mouse_hg38)
write.table(ccans_mouse_hg38,file="Mouse_ccans_cicero_hg38.txt",sep = "\t",quote = F,row.names = T,col.names = T)

mouse<-readRDS("/brahms/torkenczyk/H_M_storage/analysis//active_use_objects/Select_mouse_tissues.rds")
DefaultAssay(mouse)<-"mm9"
links <- ConnectionsToLinks(conns = conns_mouse, ccans = ccans_mouse)
Links(mouse) <- links
DefaultAssay(mouse)<-"mm10"
links <- ConnectionsToLinks(conns = conns_mouse_mm10, ccans = ccans_mouse_mm10)
Links(mouse) <- links
DefaultAssay(mouse)<-"cre"
links <- ConnectionsToLinks(conns = conns_mouse_hg38, ccans = ccans_mouse_hg38)
Links(mouse) <- links
saveRDS(mouse,file = "/brahms/torkenczyk/H_M_storage/analysis//active_use_objects/Mouse_added_links.rds")




load("/brahms/torkenczyk/H_M_storage/analysis/objects/DA_info/Shared.bw.species.DA.RData")
load("/brahms/torkenczyk/H_M_storage/analysis/objects/DA_info/Human_only.bw.species.DA.RData")
load("/brahms/torkenczyk/H_M_storage/analysis/objects/DA_info/Mouse_only.bw.species.DA.RData")


all_human_peaks_in_CCAN<-unique(StringToGRanges(ccans_human$Peak))
all_mouse_peaks_in_CCAN<-unique(StringToGRanges(ccans_mouse$Peak))
all_peaks_human<-unique(StringToGRanges(row.names(human)))
all_peaks_mouse<-unique(StringToGRanges(row.names(mouse)))


#add links to human

links <- ConnectionsToLinks(conns = conns_human, ccans = ccans_human)
Links(human) <- links
saveRDS(human,file = "/brahms/torkenczyk/H_M_storage/analysis/objects/Human_added_links.rds")



#hypergeom
#easy hypergeom
hypergeom_test=data.frame(celltype=character(),CCAN=character(),group=character(),hgtest=double())
percent_overlap=data.frame(celltype=character(),CCAN=character(),group=character(),frac=double())
for (i in names(shared_DA))
{
print(i)
  print("shared")
#shared DA sites
querypos<-unique(StringToGRanges(sort(shared_DA[[i]]$feature_hg38_universal)))

#human CCAN
mpre<-suppressWarnings(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(all_peaks_human)))
xpre<-suppressWarnings(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(querypos)))
m<-length(unique(queryHits(mpre)))
x<-length(unique(queryHits(xpre)))

print(paste0("Number of peaks in CCAN: ",m))
print(paste0("Number of peaks not in CCAN: ",length(all_peaks_human)-m))
print(paste0("Frac of peaks in CCAN/all peaks: ",m/length(all_peaks_human)))
print(paste0("enrichment in shared: ",x))
print(paste0("enrichment in shared/thousand DA peaks: ",x/1000))

shared_sites_in_CCANS<-sort(querypos)[(queryHits(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(querypos))))]
CCAN_sites_in_sh_sites<-all_human_peaks_in_CCAN[(subjectHits(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(querypos))))]
sh<-data.frame(group="shared",celltype=i,CCAN="human",hgtest=(1-phyper(x-1,m,length(all_peaks_human)-m,length(querypos))))
hypergeom_test<-rbind(hypergeom_test,sh)
sh2<-data.frame(group="shared",celltype=i,CCAN="human",frac=x/1000*length(all_peaks_human)/m)
percent_overlap<-rbind(percent_overlap,sh2)



#mouse CCAN
#shared DA sites
querypos<-unique(StringToGRanges(sort(shared_DA[[i]]$feature_mm9_mouse)))
mpre<-suppressWarnings(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(all_peaks_mouse)))
xpre<-suppressWarnings(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(querypos)))
m<-length(unique(queryHits(mpre)))
x<-length(unique(queryHits(xpre)))

print(paste0("Number of peaks in CCAN: ",m))
print(paste0("Number of peaks not in CCAN: ",length(all_peaks_mouse)-m))
print(paste0("Frac of peaks in CCAN/all peaks: ",m/length(all_peaks_mouse)))
print(paste0("enrichment in shared: ",x))
print(paste0("enrichment in shared/thousand DA peaks: ",x/1000))



shared_sites_in_CCANS<-sort(querypos)[(queryHits(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(querypos))))]
CCAN_sites_in_sh_sites<-all_mouse_peaks_in_CCAN[(subjectHits(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(querypos))))]
sh<-data.frame(group="shared",celltype=i,CCAN="mouse",hgtest=(1-phyper(x-1,m,length(all_peaks_mouse)-m,length(querypos))))
hypergeom_test<-rbind(hypergeom_test,sh)
sh2<-data.frame(group="shared",celltype=i,CCAN="mouse",frac=x/1000*length(all_peaks_human)/m)
percent_overlap<-rbind(percent_overlap,sh2)


#mouse DA sites
querypos<-unique(StringToGRanges(sort(mouse_DA[[i]]$feature_hg38_universal)))

#human CCAN
mpre<-suppressWarnings(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(all_peaks_human)))
xpre<-suppressWarnings(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(querypos)))
m<-length(unique(queryHits(mpre)))
x<-length(unique(queryHits(xpre)))
print("mouse DA human CCANS")
print(paste0("enrichment in mouse DA: ",x))
print(paste0("enrichment in mouse DA/thousand DA peaks: ",x/1000))
print(paste0("Frac of peaks in CCAN/all peaks: ",m/length(all_peaks_human)))

mouse_sites_in_CCANS<-sort(querypos)[(queryHits(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(querypos))))]
CCAN_sites_in_mouse_sites<-all_human_peaks_in_CCAN[(subjectHits(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(querypos))))]
sh<-data.frame(group="mouse",celltype=i,CCAN="human",hgtest=(1-phyper(x-1,m,length(all_peaks_human)-m,length(querypos))))
hypergeom_test<-rbind(hypergeom_test,sh)
sh2<-data.frame(group="mouse",celltype=i,CCAN="human",frac=x/1000*length(all_peaks_human)/m)
percent_overlap<-rbind(percent_overlap,sh2)



#mouse CCAN
#mouse DA sites
querypos<-unique(StringToGRanges(sort(mouse_DA[[i]]$feature_mm9_mouse)))
mpre<-suppressWarnings(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(all_peaks_mouse)))
xpre<-suppressWarnings(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(querypos)))
m<-length(unique(queryHits(mpre)))
x<-length(unique(queryHits(xpre)))

print("mouse DA mouse CCANS")
print(paste0("enrichment in mouse DA: ",x))
print(paste0("enrichment in mouse DA/thousand DA peaks: ",x/1000))
print(paste0("Frac of peaks in CCAN/all peaks: ",m/length(all_peaks_mouse)))

mouse_sites_in_CCANS<-sort(querypos)[(queryHits(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(querypos))))]
CCAN_sites_in_mouse_sites<-all_mouse_peaks_in_CCAN[(subjectHits(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(querypos))))]
sh<-data.frame(group="mouse",celltype=i,CCAN="mouse",hgtest=(1-phyper(x-1,m,length(all_peaks_mouse)-m,length(querypos))))
hypergeom_test<-rbind(hypergeom_test,sh)
sh2<-data.frame(group="mouse",celltype=i,CCAN="mouse",frac=x/1000*length(all_peaks_human)/m)
percent_overlap<-rbind(percent_overlap,sh2)


#human DA sites
querypos<-unique(StringToGRanges(sort(human_DA[[i]]$feature_hg38_universal)))

#human CCAN
mpre<-suppressWarnings(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(all_peaks_human)))
xpre<-suppressWarnings(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(querypos)))
m<-length(unique(queryHits(mpre)))
x<-length(unique(queryHits(xpre)))

print("human DA human CCANS")
print(paste0("enrichment in human DA: ",x))
print(paste0("enrichment in human DA/thousand DA peaks: ",x/1000))
print(paste0("Frac of peaks in CCAN/all peaks: ",m/length(all_peaks_human)))

human_sites_in_CCANS<-sort(querypos)[(queryHits(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(querypos))))]
CCAN_sites_in_human_sites<-all_human_peaks_in_CCAN[(subjectHits(findOverlaps(subject = all_human_peaks_in_CCAN,query = sort(querypos))))]
sh<-data.frame(group="human",celltype=i,CCAN="human",hgtest=(1-phyper(x-1,m,length(all_peaks_human)-m,length(querypos))))
hypergeom_test<-rbind(hypergeom_test,sh)
sh2<-data.frame(group="human",celltype=i,CCAN="human",frac=x/1000*length(all_peaks_human)/m)
percent_overlap<-rbind(percent_overlap,sh2)



#mouse CCAN
#human DA sites
querypos<-unique(StringToGRanges(sort(human_DA[[i]]$feature_mm9_mouse)))
mpre<-suppressWarnings(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(all_peaks_mouse)))
xpre<-suppressWarnings(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(querypos)))
m<-length(unique(queryHits(mpre)))
x<-length(unique(queryHits(xpre)))

print("human DA mouse CCANS")
print(paste0("enrichment in human DA: ",x))
print(paste0("enrichment in human DA/thousand DA peaks: ",x/1000))
print(paste0("Frac of peaks in CCAN/all peaks: ",m/length(all_peaks_mouse)))

human_sites_in_CCANS<-sort(querypos)[(queryHits(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(querypos))))]
CCAN_sites_in_human_sites<-all_mouse_peaks_in_CCAN[(subjectHits(findOverlaps(subject = all_mouse_peaks_in_CCAN,query = sort(querypos))))]
sh<-data.frame(group="human",celltype=i,CCAN="mouse",hgtest=(1-phyper(x-1,m,length(all_peaks_mouse)-m,length(querypos))))
hypergeom_test<-rbind(hypergeom_test,sh)
sh2<-data.frame(group="human",celltype=i,CCAN="mouse",frac=x/1000*length(all_peaks_human)/m)
percent_overlap<-rbind(percent_overlap,sh2)

}

percent_overlap_m<-percent_overlap[which(percent_overlap$CCAN=="mouse"),]
percent_overlap_h<-percent_overlap[which(percent_overlap$CCAN=="human"),]

M<-t(scale(t(acast(percent_overlap_m, celltype~group, value.var="frac"))))
H<-t(scale(t(acast(percent_overlap_h, celltype~group, value.var="frac"))))

HM<-Heatmap(H)+Heatmap(M)

write.table(percent_overlap,file="Percent_overlap_CCANs.txt",quote = F,sep = "\t",row.names = T,col.names = T)


p1<-ggplot(percent_overlap,aes(x=CCAN,fill=group,y=frac))+geom_boxplot()+theme_classic()+ylab("Fraction of DA peaks in a CCAN/Fraction of all peaks in a CCAN")

library(ggpubr)

percent_overlap_h<-percent_overlap[percent_overlap$CCAN=="human", ]
percent_overlap_m<-percent_overlap[percent_overlap$CCAN=="mouse",]

percent_overlap_h2<-data.frame(shared=percent_overlap_h[percent_overlap_h$group=="shared",]$frac,mouse=percent_overlap_h[percent_overlap_h$group=="mouse",]$frac,human=percent_overlap_h[percent_overlap_h$group=="human",]$frac)


# Visualize: Specify the comparisons you want
my_comparisons <- list( c("shared","human"), c("mouse", "shared"))
p_h_comp<-ggboxplot(percent_overlap_h, x = "group", y = "frac",fill = "group", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 5)+ylab("Fraction of DA peaks in a CCAN/Fraction of all peaks in a CCAN")


#should maybe test paired. same result

p_h_comp<-ggpaired(percent_overlap_h2, cond1 = "shared", cond2 = "human",line.color = "gray", line.size = 0.4)+ 
  stat_compare_means(paired = T)+ylab("Fraction of DA peaks in a CCAN/Fraction of all peaks in a CCAN")



p_m_comp<-ggboxplot(percent_overlap_m, x = "group", y = "frac",fill = "group",palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons,paired = F)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 5)+ylab("Fraction of DA peaks in a CCAN/Fraction of all peaks in a CCAN")



# Visualize: Specify the comparisons you want
my_comparisons <- list( c("human", "shared"), c("mouse", "shared"))
p_all<-ggboxplot(percent_overlap, x = "group", y = "frac",fill="group", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 11)+ylab("Fraction of DA peaks in a CCAN/Fraction of all peaks in a CCAN")

ggsave(p1,filename = "~/H_M_chromatin_atlas/human_data/plots/paper_pdfs/Comparison_CCAN_shared_unique_CCAN_breakdown.pdf")
ggsave(p2,filename = "~/H_M_chromatin_atlas/human_data/plots/paper_pdfs/Comparison_CCAN_shared_unique.pdf")

```


```{r}
#Now lets look at examples with links that are unique

human<-readRDS("/brahms/torkenczyk/H_M_storage/analysis/objects/Human_added_links.rds")
mouse<-readRDS("/brahms/torkenczyk/H_M_storage/analysis/objects/Mouse_added_links.rds")
 
human[["Deviation_DA_mouse"]]<-NULL
human[["Deviation_DA_human"]]<-NULL
human[["Deviation_DA_int"]]<-NULL
human[["chromvar"]]<-NULL
mouse[["Deviation_DA_mouse"]]<-NULL
mouse[["Deviation_DA_human"]]<-NULL
mouse[["Deviation_DA_int"]]<-NULL
mouse[["mm9"]]<-NULL
mouse[["mm10"]]<-NULL
mouse[["cre05"]]<-NULL
mouse[["cre99"]]<-NULL
mouse[["chromvar"]]<-NULL
human$species="human"
mouse$species="mouse"

H_M_merged_cre<-merge(human,mouse)


ccans_human<-read.table("Human_ccans_cicero.txt")

H_M_merged_cre$celltype<-Idents(H_M_merged_cre)

human.cicero.ccan.aggregate<-corrected_reads_CCAN(object = H_M_merged_cre,assay = "cre",group_name = "celltype",ccan=ccans_human)

hum<-human.cicero.ccan.aggregate$human
row.names(hum)<-hum$Group.1
hum$Group.1<-NULL

mus<-human.cicero.ccan.aggregate$mouse
row.names(mus)<-mus$Group.1
mus$Group.1<-NULL

H<-Heatmap(hum,show_column_names = TRUE,show_row_names = FALSE,cluster_rows = TRUE)
H2<-Heatmap(t(scale(t(hum))),show_column_names = TRUE,show_row_names = FALSE,cluster_rows = TRUE,cluster_columns = F)

M<-Heatmap(mus,show_column_names = TRUE,show_row_names = FALSE,cluster_rows = TRUE)
M2<-Heatmap(t(scale(t(mus))),show_column_names = TRUE,show_row_names = FALSE,cluster_rows = TRUE,cluster_columns = F)

s_h_ccan<-as.data.frame(t(scale(t(hum))))
s_m_ccan<-as.data.frame(t(scale(t(mus))))

write.table(s_h_ccan,file="/brahms/torkenczyk/H_M_storage/analysis//Human_ccan.table.txt")
write.table(s_m_ccan,file="/brahms/torkenczyk/H_M_storage/analysis//Mouse_ccan.table.txt")



correlates_CCAN_M_H<-cor(x=s_h_ccan,y=s_m_ccan,use = "complete.obs")

Heatmap(correlates_CCAN_M_H,show_column_names = TRUE,show_row_names = TRUE,cluster_rows = TRUE,cluster_columns = T)



#ordering on the difference on accessibility
CARD_order<-s_h_ccan$Cardiomyocytes-s_m_ccan$Cardiomyocytes
CCANs_human_unique<-row.names(s_h_ccan[order(CARD_order,decreasing = T),]-s_m_ccan[order(CARD_order,decreasing = T),])
CCAN_counts<-as.data.frame(table(ccans_human$CCAN))
row.names(CCAN_counts)<-paste0("CCAN_CCAN_",CCAN_counts$Var1)

Human_ordered_ccan<-CCAN_counts[match(CCANs_human_unique,row.names(CCAN_counts)),]

#looking just at min and max

CARD_min<-which.min(s_h_ccan$Cardiomyocytes-s_m_ccan$Cardiomyocytes)

s_h_ccan[CARD_max,]
s_m_ccan[CARD_max,]

s_h_ccan[CARD_min,]
s_m_ccan[CARD_min,]

#CCAN_CCAN_27
#CCAN_29778

CARD_reg_sub<-StringToGRanges(ccans_human[which(ccans_human$CCAN==27),1])
CARD_reg<-StringToGRanges("chr1-101216546-101218131")
 

  m_cov<-CoveragePlot(
    object = mouse, 
    region = CARD_reg,assay ="cre",
    extend.upstream = 1000,
    extend.downstream = 1000,annotation=T,peaks=T, region.highlight = CARD_reg_sub
  )
  
h_cov<-CoveragePlot(
    object = human, 
    region = CARD_reg,assay ="cre",
    extend.upstream = 1000,
    extend.downstream = 1000,annotation=T,peaks=T, region.highlight = CARD_reg_sub
  )
  

CARD_reg_sub<-StringToGRanges(ccans_human[which(ccans_human$CCAN==29778),1])
CARD_reg<-StringToGRanges("chr6-72552870-72554535")
 

  m_cov2<-CoveragePlot(
    object = mouse, 
    region = CARD_reg,assay ="cre",
    extend.upstream = 1000,
    extend.downstream = 1000,annotation=T,peaks=T, region.highlight = CARD_reg_sub
  )
  
h_cov2<-CoveragePlot(
    object = human, 
    region = CARD_reg,assay ="cre",
    extend.upstream = 1000,
    extend.downstream = 1000,annotation=T,peaks=T, region.highlight = CARD_reg_sub
  )

#go thorugh all ccans and and create an hg38 granges object

ccan_all_regions=data.frame(chr=character(),start=double(),end=double(),ccan=character())
for (ccan in levels(as.factor(ccans_human$CCAN))){

ccan_peaks<-ccans_human[ccans_human$CCAN==ccan,1]
ccan_peaks_granges<-StringToGRanges(ccan_peaks)
globseqnames<-unique(levels(seqnames(ccan_peaks_granges)))
globstart<-min(start(ccan_peaks_granges))
globend<-max(end(ccan_peaks_granges)) 

ccan_region=data.frame(chr=globseqnames,start=globstart,end=globend,ccan=ccan)
ccan_all_regions<-rbind(ccan_all_regions,ccan_region)
}

write.table(ccan_all_regions,file="CCAN_human_start_end_hg38.bed",quote = F,row.names = F,col.names = F,sep = "\t")
#lifted over
ccan_all_regions_mm10<-read.table("/brahms/torkenczyk/H_M_storage/analysis//CCAN_human_start_end_mm10.bed",header = F)
ccan_all_regions_mm9<-read.table("/brahms/torkenczyk/H_M_storage/analysis//CCAN_human_start_end_mm9.bed",header = F)


master_ccan_df=data.frame(ccan=ccan_all_regions$ccan,region_hg38=paste(ccan_all_regions$chr,ccan_all_regions$start,ccan_all_regions$end,sep = "_"))
mm10=data.frame(ccan=ccan_all_regions_mm10$V4,region_mm10=paste(ccan_all_regions_mm10$V1,ccan_all_regions_mm10$V2,ccan_all_regions_mm10$V3,sep = "_"))
mm9=data.frame(ccan=ccan_all_regions_mm9$V4,region_mm9=paste(ccan_all_regions_mm9$V1,ccan_all_regions_mm9$V2,ccan_all_regions_mm9$V3,sep = "_"))

master_ccan_df<-merge(master_ccan_df,mm10,by="ccan")
master_ccan_df<-merge(master_ccan_df,mm9,by="ccan")
write.table(master_ccan_df,file="CCAN_human_start_end_hg38_mm9_mm10.txt",quote = F,row.names = F,col.names = T,sep = "\t")
master_ccan_df<-read.table("CCAN_human_start_end_hg38_mm9_mm10.txt",header = T)


names(CCAN_counts)<-c("ccan","number of peaks")
master_ccan_df<-merge(master_ccan_df,CCAN_counts,by="ccan")

#add in order
for (ct in names(s_h_ccan)){
print(ct)
s_m_diff<-s_h_ccan[[ct]]-s_m_ccan[[ct]]
CCANs_human_unique<-row.names(s_h_ccan[order(s_m_diff,decreasing = T),]-s_m_ccan[order(s_m_diff,decreasing = T),])
ccan_order=data.frame(ccan=as.numeric(str_split_fixed(CCANs_human_unique,n=3,pattern = "_")[,3]),ct=1:length(str_split_fixed(CCANs_human_unique,n=3,pattern = "_")[,3]))
names(ccan_order)<-c("ccan",ct)
master_ccan_df<-merge(master_ccan_df,ccan_order,by="ccan",all.x=T)


}

#with TSS approach
tss.positions_hg38 <- GetTSSPositions(ranges = Annotation(human))
DefaultAssay(mouse)<-"mm10"
#this is mm10 do not forget 
tss.positions_mouse <- GetTSSPositions(ranges = Annotation(mouse))
master_ccan_df$human_TSS_dist<-distanceToNearest(StringToGRanges(master_ccan_df$region_hg38,sep = c("_","_")),tss.positions_hg38)@elementMetadata$distance
master_ccan_df$mouse_TSS_dist<-distanceToNearest(StringToGRanges(master_ccan_df$region_mm10,sep = c("_","_")),tss.positions_mouse)@elementMetadata$distance

write.table(master_ccan_df,file="CCAN_human_start_end_hg38_mm9_mm10_lift02mathc.txt",quote = F,row.names = F,col.names = T,sep = "\t")

#now look at ordered cases
#CARD_order<-s_h_ccan$Cardiomyocytes-s_m_ccan$Cardiomyocytes
#CARD_reg_sub<-StringToGRanges(ccans_human[which(ccans_human$CCAN==35177),1])
#CARD_reg<-StringToGRanges("chr9-73445370-73584010")

master_ccan_df[order(master_ccan_df$Hepatocytes),2][2]

#larger CCAN less specific?

ggplot(master_ccan_df[master_ccan_df$human_TSS_dist==0,],aes(x=number.of.peaks,y=Enterocytes))+geom_point()

#not enough proof

A<-umap::umap(d = na.omit(((s_h_ccan-s_m_ccan+s_m_ccan))))
B<-umap::umap(d = na.omit(s_h_ccan-s_m_ccan+2*s_m_ccan-s_h_ccan))

Acomb<-s_h_ccan
Bcomb<-s_m_ccan

names(Acomb)<-paste0("human_",names(s_h_ccan))
names(Bcomb)<-paste0("mouse_",names(s_m_ccan))

AB<-cbind(Acomb,Bcomb)

AB_u<-umap::umap(d = na.omit(AB))

pmax(na.omit(s_h_ccan-s_m_ccan+s_m_ccan))

h_max<-apply(na.omit(s_h_ccan-s_m_ccan+s_m_ccan), 1, max)
h_max_names<-apply(na.omit(s_h_ccan-s_m_ccan+s_m_ccan), 1, FUN=function(x){names(s_h_ccan)[which.max(x)]})
m_max<-apply(na.omit(s_h_ccan-s_m_ccan+2*s_m_ccan-s_h_ccan), 1, max)
m_max_names<-apply(na.omit(s_h_ccan-s_m_ccan+2*s_m_ccan-s_h_ccan), 1, FUN=function(x){names(s_h_ccan)[which.max(x)]})




AB<-as.data.frame(AB_u$layout)

AB$h_max<-h_max
AB$m_max<-m_max
AB$h_max_names<-h_max_names
AB$m_max_names<-m_max_names
AB$ccan<-row.names(AB)


G<-data.frame(ccan=paste0("CCAN_CCAN_",master_ccan_df$ccan),number.of.peaks=master_ccan_df$number.of.peaks)
AB_merg<-merge(AB,G,by="ccan")


C<-cbind(as.data.frame(A$layout),na.omit(s_h_ccan-s_m_ccan+s_m_ccan))
C2<-cbind(as.data.frame(B$layout),na.omit(s_h_ccan-s_m_ccan+s_m_ccan))
D<-cbind(as.data.frame(B$layout),na.omit(s_h_ccan-s_m_ccan+2*s_m_ccan-s_h_ccan))
D2<-cbind(as.data.frame(A$layout),na.omit(s_h_ccan-s_m_ccan+2*s_m_ccan-s_h_ccan))

plot_list<-list()
for (ct in names(s_h_ccan)) {
print(ct)
p1h<-ggplot(C2,aes_string(x="V1",y="V2",color=paste0(ct)))+geom_point()+theme_classic()+ylab("UMAP2")+xlab("UMAP1")+theme(legend.position = "bottom")
p2m<-ggplot(D,aes_string(x="V1",y="V2",color=paste0(ct)))+geom_point()+theme_classic()+ylab("UMAP2")+xlab("UMAP1")+theme(legend.position = "bottom")
plot_list[[ct]]<-p1h|p2m

}

ph_names<-ggplot(AB,aes(x=V1,y=V2,color=h_max_names))+geom_point()+theme_classic()+ylab("UMAP2")+xlab("UMAP1")+theme(legend.position = "bottom")
pm_names<-ggplot(AB,aes(x=V1,y=V2,color=m_max_names))+geom_point()+theme_classic()+ylab("UMAP2")+xlab("UMAP1")+theme(legend.position = "bottom")

ph_vals<-ggplot(AB,aes(x=V1,y=V2,color=h_max))+geom_point()+theme_classic()+ylab("UMAP2")+xlab("UMAP1")+theme(legend.position = "bottom")
pm_vals<-ggplot(AB,aes(x=V1,y=V2,color=m_max))+geom_point()+theme_classic()+ylab("UMAP2")+xlab("UMAP1")+theme(legend.position = "bottom")

wrap_plots(list(ph_names,ph_vals,pm_names,pm_vals),ncol = 2,nrow = 2)

pnumb<-ggplot(AB_merg,aes(x=V1,y=V2,color=log10(number.of.peaks)))+geom_point()+theme_classic()+ylab("UMAP2")+xlab("UMAP1")+theme(legend.position = "bottom")

pnumb2<-ggplot(AB_merg,aes(x=h_max,y=number.of.peaks))+geom_point()+theme_classic()+xlab("human max score")+ylab("log10(number of peaks)")+theme(legend.position = "bottom")

#plot some of them

master_ccan_df<-read.delim(file = "/brahms/torkenczyk/H_M_storage/analysis//CCAN_human_start_end_hg38_mm9_mm10_lift02mathc.txt")
human<-readRDS("/brahms/torkenczyk/H_M_storage/analysis/objects/Human_added_links.rds")
mouse<-readRDS("/brahms/torkenczyk/H_M_storage/analysis/objects/Mouse_added_links.rds")
ccans_mouse_hg38<-read.delim(file="/home/torkenczyk/git_repos/H_M_analysis/Mouse_ccans_cicero_hg38.txt")
ccans_mouse_mm10<-read.delim(file="/home/torkenczyk/git_repos/H_M_analysis/Mouse_ccans_cicero_mm10.txt")
ccans_human<-read.delim(file="/brahms/torkenczyk/H_M_storage/analysis//Human_ccans_cicero.txt")

plot_CCAN_species<-function(number=2,ct="Hepatocytes",mincut=0.1,upstr=1000,dstr=1000)
{

human_ccan_select<-master_ccan_df[order(master_ccan_df[ct]),1][number]
human_region<- StringToGRanges(master_ccan_df[order(master_ccan_df[ct]),2][number],sep = c("_","_"))
mouse_region<-StringToGRanges(master_ccan_df[order(master_ccan_df[ct]),3][number],sep = c("_","_"))

human_highlight<-StringToGRanges(ccans_human$Peak[ccans_human$CCAN==human_ccan_select],sep = c("-","-"))
m_peaks<-StringToGRanges(ccans_mouse_mm10$Peak,sep = c("-","-"))
mouse_highlight<-m_peaks[queryHits(findOverlaps(query = m_peaks,subject = mouse_region))]

if (length(mouse_highlight)==0)
{
    m_cov1<-CoveragePlot(
    object = mouse, 
    region = mouse_region,assay ="mm10",
    extend.upstream = upstr,
    extend.downstream = dstr,annotation=T,peaks=T,links = FALSE
  )
} else {
  m_cov1<-CoveragePlot(
    object = mouse, 
    region = mouse_region,assay ="mm10",
    extend.upstream = upstr,
    extend.downstream = dstr,annotation=T,peaks=T,region.highlight = mouse_highlight,links = F
  )  
  
}
  
  

  
  mregion_extend<-FindRegion(
    object = mouse,
    region = mouse_region,
    assay ="mm10",
    extend.upstream = upstr,
    extend.downstream = dstr
  )
  
  mouse_link_mm10<-LinkPlot(object = mouse[["mm10"]], region = mregion_extend,min.cutoff = mincut)

  
  m_cov2<-CoveragePlot(
    object = mouse, 
    region =  human_region,assay ="cre",region.highlight = human_highlight,
    extend.upstream = upstr,
    extend.downstream = dstr,annotation=T,peaks=T,links = F
  )
  
  hregion_extend<-FindRegion(
    object = mouse,
    region = human_region,
    assay ="cre",
    extend.upstream = upstr,
    extend.downstream = dstr
  )
  
  mouse_link_hg38<-LinkPlot(object = mouse[["cre"]], region = hregion_extend,min.cutoff = mincut)

  
h_cov2<-CoveragePlot(
    object = human, 
    region =  human_region,assay ="cre",region.highlight = human_highlight,
    extend.upstream = upstr,
    extend.downstream = dstr,annotation=T,peaks=T,links = F
  )
  
human_link_hg38<-LinkPlot(object = human[["cre"]], region = hregion_extend,min.cutoff = mincut)

p<-wrap_plots(list(m_cov1,m_cov2,h_cov2,mouse_link_mm10,mouse_link_hg38,human_link_hg38),ncol = 3,nrow = 2)  
return(p)
}

C<-plot_CCAN_species(number=5,mincut = 0.1,ct = "Cardiomyocytes")
MY<-plot_CCAN_species(number=1,mincut = 0.1,ct = "Myofibroblasts")
B<-plot_CCAN_species(number=1,mincut = 0.1,ct = "B-lymphocytes")
E<-plot_CCAN_species(number=1,mincut = 0.1,ct = "Enterocytes")
Fi<-plot_CCAN_species(number=1,mincut = 0.1,ct = "Fibroblasts")
T<-plot_CCAN_species(number=2,mincut = 0.1,ct = "T-lymphocytes")
End<-plot_CCAN_species(number=1,mincut = 0.1,ct = "Endothelial cells")
Gob<-plot_CCAN_species(number=5,mincut = 0.15,ct = "Goblet.cells",upstr = 100000,dstr = 100000)
PN<-plot_CCAN_species(number=1,mincut = 0.15,ct = "Pneumocytes")
Sch<-plot_CCAN_species(number=1,mincut = 0.15,ct = "Schwann.cells")
M<-plot_CCAN_species(number=6,mincut = 0.15,ct = "Macrophages")
S<-plot_CCAN_species(number=1,mincut = 0.15,ct = "Smooth.muscle.cells")



png(filename = paste0("/brahms/torkenczyk/H_M_storage/analysis//plots/",gene,"_combined.wcons_yaxe_human_VAL=.png"),width = 20,height=28,units = "cm",res=300)



```

Might be easy to show GA across mouse and human as shared or unique

```{r}

human<-readRDS("/brahms/torkenczyk/H_M_storage/analysis//active_use_objects/Select_human_tissues.rds")
mouse<-readRDS("/brahms/torkenczyk/H_M_storage/analysis//active_use_objects/Select_mouse_tissues.rds")

gene.activities<-GeneActivity(object = human,assay = "cre")

# add the gene activity matrix to the Seurat object as a new assay
human[['RNA']] <- CreateAssayObject(counts = gene.activities)
human <- NormalizeData(
  object = human,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(human$nCount_RNA)
)

gene.activities<-GeneActivity(object = mouse,assay = "cre")

# add the gene activity matrix to the Seurat object as a new assay
mouse[['RNA']] <- CreateAssayObject(counts = gene.activities)
mouse <- NormalizeData(
  object = mouse,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(mouse$nCount_RNA)
)

DefaultAssay(mouse)<-"RNA"
DefaultAssay(human)<-"RNA"

H_genes<-FindAllMarkers(human)

M_genes<-FindAllMarkers(mouse)

sh_genes<-list()
genes_to_plot<-c()
for (ct in levels(as.factor(H_genes$cluster)))
{
H_pos<-H_genes[H_genes$avg_log2FC>0 & H_genes$cluster==ct,]
M_pos<-M_genes[M_genes$avg_log2FC>0 & M_genes$cluster==ct,]
sh_genes[[ct]]<-merge(H_pos,M_pos,by="gene")
genes_to_plot=c(genes_to_plot,head(sh_genes[[ct]]$gene,n=5))
}




p1<-DotPlot(human,features = unique(genes_to_plot))+theme(legend.position="top")+ggtitle(label = "Human")+xlab(label = "")+coord_flip()
p2<-DotPlot(mouse,features = unique(genes_to_plot))+theme()+theme(legend.position="top")+ggtitle(label = "Mouse") +theme(axis.title.y=element_blank(),axis.text.y=element_blank(), axis.ticks.y=element_blank())+xlab(label = "")



DotPlot(human,features = c("MYH7","MYH6","SHANK3","SPARC","MPRIP","CNN1","MYH11","ZFHX3","TAGLN","SEPT9","CARMN","MYH11"))
DotPlot(mouse,features = c("MYH7","MYH6","SHANK3","SPARC","MPRIP","CNN1","ZFHX3","TAGLN","SEPT9","CARMN","MYH11"))


```