Skip to content

Instantly share code, notes, and snippets.

@cplaisier
Created May 8, 2020 17:42
Show Gist options
  • Select an option

  • Save cplaisier/35aa2a89c67dfdf52fdf1d5653be0bc4 to your computer and use it in GitHub Desktop.

Select an option

Save cplaisier/35aa2a89c67dfdf52fdf1d5653be0bc4 to your computer and use it in GitHub Desktop.
Example code to conduct AM and MDM analyses:
# Set working directory
#setwd('C:/Users/cplaisie/Dropbox/Urdahl/Ag85B_ESAT6_ts/network')
#setwd('/home/cplaisier/Dropbox (ASU)/Urdahl/macs/cMonkey')
setwd('C:/Users/plais/Dropbox (ASU)/Urdahl/macs/cMonkey')
#setwd('C:/Users/Christopher/Dropbox/Urdahl/Ag85B_ESAT6_ts/network')
# Get eigengenes
d1 = read.csv('eigengenes/biclusterEigengenes_tfbs_db.csv',header=T, row.names=1)
rownames(d1) = paste('TFBS_DB_',rownames(d1),sep='')
d2 = read.csv('eigengenes/biclusterEigengenes_pita.csv',header=T, row.names=1)
rownames(d2) = paste('PITA_',rownames(d2),sep='')
d3 = read.csv('eigengenes/biclusterEigengenes_targetscan.csv',header=T, row.names=1)
rownames(d3) = paste('TargetScan_',rownames(d3),sep='')
dAll = rbind(d1,d2,d3)
# Get expression data
ge1 = read.table('../gexp_cMonkey_norm_cv8000_preProcessed.tsv', header=T, row.names=1, sep='\t')
# Ordering samples
#os1 = c(c(1:4,9:12,16:19,23:25),c(5:8,13:15,20:22,26:28))
#ge2 = ge1[,os1]
#dAll2 = dAll[,os1]
ge2 = ge1
dAll2 = dAll[,colnames(ge1)]
# Plot Bicluster vs. conditions
library(pheatmap)
library(RColorBrewer)
library(gplots)
sets = c(rep(1,4),rep(2,3),rep(3,4),rep(4,3))
names1 = c(rep('Uninf_AM',4),rep('Inf_AM',3),rep('Uninf_Mono',4),rep('Inf_Mono',4))
#tfs1 = c('uc007ldr.2', rep('uc007xlk.1',2),rep('uc009orz.1',7),rep('uc007ifo.1',2))
#tfs2 = c(c('TargetScan_331'), c('TargetScan_82', 'PITA_56'), c('TargetScan_2','PITA_241','PITA_312','TFBS_DB_169','PITA_234','TFBS_DB_274','PITA_47'),c('TargetScan_188','PITA_334'))
tfNames2 = tfs2
pdf('tfBiclusters_noScale_lines.pdf')
#par(mar=c(7,4,4,2),mfrow=c(1,2))
for(tf2 in 1:length(tfs2)) {
par(mar=c(7,4,4,2),mfrow=c(1,2))
tmp2 = t(sapply(1:max(sets), function(x) { mean(as.numeric(ge2[tfs1[tf2],which(sets==x)])) } ))
sdTmp2 = t(sapply(1:max(sets), function(x) { sd(as.numeric(ge2[tfs1[tf2],which(sets==x)])) } ))
tmp1 = t(sapply(1:max(sets), function(x) { mean(as.numeric(dAll2[tfs2[tf2],which(sets==x)])) } ))
sdTmp1 = t(sapply(1:max(sets), function(x) { sd(as.numeric(dAll2[tfs2[tf2],which(sets==x)])) } ))
colnames(tmp1) = names1
ylim1 = c(min(c(tmp1,tmp2)),max(c(tmp1,tmp2)))
plotCI(y=tmp1[1:4],x=c(0,2,4,8),uiw=sdTmp1[o1],gap=0,ylab='Bicluster Relative Expression',pch=19,col='black',main=paste(tfs1[tf2],' vs. ',tfs2[tf2],'\nESAT6',sep=''),xlab='Time Post Stim.',ylim=ylim1)
lines(c(0,2,4,8),tmp1[1:4],lwd=1,col=rgb(0,0,0,0.5),lty=2)
par(new=T)
plotCI(y=tmp2[1:4],x=c(0,2,4,8),uiw=sdTmp2[o1],gap=0,ylab='',xaxt='n',xlab='',pch=19,col='red',yaxt='n')
lines(c(0,2,4,8),tmp2[1:4],lwd=2,col=rgb(1,0,0,0.5))
plotCI(y=tmp1[5:8],x=c(0,2,4,8),uiw=sdTmp1[o1],gap=0,ylab='Bicluster Relative Expression',pch=19,col=rgb(0,0,0,1),main=paste(tfs1[tf2],' vs. ',tfs2[tf2],'\nAg85B',sep=''),xlab='Time Post Stim.',ylim=ylim1)
lines(c(0,2,4,8),tmp1[5:8],lwd=1,col=rgb(0,0,0,0.5),lty=2)
par(new=T)
plotCI(y=tmp2[5:8],x=c(0,2,4,8),uiw=sdTmp2[o1],gap=0,ylab='',xaxt='n',xlab='',pch=19,col=rgb(1,0,0,1),yaxt='n')
lines(c(0,2,4,8),tmp2[5:8],lwd=2,col=rgb(1,0,0,0.5))
}
dev.off()
#### GOOD ONE ####
library(pheatmap)
library(RColorBrewer)
sets1 = c(rep(1,4),rep(2,4),rep(3,4),rep(4,3),rep(5,4),rep(6,3),rep(7,3),rep(8,3))
names1 = c('ESAT6_0H','ESAT6_2H','ESAT6_4H','ESAT6_8H','Ag85B_0H','Ag85B_2H','Ag85B_4H','Ag85B_8H')
tfs1 = c('uc007ldr.2', rep('uc007xlk.1',2),rep('uc009orz.1',7),rep('uc007ifo.1',2),rep('uc012gkc.1',3))
tfs2 = c(c('TargetScan_331'), c('TargetScan_82', 'PITA_56'), c('TargetScan_2','PITA_241','PITA_312','TFBS_DB_169','PITA_234','TFBS_DB_274','PITA_47'),c('TargetScan_188','PITA_334'),c('TFBS_DB_71','PITA_50','TFBS_DB_309'))
tfNames2 = tfs2
pdf('tfBiclusters_noScale.pdf')
#par(mar=c(7,4,4,2))
for(tf2 in 1:length(tfs2)) {
par(mfrow=c(2,1),mar=c(7,4,4,2))
tmp2 = t(sapply(1:max(sets1), function(x) { mean(as.numeric(ge2[tfs1[tf2],which(sets1==x)])) } ))
colnames(tmp2) = names1
c1 = cor.test(as.numeric(ge2[tfs1[tf2],]),as.numeric(dAll2[tfs2[tf2],]))
c2_ESAT6 = cor.test(as.numeric(ge2[tfs1[tf2],])[1:15],as.numeric(dAll2[tfs2[tf2],])[1:15])
c2_Ag85B = cor.test(as.numeric(ge2[tfs1[tf2],])[16:28],as.numeric(dAll2[tfs2[tf2],])[16:28])
boxplot(as.numeric(dAll2[tfs2[tf2],]) ~ sets1, names=names1, las=2,main=paste(tfNames2[tf2],':',tfs1[tf2],'\nR.all = ',round(c1$estimate,2),'; R.ESAT6 = ',round(c2_ESAT6$estimate,2),'; R.Ag85B = ',round(c2_Ag85B$estimate,2)), col=brewer.pal(8, "Paired"))
par(new=T)
plot(1:4,tmp2[1:4],col=rgb(1,0,0,0.3),type='l',lty=1,lwd=2,axes=F,xlim=c(0.5,8.5),ylab='',xlab='')
par(new=T)
points(1:4,tmp2[1:4],col=rgb(1,0,0,0.8),pch=19,xlim=c(0.5,8.5),ylab='',xlab='')
par(new=T)
plot(5:8,tmp2[5:8],col=rgb(1,0,0,0.3),type='l',lty=1,lwd=2,axes=F,xlim=c(0.5,8.5),ylab='',xlab='')
par(new=T)
points(5:8,tmp2[5:8],col=rgb(1,0,0,0.8),pch=19,xlim=c(0.5,8.5),ylab='',xlab='')
abline(v=4.5,col='black')
boxplot(t(as.matrix(ge2[biclustMembership[[tfs2[tf2]]],]))~sets1,names=names1,las=2,main=paste('Genes: ',tfNames2[tf2]),col=brewer.pal(8, "Paired"))
tmp1 = t(sapply(1:max(sets1), function(x) { rowMeans(ge2[biclustMembership[[tfs2[tf2]]],which(sets1==x)]) } ))
rownames(tmp1) = names1
matlines(1:4,tmp1[1:4,],type='l',col=rgb(0.3,0.3,0.3,0.5),lty=1,lwd=1)
matlines(5:8,tmp1[5:8,],type='l',col=rgb(0.3,0.3,0.3,0.5),lty=1,lwd=1)
par(new=T)
plot(1:4,tmp2[1:4],col=rgb(1,0,0,0.3),type='l',lty=1,lwd=2,axes=F,xlim=c(0.5,8.5),ylab='',xlab='')
par(new=T)
points(1:4,tmp2[1:4],col=rgb(1,0,0,0.8),pch=19,xlim=c(0.5,8.5),ylab='',xlab='')
par(new=T)
plot(5:8,tmp2[5:8],col=rgb(1,0,0,0.3),type='l',lty=1,lwd=2,axes=F,xlim=c(0.5,8.5),ylab='',xlab='')
par(new=T)
points(5:8,tmp2[5:8],col=rgb(1,0,0,0.8),pch=19,xlim=c(0.5,8.5),ylab='',xlab='')
abline(v=4.5,col='black')
pheatmap(t(tmp1),cluster_cols=F,scale='row')
}
dev.off()
# Select out clusters of
library(GMD)
tm1 = sapply(unique(sets), function(x) { rowMeans(dAll[,which(sets==x)]) })
hc1 = hclust(dist(dAll))
css1 = css.hclust(dist(dAll),hc1)
#hc1 = hclust(dist(cor(t(d1[,names(col.rug)[f1]]))))
#css1 = css.hclust(dist(cor(t(d1[,names(col.rug)[f1]]))),hc1)
elbow1 = elbow.batch(css1)
elbow1 = elbow.batch(css1,inc.thres=0.1,ev.thres=c(0.9,0.8))
cut1 = cutree(hc1, k=elbow1$k)
#cut1 = cutree(hc1, k=5)
#write.csv(cbind('Cluster Index'=cut1),'clusters.csv')
#
tm1 = sapply(unique(sets), function(x) { rowMeans(dAll[,which(sets==x)]) })
hc1 = hclust(dist(tm1))
css1 = css.hclust(dist(tm1),hc1)
#hc1 = hclust(dist(cor(t(d1[,names(col.rug)[f1]]))))
#css1 = css.hclust(dist(cor(t(d1[,names(col.rug)[f1]]))),hc1)
elbow1 = elbow.batch(css1)
elbow1 = elbow.batch(css1,inc.thres=0.1,ev.thres=c(0.9,0.8))
cut1 = cutree(hc1, k=elbow1$k)
library(pheatmap)
#ann1 = data.frame(Clusters=as.factor(cut2), Subtypes=Subtypes)
#ann1 = data.frame(Subtypes=Subtypes)
#rownames(ann1) = names(cut2)
#Clusters = c('orange','red','black','cadetblue')
#names(Clusters) = c(1,2,3,4)
#SubColors = c('black','orange','cadetblue','red','green','white')
#names(SubColors) = c('classical','mesenchymal','neural','proneural','control','Unknown')
#ann_colors = list(Clusters=Clusters,Subtypes=SubColors)
#ann_colors = list(Subtypes=SubColors)
#pheatmap(d2[c('tfbs_db_230','pita_157','targetscan_530','hsa_miR_34a','hsa_miR_210','CEBPE','ETV6','FOSL1','IRF1','MAFF','NFKB1','TSC22D1'),names(cut2)[order(Subtypes)]],color=colorpanel(256,'blue','black','yellow'),annotation=ann1,show_rownames=T,show_colnames=F, annotation_colors=ann_colors,scale='row',cluster_cols=F)
#pheatmap(dAll,color=colorpanel(256,'blue','black','yellow'),annotation=ann1,show_rownames=T,show_colnames=F, annotation_colors=ann_colors,scale='row',cluster_cols=F)
pheatmap(dAll,color=colorpanel(256,'blue','black','yellow'),show_rownames=F,show_colnames=F,cluster_cols=F)
tm1 = sapply(unique(sets), function(x) { rowMeans(dAll[,which(sets==x)]) })
pheatmap(tm1,color=colorpanel(256,'blue','black','yellow'),show_rownames=F,show_colnames=F,cluster_cols=F)
m3 = t(sapply(1:length(unique(cut1)), function(x) {apply(dAll[names(which(cut1==x)),],2,mean)}))
heatmap.2(as.matrix(m3),trace='none',col=colorpanel(256,'blue','black','yellow'),Colv=F,Rowv=T,dendrogram='row',density.info='none') #,ColSideColors=as.character(col.rug),scale='row')
#dev.off()
library(miscTools)
tm3 = t(sapply(unique(cut1), function(x) { colMedians(tm1[names(which(cut1==x)),]) } ))
heatmap.2(as.matrix(tm3),trace='none',col=colorpanel(256,'blue','black','yellow'),Colv=F,Rowv=T,dendrogram='row',density.info='none') #,ColSideColors=as.character(col.rug),scale='row')
o1 = order(table(cut1))
#o1 = order(tm3[,3],decreasing=T)
heatmap.2(as.matrix(m3)[o1,],trace='none',col=colorpanel(256,'blue','black','yellow'),Colv=F,Rowv=F,dendrogram='row',density.info='none') #,ColSideColors=as.character(col.rug),scale='row')
# Resorting clusters
cut2 = cut1
cut2[which(cut1==o1[1])] = 1
cut2[which(cut1==o1[2])] = 2
cut2[which(cut1==o1[3])] = 3
cut2[which(cut1==o1[4])] = 4
cut2[which(cut1==o1[5])] = 5
cut2[which(cut1==o1[6])] = 6
cut2[which(cut1==o1[7])] = 7
cut2[which(cut1==o1[8])] = 8
cut2[which(cut1==o1[9])] = 9
cut2[which(cut1==o1[10])] = 10
cut2[which(cut1==o1[11])] = 11
cut2[which(cut1==o1[12])] = 12
cut2[which(cut1==o1[13])] = 13
cut2[which(cut1==o1[14])] = 14
# Write out bicluster memberhips
write.csv(cut2,'clustersOfBiclusters_VALL.csv')
########################################
# Heatmap of profiles for each cluster #
########################################
m3 = t(sapply(1:14, function(x) {apply(dAll[names(which(cut2==x)),],2,mean)}))
pdf('tCellNetworkClustersOfBiclusters_trends_VALL.pdf')
heatmap.2(as.matrix(m3),trace='none',col=colorpanel(256,'blue','black','yellow'),Colv=F,Rowv=F,dendrogram='none',RowSideColors=rainbow(length(unique(cut2))),density.info='none') #,scale='row')
dev.off()
#######################
## Getting Gene Sets ##
#######################
g1 = read.csv('genesInBiclusters/cluster.members.genes_tfbs_db.txt',header=F)
biclustMembership = list()
for(j in 1:length(g1[,1])) {
biclustMembership[[paste('TFBS_DB_',j,sep='')]] = strsplit(as.character(g1[j,]),split=' ')[[1]][-1]
}
g2 = read.csv('genesInBiclusters/cluster.members.genes_pita.txt',header=F)
for(j in 1:length(g2[,1])) {
biclustMembership[[paste('PITA_',j,sep='')]] = strsplit(as.character(g2[j,]),split=' ')[[1]][-1]
}
g3 = read.csv('genesInBiclusters/cluster.members.genes_targetScan.txt',header=F)
for(j in 1:length(g3[,1])) {
biclustMembership[[paste('TargetScan_',j,sep='')]] = strsplit(as.character(g3[j,]),split=' ')[[1]][-1]
}
# Take out
clusterMembership = list()
for(i in 1:18) {
tmp = c()
for(j in names(cut2[which(cut2==i)])) {
tmp = c(tmp,biclustMembership[[j]])
}
clusterMembership[[i]] = unique(tmp)
}
# Write out clusterMembership
dClustMemb = cbind(meta_cluster=1:length(clusterMembership),genes=unlist(lapply(clusterMembership,paste,collapse=';')))
write.csv(dClustMemb,'metaClusterMembership.csv',row.names=F)
# Number of genes per cluster
# 105 686 602 376 48 297 247 1118 470 929
sapply(1:length(clusterMembership), function(x) { length(clusterMembership[[x]]) })
# Functional enrichment
library(topGO)
# Read in GO mappings to affymetrix probe ids
library('org.Mm.eg.db')
# Make a map of Entrez IDs to UCSC transcript IDs
d1 = read.csv('entrez2ucsc.csv',header=F)
ucsc2entrez = list()
for(j in 1:length(d1[,1])) {
tmp = strsplit(as.character(d1[j,2]),split=';')[[1]]
for(k in tmp) {
ucsc2entrez[[as.character(k)]] = as.character(d1[j,1])
}
}
# Convert IDs to entrez
biclustMembership.eg = list()
for(j in 1:length(clusterMembership)) {
biclustMembership.eg[[j]] = as.character(unlist(ucsc2entrez[clusterMembership[[j]]]))
}
biclustMembership.eg = clusterMembership
# Setup for analysis
xx <- annFUN.org("BP", mapping = "org.Mm.eg.db", ID = "entrez")
#geneNames <- unique(unlist(xx))
geneNames <- intersect(unique(unlist(biclustMembership.eg)), unique(unlist(xx)))
tmp1 <- biclustMembership.eg[[1]]
geneList <- factor(as.integer(geneNames %in% tmp1))
names(geneList) <- geneNames
# Make Biological Process GOData object
GOdata.BP <- new("topGOdata", ontology='BP', allGenes = geneList, annot = annFUN.org, mapping = 'org.Mm.eg.db', ID = 'entrez')
m1.BP <- matrix(nrow = length(GOdata.BP@graph@nodes), ncol = length(biclustMembership.eg), dimnames = list(GOdata.BP@graph@nodes, 1:length(biclustMembership.eg)))
m2.BP <- matrix(nrow = length(biclustMembership.eg), ncol=2, dimnames = list(1:length(biclustMembership.eg),c('Top10.Terms.BP','BH.sig.GO.Ids.BP')))
# Make Molecular Function GOData object
GOdata.MF <- new("topGOdata", ontology='MF', allGenes = geneList, annot = annFUN.org, mapping = 'org.Mm.eg.db')
m1.MF <- matrix(nrow = length(GOdata.MF@graph@nodes), ncol = length(biclustMembership.eg), dimnames = list(GOdata.MF@graph@nodes, 1:length(biclustMembership.eg)))
m2.MF <- matrix(nrow = length(biclustMembership.eg), ncol=2, dimnames = list(1:length(biclustMembership.eg),c('Top10.Terms.MF','BH.sig.GO.Ids.MF')))
# Make Cellular Component GOData object
GOdata.CC <- new("topGOdata", ontology='CC', allGenes = geneList, annot = annFUN.org, mapping = 'org.Mm.eg.db')
m1.CC <- matrix(nrow = length(GOdata.CC@graph@nodes), ncol = length(biclustMembership.eg), dimnames = list(GOdata.CC@graph@nodes, 1:length(biclustMembership.eg)))
m2.CC <- matrix(nrow = length(biclustMembership.eg), ncol=2, dimnames = list(1:length(biclustMembership.eg),c('Top10.Terms.CC','BH.sig.GO.Ids.CC')))
for( clust in (1:length(biclustMembership.eg)) ) {
# Expand gene list and change factor in GOdata
clusterGenes <- biclustMembership.eg[[clust]]
GOdata.BP@allScores <- factor(as.integer(geneNames %in% clusterGenes))
GOdata.MF@allScores <- factor(as.integer(geneNames %in% clusterGenes))
GOdata.CC@allScores <- factor(as.integer(geneNames %in% clusterGenes))
# Biological process
r1.BP = runTest(GOdata.BP, algorithm = 'classic', statistic = 'fisher')
m1.BP[,clust] = r1.BP@score
m2.BP[clust,1] = gsub(',','_',paste(GenTable(GOdata.BP, r1.BP)[,2],collapse=';'))
m2.BP[clust,2] = paste(names(which(p.adjust(r1.BP@score,method='BH')<=0.05)),collapse=';')
#m2.BP[clust,2] = paste(names(which(r1.BP@score<=0.05)),collapse=';')
r1.MF = runTest(GOdata.MF, algorithm = 'classic', statistic = 'fisher')
m1.MF[,clust] = r1.MF@score
m2.MF[clust,1] = paste(GenTable(GOdata.MF, r1.MF)[,2],collapse=':')
m2.MF[clust,2] = paste(names(which(p.adjust(r1.MF@score,method='BH')<=0.05)),collapse=';')
#m2.MF[clust,2] = paste(names(which(r1.MF@score<=0.05)),collapse=';')
r1.CC = runTest(GOdata.CC, algorithm = 'classic', statistic = 'fisher')
m1.CC[,clust] = r1.CC@score
m2.CC[clust,1] = paste(GenTable(GOdata.CC, r1.CC)[,2],collapse=';')
m2.CC[clust,2] = paste(names(which(p.adjust(r1.CC@score,method='BH')<=0.05)),collapse=';')
#m2.CC[clust,2] = paste(names(which(r1.CC@score<=0.05)),collapse=';')
}
write.csv(m2.BP,'clusterEnrichment_GOBP_0.05.csv')
write.csv(m2.MF,'clusterEnrichment_GOMF_0.05.csv')
write.csv(m2.CC,'clusterEnrichment_GOCC_0.05.csv')
# Read in clusters of biclusters
clusters = {}
inFile = open('clustersOfBiclusters_VALL.csv','r')
inFile.readline() # Get rid of header
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.strip().split(',')
if not splitUp[1] in clusters:
clusters[splitUp[1]] = []
clusters[splitUp[1]].append(splitUp[0].strip('"').lower())
inFile.close()
# Read in post processed files
#Bicluster - 0 Genes - 1 Conditions - 2 Norm. Residual - 3 Norm. Residual Perm. P-Value - 4 Var. Exp. First PC - 5 Var. Exp. First PC Perm. P-Value - 6 MEME Motif1 E-Value - 7 Up.MEME Motif1 Consensus - 8 Up.MEME Motif1 Matches - 9 Up.MEME Motif1 Expanded Matches - 10 Up.MEME Motif1 Correlated Matches - 11 Up.MEME Motif1 Original/Expanded - 12 Up.MEME Motif1 Minimum Correlated - 13 Up.MEME Motif2 E-Value - 14 Up.MEME Motif2 Consensus - 15 Up.MEME Motif2 Matches - 16 Up.MEME Motif2 Expanded Matches - 17 Up.MEME Motif2 Correlated Matches - 18 Up.MEME Motif2 Original/Expanded - 19 Up.MEME Motif2 Minimum Correlated - 20 Up.WEEDER Motif1 Score - 21 Up.WEEDER Motif1 Consensus - 22 Up.WEEDER Motif1 Matches - 23 Up.WEEDER Motif1 Expanded Matches - 24 Up.WEEDER Motif1 Correlated Matches - 25 Up.WEEDER Motif1 Original/Expanded - 26 Up.WEEDER Motif1 Minimum Correlated - 27 Up.WEEDER Motif2 Score - 28 Up.WEEDER Motif2 Consensus - 29 Up.WEEDER Motif2 Matches - 30 Up.WEEDR Motif2 Expanded Matches - 31 Up.WEEDER Motif2 Correlated Matches - 32 Up.WEEDER Motif2 Original/Expanded - 33 Up.WEEDER Motif2 Minimum Correlated - 34 TFBS_DB.TFs - 35 TFBS_DB.percTargets - 36 TFBS_DB.pValue - 37 TFBS_DB.Exapnded Matches - 38 TFBS_DB.Correlated Matches - 39 TFBS_DB.Original/Expanded - 40 TFBS_DB.Minimum Correlated - 41 3pUTR.WEEDER Motif1 E-Value - 42 3pUTR.WEEDER Motif1 Consensus - 43 3pUTR.WEEDER Motif1 Matches - 44 3pUTR.WEEDER Motif1 Model - 45 3pUTR.WEEDER Motif2 E-Value - 46 3pUTR.WEEDER Motif2 Consensus - 47 3pUTR.WEEDER Motif2 Matches - 48 3pUTR.WEEDER Motif2 Model - 49 3pUTR_pita.miRNAs - 50 3pUTR_pita.percTargets - 51 3pUTR_pita.pValue - 52 3pUTR_targetScan.miRNAs - 53 3pUTR_targetScan.percTargets - 54 3pUTR_targetScan.pValue - 55 GO_Term_BP - 56
biclusters = {}
biclusterHeader = ''
for i in ['pita','targetscan','tfbs_db']:
inFile = open('postProcessed/postProcessed_urdahl_alvMac_'+i+'.csv','r')
biclusterHeader = inFile.readline() # Get rid of header
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.split(',')
biclusters[i+'_'+splitUp[0]] = line.strip()
inFile.close()
inFile = open('genesInBiclusters/cluster.members.genes_'+i+'.txt','r')
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.strip().split(' ')
j = splitUp.pop(0)
biclusters[i+'_'+j] = biclusters[i+'_'+j] + ',' + ';'.join(splitUp)
inFile.close()
inFile = open('genesInBiclusters/cluster.members.conditions_'+i+'.txt','r')
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.strip().split(' ')
j = splitUp.pop(0)
biclusters[i+'_'+j] = biclusters[i+'_'+j] + ',' + ';'.join(splitUp)
inFile.close()
# Dump combined postProcessed file with clusters IDs
#Bicluster - 0 Genes - 1 Conditions - 2 Norm. Residual - 3 Norm. Residual Perm. P-Value - 4 Var. Exp. First PC - 5 Var. Exp. First PC Perm. P-Value - 6 MEME Motif1 E-Value - 7 Up.MEME Motif1 Consensus - 8 Up.MEME Motif1 Matches - 9 Up.MEME Motif1 Expanded Matches - 10 Up.MEME Motif1 Correlated Matches - 11 Up.MEME Motif1 Original/Expanded - 12 Up.MEME Motif1 Minimum Correlated - 13 Up.MEME Motif2 E-Value - 14 Up.MEME Motif2 Consensus - 15 Up.MEME Motif2 Matches - 16 Up.MEME Motif2 Expanded Matches - 17 Up.MEME Motif2 Correlated Matches - 18 Up.MEME Motif2 Original/Expanded - 19 Up.MEME Motif2 Minimum Correlated - 20 Up.WEEDER Motif1 Score - 21 Up.WEEDER Motif1 Consensus - 22 Up.WEEDER Motif1 Matches - 23 Up.WEEDER Motif1 Expanded Matches - 24 Up.WEEDER Motif1 Correlated Matches - 25 Up.WEEDER Motif1 Original/Expanded - 26 Up.WEEDER Motif1 Minimum Correlated - 27 Up.WEEDER Motif2 Score - 28 Up.WEEDER Motif2 Consensus - 29 Up.WEEDER Motif2 Matches - 30 Up.WEEDR Motif2 Expanded Matches - 31 Up.WEEDER Motif2 Correlated Matches - 32 Up.WEEDER Motif2 Original/Expanded - 33 Up.WEEDER Motif2 Minimum Correlated - 34 TFBS_DB.TFs - 35 TFBS_DB.percTargets - 36 TFBS_DB.pValue - 37 TFBS_DB.Exapnded Matches - 38 TFBS_DB.Correlated Matches - 39 TFBS_DB.Original/Expanded - 40 TFBS_DB.Minimum Correlated - 41 3pUTR.WEEDER Motif1 E-Value - 42 3pUTR.WEEDER Motif1 Consensus - 43 3pUTR.WEEDER Motif1 Matches - 44 3pUTR.WEEDER Motif1 Model - 45 3pUTR.WEEDER Motif2 E-Value - 46 3pUTR.WEEDER Motif2 Consensus - 47 3pUTR.WEEDER Motif2 Matches - 48 3pUTR.WEEDER Motif2 Model - 49 3pUTR_pita.miRNAs - 50 3pUTR_pita.percTargets - 51 3pUTR_pita.pValue - 52 3pUTR_targetScan.miRNAs - 53 3pUTR_targetScan.percTargets - 54 3pUTR_targetScan.pValue - 55 GO_Term_BP - 56
tfs = {}
correlatedTfs = {}
bestTfs = {}
miRNAs = {}
goBP = {}
writeMe = ['Meta-Cluster,bicluster,'+biclusterHeader.strip()+',Gene Names,Condition Names']
for i in clusters:
tfs[i] = {'MEME':[], 'WEEDER':[], 'TFBS_DB':[]}
correlatedTfs[i] = {'MEME':[], 'WEEDER':[], 'TFBS_DB':[]}
bestTfs[i] = {'MEME':[], 'WEEDER':[], 'TFBS_DB':[]}
miRNAs[i] = []
goBP[i] = []
for j in clusters[i]:
writeMe.append(i+','+j+','+biclusters[j])
# Process out tfs and miRNAs
splitUp = dict(zip(biclusterHeader.strip().split(','),biclusters[j].split(',')))
# MEME1
if (not splitUp['MEME Motif1 E-Value']=='NA') and float(splitUp['MEME Motif1 E-Value'])<=0.05 and (not splitUp['Up.MEME Motif1 Matches']=='NA'): # and len(splitUp[10].split(' '))<10:
tfs[i]['MEME'] += splitUp['Up.MEME Motif1 Matches'].split(' ')
if not splitUp['Up.MEME Motif1 Correlated Matches']=='NA':
correlatedTfs[i]['MEME'] += [k.split(':')[0] for k in splitUp['Up.MEME Motif1 Correlated Matches'].split(' ')]
bestTfs[i]['MEME'] += [splitUp['Up.MEME Motif1 Minimum Correlated'].split(':')[0]]
# MEME2
if (not splitUp['Up.MEME Motif2 E-Value']=='NA') and float(splitUp['Up.MEME Motif2 E-Value'])<=0.05 and (not splitUp['Up.MEME Motif2 Matches']=='NA'): # and len(splitUp[14].split(' '))<10:
tfs[i]['MEME'] += splitUp['Up.MEME Motif2 Matches'].split(' ')
if not splitUp['Up.MEME Motif2 Correlated Matches']=='NA':
correlatedTfs[i]['MEME'] += [k.split(':')[0] for k in splitUp['Up.MEME Motif2 Correlated Matches'].split(' ')]
bestTfs[i]['MEME'] += [splitUp['Up.MEME Motif2 Minimum Correlated'].split(':')[0]]
# Weeder
if (not splitUp['Up.WEEDER Motif1 Matches']=='NA'): # and len(splitUp[17].split(' '))<10:
tfs[i]['WEEDER'] += splitUp['Up.WEEDER Motif1 Matches'].split(' ')
if not splitUp['Up.WEEDER Motif1 Correlated Matches']=='NA':
correlatedTfs[i]['WEEDER'] += [k.split(':')[0] for k in splitUp['Up.WEEDER Motif1 Correlated Matches'].split(' ')]
bestTfs[i]['WEEDER'] += [splitUp['Up.WEEDER Motif1 Minimum Correlated'].split(':')[0]]
# Weeder
if (not splitUp['Up.WEEDER Motif2 Matches']=='NA'): # and len(splitUp[20].split(' '))<10:
tfs[i]['WEEDER'] += splitUp['Up.WEEDER Motif2 Matches'].split(' ')
if not splitUp['Up.WEEDER Motif2 Correlated Matches']=='NA':
correlatedTfs[i]['WEEDER'] += [k.split(':')[0] for k in splitUp['Up.WEEDER Motif2 Correlated Matches'].split(' ')]
bestTfs[i]['WEEDER'] += [splitUp['Up.WEEDER Motif2 Minimum Correlated'].split(':')[0]]
# TFBS_DB
if (not splitUp['TFBS_DB.percTargets']=='NA') and float((splitUp['TFBS_DB.percTargets'].split(' '))[0])>=0.1 and float((splitUp['TFBS_DB.pValue'].split(' '))[0])<=0.05:
tfs[i]['TFBS_DB'] += splitUp['TFBS_DB.TFs'].split(' ')
if not splitUp['TFBS_DB.Correlated Matches']=='NA':
correlatedTfs[i]['TFBS_DB'] += [k.split(':')[0] for k in splitUp['TFBS_DB.Correlated Matches'].split(' ')]
bestTfs[i]['TFBS_DB'] += [splitUp['TFBS_DB.Minimum Correlated'].split(':')[0]]
# WEEDER:miRvestigator miRNA?
if splitUp['3pUTR.WEEDER Motif1 Model']=='8mer':
miRNAs[i] += [splitUp['3pUTR.WEEDER Motif1 Matches']]
if splitUp['3pUTR.WEEDER Motif2 Model']=='8mer':
miRNAs[i] += [splitUp['3pUTR.WEEDER Motif2 Matches']]
# PITA miRNA?
if (not splitUp['3pUTR_pita.percTargets']=='NA') and float((splitUp['3pUTR_pita.percTargets'].split(' '))[0])>=0.1 and float((splitUp['3pUTR_pita.pValue'].split(' '))[0])<=0.05: # and len(splitUp[36].split(' '))<10:
miRNAs[i] += splitUp['3pUTR_pita.miRNAs'].split(' ')
# TargetScan miRNA?
if (not splitUp['3pUTR_targetScan.percTargets']=='NA') and float((splitUp['3pUTR_targetScan.percTargets'].split(' '))[0])>=0.1 and float((splitUp['3pUTR_targetScan.pValue'].split(' '))[0])<=0.05: # and len(splitUp[39].split(' '))<10:
miRNAs[i] += splitUp['3pUTR_targetScan.miRNAs'].split(' ')
# GO Biological Processes
if not splitUp['GO_Term_BP']=='NA':
goBP[i] += splitUp['GO_Term_BP'].split(';')
outFile = open('postProcessed_clustersOfBiclusters.csv','w')
outFile.write('\n'.join(writeMe))
outFile.close()
# Dump out TFs to test
outFile = open('tf_regulators.csv','w')
outFile.write('cluster,MEME,WEEDER,TFBS_DB,Number Predicted\n')
outFile.write('\n'.join([i+','+';'.join(list(set(tfs[i]['MEME'])))+','+';'.join(list(set(tfs[i]['WEEDER'])))+','+';'.join(list(set(tfs[i]['TFBS_DB'])))+','+str(len(list(set(list(set(tfs[i]['MEME']))+list(set(tfs[i]['WEEDER']))+list(set(tfs[i]['TFBS_DB'])))))) for i in tfs]))
outFile.close()
# Dump out TFs to test
outFile = open('correlated_tf_regulators.csv','w')
outFile.write('cluster,MEME,WEEDER,TFBS_DB,Number Predicted\n')
outFile.write('\n'.join([i+','+';'.join(list(set(correlatedTfs[i]['MEME'])))+','+';'.join(list(set(correlatedTfs[i]['WEEDER'])))+','+';'.join(list(set(correlatedTfs[i]['TFBS_DB'])))+','+str(len(list(set(list(set(correlatedTfs[i]['MEME']))+list(set(correlatedTfs[i]['WEEDER']))+list(set(correlatedTfs[i]['TFBS_DB'])))))) for i in correlatedTfs]))
outFile.close()
# Dump out TFs to test
outFile = open('best_tf_regulators.csv','w')
outFile.write('cluster,MEME,WEEDER,TFBS_DB,Number Predicted\n')
outFile.write('\n'.join([i+','+';'.join(list(set(bestTfs[i]['MEME'])))+','+';'.join(list(set(bestTfs[i]['WEEDER'])))+','+';'.join(list(set(bestTfs[i]['TFBS_DB'])))+','+str(len(list(set(list(set(bestTfs[i]['MEME']))+list(set(bestTfs[i]['WEEDER']))+list(set(bestTfs[i]['TFBS_DB'])))))) for i in bestTfs]))
outFile.close()
# Dump out miRNAs to test
outFile = open('miRNA_regulators.csv','w')
outFile.write('\n'.join([i+','+';'.join(list(set(miRNAs[i])))+','+str(len(list(set(miRNAs[i])))) for i in miRNAs]))
outFile.close()
# Dump out GO BPs to test
outFile = open('GO_BPs_metaClusters.csv','w')
outFile.write('\n'.join([i+','+';'.join(list(set(goBP[i])))+','+str(len(list(set(goBP[i])))) for i in goBP]))
outFile.close()
# Read in clusters of biclusters
clusters = {}
inFile = open('clustersOfBiclusters_VALL.csv','r')
inFile.readline() # Get rid of header
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.strip().split(',')
if not splitUp[1] in clusters:
clusters[splitUp[1]] = []
clusters[splitUp[1]].append(splitUp[0].strip('"').lower())
inFile.close()
# Read in post processed files
#Bicluster - 0 Genes - 1 Conditions - 2 Norm. Residual - 3 Norm. Residual Perm. P-Value - 4 Var. Exp. First PC - 5 Var. Exp. First PC Perm. P-Value - 6 MEME Motif1 E-Value - 7 Up.MEME Motif1 Consensus - 8 Up.MEME Motif1 Matches - 9 Up.MEME Motif1 Expanded Matches - 10 Up.MEME Motif1 Correlated Matches - 11 Up.MEME Motif1 Original/Expanded - 12 Up.MEME Motif1 Minimum Correlated - 13 Up.MEME Motif2 E-Value - 14 Up.MEME Motif2 Consensus - 15 Up.MEME Motif2 Matches - 16 Up.MEME Motif2 Expanded Matches - 17 Up.MEME Motif2 Correlated Matches - 18 Up.MEME Motif2 Original/Expanded - 19 Up.MEME Motif2 Minimum Correlated - 20 Up.WEEDER Motif1 Score - 21 Up.WEEDER Motif1 Consensus - 22 Up.WEEDER Motif1 Matches - 23 Up.WEEDER Motif1 Expanded Matches - 24 Up.WEEDER Motif1 Correlated Matches - 25 Up.WEEDER Motif1 Original/Expanded - 26 Up.WEEDER Motif1 Minimum Correlated - 27 Up.WEEDER Motif2 Score - 28 Up.WEEDER Motif2 Consensus - 29 Up.WEEDER Motif2 Matches - 30 Up.WEEDR Motif2 Expanded Matches - 31 Up.WEEDER Motif2 Correlated Matches - 32 Up.WEEDER Motif2 Original/Expanded - 33 Up.WEEDER Motif2 Minimum Correlated - 34 TFBS_DB.TFs - 35 TFBS_DB.percTargets - 36 TFBS_DB.pValue - 37 TFBS_DB.Exapnded Matches - 38 TFBS_DB.Correlated Matches - 39 TFBS_DB.Original/Expanded - 40 TFBS_DB.Minimum Correlated - 41 3pUTR.WEEDER Motif1 E-Value - 42 3pUTR.WEEDER Motif1 Consensus - 43 3pUTR.WEEDER Motif1 Matches - 44 3pUTR.WEEDER Motif1 Model - 45 3pUTR.WEEDER Motif2 E-Value - 46 3pUTR.WEEDER Motif2 Consensus - 47 3pUTR.WEEDER Motif2 Matches - 48 3pUTR.WEEDER Motif2 Model - 49 3pUTR_pita.miRNAs - 50 3pUTR_pita.percTargets - 51 3pUTR_pita.pValue - 52 3pUTR_targetScan.miRNAs - 53 3pUTR_targetScan.percTargets - 54 3pUTR_targetScan.pValue - 55 GO_Term_BP - 56
biclusters = {}
biclusterHeader = ''
for i in ['pita','targetscan','tfbs_db']:
inFile = open('postProcessed/postProcessed_urdahl_alvMac_'+i+'.csv','r')
biclusterHeader = inFile.readline() # Get rid of header
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.split(',')
biclusters[i+'_'+splitUp[0]] = line.strip()
inFile.close()
inFile = open('genesInBiclusters/cluster.members.genes_'+i+'.txt','r')
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.strip().split(' ')
j = splitUp.pop(0)
biclusters[i+'_'+j] = biclusters[i+'_'+j] + ',' + ';'.join(splitUp)
inFile.close()
inFile = open('genesInBiclusters/cluster.members.conditions_'+i+'.txt','r')
while 1:
line = inFile.readline()
if not line:
break
splitUp = line.strip().split(' ')
j = splitUp.pop(0)
biclusters[i+'_'+j] = biclusters[i+'_'+j] + ',' + ';'.join(splitUp)
inFile.close()
# Dump combined postProcessed file with clusters IDs
#Bicluster - 0 Genes - 1 Conditions - 2 Norm. Residual - 3 Norm. Residual Perm. P-Value - 4 Var. Exp. First PC - 5 Var. Exp. First PC Perm. P-Value - 6 MEME Motif1 E-Value - 7 Up.MEME Motif1 Consensus - 8 Up.MEME Motif1 Matches - 9 Up.MEME Motif1 Expanded Matches - 10 Up.MEME Motif1 Correlated Matches - 11 Up.MEME Motif1 Original/Expanded - 12 Up.MEME Motif1 Minimum Correlated - 13 Up.MEME Motif2 E-Value - 14 Up.MEME Motif2 Consensus - 15 Up.MEME Motif2 Matches - 16 Up.MEME Motif2 Expanded Matches - 17 Up.MEME Motif2 Correlated Matches - 18 Up.MEME Motif2 Original/Expanded - 19 Up.MEME Motif2 Minimum Correlated - 20 Up.WEEDER Motif1 Score - 21 Up.WEEDER Motif1 Consensus - 22 Up.WEEDER Motif1 Matches - 23 Up.WEEDER Motif1 Expanded Matches - 24 Up.WEEDER Motif1 Correlated Matches - 25 Up.WEEDER Motif1 Original/Expanded - 26 Up.WEEDER Motif1 Minimum Correlated - 27 Up.WEEDER Motif2 Score - 28 Up.WEEDER Motif2 Consensus - 29 Up.WEEDER Motif2 Matches - 30 Up.WEEDR Motif2 Expanded Matches - 31 Up.WEEDER Motif2 Correlated Matches - 32 Up.WEEDER Motif2 Original/Expanded - 33 Up.WEEDER Motif2 Minimum Correlated - 34 TFBS_DB.TFs - 35 TFBS_DB.percTargets - 36 TFBS_DB.pValue - 37 TFBS_DB.Exapnded Matches - 38 TFBS_DB.Correlated Matches - 39 TFBS_DB.Original/Expanded - 40 TFBS_DB.Minimum Correlated - 41 3pUTR.WEEDER Motif1 E-Value - 42 3pUTR.WEEDER Motif1 Consensus - 43 3pUTR.WEEDER Motif1 Matches - 44 3pUTR.WEEDER Motif1 Model - 45 3pUTR.WEEDER Motif2 E-Value - 46 3pUTR.WEEDER Motif2 Consensus - 47 3pUTR.WEEDER Motif2 Matches - 48 3pUTR.WEEDER Motif2 Model - 49 3pUTR_pita.miRNAs - 50 3pUTR_pita.percTargets - 51 3pUTR_pita.pValue - 52 3pUTR_targetScan.miRNAs - 53 3pUTR_targetScan.percTargets - 54 3pUTR_targetScan.pValue - 55 GO_Term_BP - 56
tfs = {}
correlatedTfs = {}
bestTfs = {}
miRNAs = {}
goBP = {}
writeMe = ['Meta-Cluster,bicluster,'+biclusterHeader.strip()+',Gene Names,Condition Names']
for i in clusters:
tfs[i] = {'MEME':[], 'WEEDER':[], 'TFBS_DB':[]}
correlatedTfs[i] = {'MEME':[], 'WEEDER':[], 'TFBS_DB':[]}
bestTfs[i] = {'MEME':[], 'WEEDER':[], 'TFBS_DB':[]}
miRNAs[i] = []
goBP[i] = []
for j in clusters[i]:
writeMe.append(i+','+j+','+biclusters[j])
# Process out tfs and miRNAs
splitUp = dict(zip(biclusterHeader.strip().split(','),biclusters[j].split(',')))
# MEME1
if (not splitUp['MEME Motif1 E-Value']=='NA') and float(splitUp['MEME Motif1 E-Value'])<=0.05 and (not splitUp['Up.MEME Motif1 Matches']=='NA'): # and len(splitUp[10].split(' '))<10:
tfs[i]['MEME'] += splitUp['Up.MEME Motif1 Matches'].split(' ')
if not splitUp['Up.MEME Motif1 Correlated Matches']=='NA':
correlatedTfs[i]['MEME'] += [k.split(':')[0] for k in splitUp['Up.MEME Motif1 Correlated Matches'].split(' ')]
bestTfs[i]['MEME'] += [splitUp['Up.MEME Motif1 Minimum Correlated'].split(':')[0]]
# MEME2
if (not splitUp['Up.MEME Motif2 E-Value']=='NA') and float(splitUp['Up.MEME Motif2 E-Value'])<=0.05 and (not splitUp['Up.MEME Motif2 Matches']=='NA'): # and len(splitUp[14].split(' '))<10:
tfs[i]['MEME'] += splitUp['Up.MEME Motif2 Matches'].split(' ')
if not splitUp['Up.MEME Motif2 Correlated Matches']=='NA':
correlatedTfs[i]['MEME'] += [k.split(':')[0] for k in splitUp['Up.MEME Motif2 Correlated Matches'].split(' ')]
bestTfs[i]['MEME'] += [splitUp['Up.MEME Motif2 Minimum Correlated'].split(':')[0]]
# Weeder
if (not splitUp['Up.WEEDER Motif1 Matches']=='NA'): # and len(splitUp[17].split(' '))<10:
tfs[i]['WEEDER'] += splitUp['Up.WEEDER Motif1 Matches'].split(' ')
if not splitUp['Up.WEEDER Motif1 Correlated Matches']=='NA':
correlatedTfs[i]['WEEDER'] += [k.split(':')[0] for k in splitUp['Up.WEEDER Motif1 Correlated Matches'].split(' ')]
bestTfs[i]['WEEDER'] += [splitUp['Up.WEEDER Motif1 Minimum Correlated'].split(':')[0]]
# Weeder
if (not splitUp['Up.WEEDER Motif2 Matches']=='NA'): # and len(splitUp[20].split(' '))<10:
tfs[i]['WEEDER'] += splitUp['Up.WEEDER Motif2 Matches'].split(' ')
if not splitUp['Up.WEEDER Motif2 Correlated Matches']=='NA':
correlatedTfs[i]['WEEDER'] += [k.split(':')[0] for k in splitUp['Up.WEEDER Motif2 Correlated Matches'].split(' ')]
bestTfs[i]['WEEDER'] += [splitUp['Up.WEEDER Motif2 Minimum Correlated'].split(':')[0]]
# TFBS_DB
if (not splitUp['TFBS_DB.percTargets']=='NA') and float((splitUp['TFBS_DB.percTargets'].split(' '))[0])>=0.1 and float((splitUp['TFBS_DB.pValue'].split(' '))[0])<=0.05:
tfs[i]['TFBS_DB'] += splitUp['TFBS_DB.TFs'].split(' ')
if not splitUp['TFBS_DB.Correlated Matches']=='NA':
correlatedTfs[i]['TFBS_DB'] += [k.split(':')[0] for k in splitUp['TFBS_DB.Correlated Matches'].split(' ')]
bestTfs[i]['TFBS_DB'] += [splitUp['TFBS_DB.Minimum Correlated'].split(':')[0]]
# WEEDER:miRvestigator miRNA?
if splitUp['3pUTR.WEEDER Motif1 Model']=='8mer':
miRNAs[i] += [splitUp['3pUTR.WEEDER Motif1 Matches']]
if splitUp['3pUTR.WEEDER Motif2 Model']=='8mer':
miRNAs[i] += [splitUp['3pUTR.WEEDER Motif2 Matches']]
# PITA miRNA?
if (not splitUp['3pUTR_pita.percTargets']=='NA') and float((splitUp['3pUTR_pita.percTargets'].split(' '))[0])>=0.1 and float((splitUp['3pUTR_pita.pValue'].split(' '))[0])<=0.05: # and len(splitUp[36].split(' '))<10:
miRNAs[i] += splitUp['3pUTR_pita.miRNAs'].split(' ')
# TargetScan miRNA?
if (not splitUp['3pUTR_targetScan.percTargets']=='NA') and float((splitUp['3pUTR_targetScan.percTargets'].split(' '))[0])>=0.1 and float((splitUp['3pUTR_targetScan.pValue'].split(' '))[0])<=0.05: # and len(splitUp[39].split(' '))<10:
miRNAs[i] += splitUp['3pUTR_targetScan.miRNAs'].split(' ')
# GO Biological Processes
if not splitUp['GO_Term_BP']=='NA':
goBP[i] += splitUp['GO_Term_BP'].split(';')
outFile = open('postProcessed_clustersOfBiclusters.csv','w')
outFile.write('\n'.join(writeMe))
outFile.close()
# Dump out TFs to test
outFile = open('tf_regulators.csv','w')
outFile.write('cluster,MEME,WEEDER,TFBS_DB,Number Predicted\n')
outFile.write('\n'.join([i+','+';'.join(list(set(tfs[i]['MEME'])))+','+';'.join(list(set(tfs[i]['WEEDER'])))+','+';'.join(list(set(tfs[i]['TFBS_DB'])))+','+str(len(list(set(list(set(tfs[i]['MEME']))+list(set(tfs[i]['WEEDER']))+list(set(tfs[i]['TFBS_DB'])))))) for i in tfs]))
outFile.close()
# Dump out TFs to test
outFile = open('correlated_tf_regulators.csv','w')
outFile.write('cluster,MEME,WEEDER,TFBS_DB,Number Predicted\n')
outFile.write('\n'.join([i+','+';'.join(list(set(correlatedTfs[i]['MEME'])))+','+';'.join(list(set(correlatedTfs[i]['WEEDER'])))+','+';'.join(list(set(correlatedTfs[i]['TFBS_DB'])))+','+str(len(list(set(list(set(correlatedTfs[i]['MEME']))+list(set(correlatedTfs[i]['WEEDER']))+list(set(correlatedTfs[i]['TFBS_DB'])))))) for i in correlatedTfs]))
outFile.close()
# Dump out TFs to test
outFile = open('best_tf_regulators.csv','w')
outFile.write('cluster,MEME,WEEDER,TFBS_DB,Number Predicted\n')
outFile.write('\n'.join([i+','+';'.join(list(set(bestTfs[i]['MEME'])))+','+';'.join(list(set(bestTfs[i]['WEEDER'])))+','+';'.join(list(set(bestTfs[i]['TFBS_DB'])))+','+str(len(list(set(list(set(bestTfs[i]['MEME']))+list(set(bestTfs[i]['WEEDER']))+list(set(bestTfs[i]['TFBS_DB'])))))) for i in bestTfs]))
outFile.close()
# Dump out miRNAs to test
outFile = open('miRNA_regulators.csv','w')
outFile.write('\n'.join([i+','+';'.join(list(set(miRNAs[i])))+','+str(len(list(set(miRNAs[i])))) for i in miRNAs]))
outFile.close()
# Dump out GO BPs to test
outFile = open('GO_BPs_metaClusters.csv','w')
outFile.write('\n'.join([i+','+';'.join(list(set(goBP[i])))+','+str(len(list(set(goBP[i])))) for i in goBP]))
outFile.close()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment