<commit_before>library(slam)
library(pheatmap)
library(amap)
library(RColorBrewer)

setwd("~/Desktop/A1-2_Project/Brg1_ChIP/NS50200/Peaks/Peaks_by_JH/")
brg1LD<-read.table("motifs/LabBrg1_DexPeaks_motifAnnotate.txt", header=T, 
                   quote="", allowEscapes=F, stringsAsFactors=F, sep="\t")

colnames(brg1LD)<-gsub("\\.Distance\\.From.+","",colnames(brg1LD), perl=T)
colnames(brg1LD)[1]<-"PeakID"
brg1LDmotifs<-brg1LD[,-c(1:21)]
brg1LDmotifs<-apply(brg1LDmotifs,1, gsub, 
                    pattern="^-*(\\d{1,3}).+", replacement="\\1", perl=T)
brg1LDmotifs<-t(brg1LDmotifs)
class(brg1LDmotifs)<-"numeric"
brg1LDmotifs<-ifelse(is.na(brg1LDmotifs), 0, ifelse(brg1LDmotifs>250,1,0))
brg1LDmotifs<-as.simple_triplet_matrix(brg1LDmotifs)

brg1LDcoOccur<-crossprod_simple_triplet_matrix(brg1LDmotifs)
brg1Hclust<-hcluster(brg1LDcoOccur, "euclid")

png("LabBrg1_Dexonly_Motif_coOccurence_transform.png", height=2000, width=2500, units="px")
pheatmap(sqrt(brg1LDcoOccur),color=colorRampPalette(brewer.pal(9,"Greens"))(100),
         cluster_rows = brg1Hclust, cluster_cols = brg1Hclust)
dev.off()

tree<-cutree(brg1Hclust,k = 5)
treeHeat<-tree %in% c(1,3,4)

png("LabBrg1_Dexonly_Motif_coOccurence_subset.png", height=1500, width=2500, units="px")
pheatmap(brg1LDcoOccur[treeHeat,treeHeat],border_color = NA,
         color=colorRampPalette(brewer.pal(9,"Greens"))(100))
dev.off()


brg1Peaks<-t(brg1LDmotifs)
peaksCoOccur<-crossprod_simple_triplet_matrix(brg1Peaks)
peakClust<-hcluster(peaksCoOccur, "eucli",nbproc = 16)

png("Brg1_Peak_adj_matrix.png",heigh=2500, width=3500, units="px")
pheatmap(peaksCoOccur, border_color=NA,color=colorRampPalette(brewer.pal(9,"Greens"))(100))


<commit_msg>major update to speed motif cooccurence analysis from peak data<commit_after>library(reshape2)
library(igraph)
library(slam)
library(pheatmap)
library(amap)
library(RColorBrewer)

files<-dir(pattern=".txt")

annot<-lapply(files, function(x) read.table(x, header=T,quote="", 
                   allowEscapes=F, stringsAsFactors=F, sep="\t", fill=T))


names(annot)<-gsub("_motifAnnotations.txt", "", files)

for(i in seq_along(annot)){
    colnames(annot[[i]])<-gsub("\\.Distance\\.From.+","",colnames(annot[[i]]), perl=T)
    colnames(annot[[i]])[1]<-"PeakID"
}

annotmotifs<-lapply(annot, function(x) x[,-c(1:21)])
annotmotifs<-lapply(annotmotifs, function(x) {
                                apply(x,1, gsub, 
                                        pattern="\\([ACTG]+,[\\+-],\\d\\.\\d\\d\\)", replacement="", perl=T)
                })
annotmotifs<-lapply(annotmotifs, t)
annotmotifs<-lapply(annotmotifs, function(x) ifelse(x=="", 0,ifelse(is.na(x)==T, 0, 1)))

for(i in seq_along(annotmotifs)){
    class(annotmotifs[[i]])<-"numeric"
    annotmotifs[[i]]<-annotmotifs[[i]][,which(!is.na(col_sums(annotmotifs[[i]])))]
}
matrixD<-lapply( annotmotifs, function(x) {
    crossprod_simple_triplet_matrix(as.simple_triplet_matrix(!x))
})
# annotmotifs<-as.simple_triplet_matrix(annotmotifs)

annotcoOccur<-lapply( annotmotifs, function(x){
    crossprod_simple_triplet_matrix(as.simple_triplet_matrix(x)) })

annotJaccard<-lapply(seq_along(annot), function(x) {
    annotcoOccur[[x]]/(nrow(annot[[x]])-matrixD[[x]])
})

motifJaccardclust<-lapply(annotJaccard, hcluster,  method="euclid",nbproc = 8)
motifAdjMatClust<-lapply( annotcoOccur,  hcluster, method="euclid",nbproc = 8)

lapply(seq_along(annotcoOccur), function(x) {
    png(paste0(names(annot)[x],"_minOdds_Motif_coOccurence.png"), height=2000, width=2500, units="px")
        pheatmap(sqrt(annotcoOccur[[x]]),color=colorRampPalette(brewer.pal(9,"Greens"))(100),
            cluster_rows = motifAdjMatClust[[x]], cluster_cols = motifAdjMatClust[[x]])
    dev.off()

    png(paste0(names(annot)[x],"_minOdds_Motif_Jaccard.png"), height=2000, width=2500, units="px")
        pheatmap(annotJaccard[[x]],color=colorRampPalette(brewer.pal(9,"Greens"))(100),
            cluster_rows = motifJaccardclust[[x]], cluster_cols = motifJaccardclust[[x]])
    dev.off()
})
###Remove Motifs that have perfect concordance...

edges<-lapply(annotJaccard, function(x){ 
    e<-melt(x)
    e[as.integer(e$Var2)>as.integer(e$Var1),]
    })

nodesToRemove<-lapply(seq_along(edges), function(x)
    rownames(annotJaccard[[x]]) %in% edges[[x]][edges[[x]]$value==1,2]
)
annotJaccard2<-mapply(x=annotJaccard, y=nodesToRemove,
                      function(x,y) x[!y,!y]
)
motifJaccardclust<-lapply(annotJaccard2, hcluster,method="euclid",nbproc = 8)
motifAdjMatClust<-lapply(annotJaccard2, hcluster, method="euclid",nbproc = 8)

lapply(seq_along(annotJaccard2), function(x) {
    png(paste0(names(annot)[x],"_minOdds_Motif_Jaccard2.png"), height=2000, width=2500, units="px")
    pheatmap(annotJaccard2[[x]],color=colorRampPalette(brewer.pal(9,"Greens"))(100),
             cluster_rows = motifJaccardclust[[x]], cluster_cols = motifJaccardclust[[x]])
    dev.off()
})


edgelist1<-lapply(annotJaccard2, function(x){
    e<-melt(x)
    e[as.integer(e$Var1)<as.integer(e$Var2),]
})

edgesIDs<-lapply(edgelist1, function(x)
    paste0(as.character(x$Var1),"_",as.character(x$Var2)))

edgesUnion<-(edgesIDs[[1]] %in% edgesIDs[[2]]) & (edgesIDs[[1]] %in% edgesIDs[[3]])

edgeDF<-lapply(edgesIDs, function(x) x %in% edgesIDs[[1]][edgesUnion])
edgesDF<-mapply(x=edgelist1, y=edgeDF, function(x,y) x[y,], SIMPLIFY = F)

edgeDF<-do.call(cbind, edgesDF)
edgeDF<-edgeDF[,-c(4,5,7,8)]
colnames(edgeDF)<-c("to", "from",names(annot))

kclust<-kmeans(edgeDF[,3:5],centers = 8)
ids<-lapply(1:8, function(x) sample(which(kclust$cluster==x),100))
ids<-unlist(ids)

png("Motif_agg1.png", height=3000, width=1000, units="px")
pheatmap(edgeDF[ids,3:5], annotation_row = edgeDF[ids,1:2],border_color = NA,
         cluster_rows = F, cluster_cols = F,
         color=colorRampPalette(brewer.pal(9,"Greens"))(100))
dev.off()

png("Motif_agg2.png", height=3000, width=1000, units="px")
    pheatmap(edgeDF[,3:5], annotation_row = edgeDF[,1:2],k_means=kclust,
             cluster_rows = F, cluster_cols = F,border_color=NA,
             color=colorRampPalette(brewer.pal(9,"Greens"))(100))
dev.off()

edgelist<-lapply(edgelist1, function(x)
    x[x$value>quantile(x$value,0.9),])

for(i in seq_along(edgelist)){
    colnames(edgelist[[i]])<-c("from","to","Weight")
}

graph<-lapply(edgelist, graph_from_data_frame ,directed=FALSE)

nodeSize<-lapply(seq_along(annotJaccard2), function(x){
    col_sums(annotmotifs[[x]])[colnames(annotmotifs[[x]]) %in% V(graph[[x]])$name]
})

for(i in seq_along(graph)){
    V(graph[[i]])$Size<-nodeSize[[i]]
}

lapply(seq_along(graph), function(x) {
    write_graph(graph[[x]], file=paste0(names(annotmotifs)[x],"_network.graphml")
                ,format="graphml")
    })
graphLay<-layout_with_fr(graph,weights = 1-edgelist$Weight)
graphLay<-layout_with_kk(graph,weights = -log(edgelist$Weight))
graphLay<-layout_with_dh(graph)

png("motif_network.png", height=2500, width=2500, units="px")
plot(graph, vertex.size=2)
dev.off()

png("motif_network_v7.png", height=2500, width=2500, units="px")
plot(graph, vertex.size=2, layout=graphLay)
dev.off()

{'fraction_non_alphanumeric': 0.19134845709454698, 'fraction_numerical': 0.0342398196421023, 'mean_word_length': 6.220752797558494, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6569112', 'n_tokens_mistral': 3124, 'n_tokens_neox': 2889, 'n_words': 304}