###

source('./utilities.R')

num_topic=100

version='rna_atac'

#read gene by topic matrix, and all gene names

rna_topic=read.csv(paste0('../Result/',version,'/topic_rna.csv'),header=T)

rna_names=rna_topic[,1]

rna_topic=rna_topic[,-1]

rownames(rna_topic)=rna_names

colnames(rna_topic)=paste0('topic',1:num_topic)

#read protein by topic matrix, and all protein names

atac_topic=read.csv(paste0('../Result/',version,'/topic_atac.csv'),header=T)

peak_names=atac_topic[,1]

atac_topic=atac_topic[,-1]

rownames(atac_topic)=peak_names

colnames(atac_topic)=paste0('topic',1:num_topic)

##read all gene locations

gene2loci=read.csv('../useful_file/rna_atac/gene2loci2.csv',header=T)

##split peak name

peak_split=c()

for (i in peak_names){peak_split=rbind(peak_split,strsplit(i,'-')[[1]])}

peak_split=data.frame(peak_split)

colnames(peak_split)=c('chr_name','start','end')

peak_split$start=as.numeric(peak_split$start)

peak_split$end=as.numeric(peak_split$end)

##

library('GenomicRanges')

peak_gr=GRanges(seqnames = peak_split$chr_name,

                ranges = IRanges(peak_split$start, end = peak_split$end))

rna_gr=GRanges(seqnames = gene2loci$chr_name,

               ranges = IRanges(gene2loci$start_position, 

                                end = gene2loci$end_position,

                                names = gene2loci$gene_name))

simu_topic_rna=c()

for (i in 1:100){#take about 1hr

  topic=c()#converted topic_i x gene

  peak_topic_v=atac_topic[,i]#1 x num_peak

  for (j in rna_names){#gene names

    d=distance(rna_gr[j],peak_gr)

    d=ifelse(d<=150000,1,0)#distance constraint

    ind=which(d==1)#index for peak

    for (k in ind){#check non-0 peak correlation

      corr=cor(as.numeric(rna_topic[j,]),as.numeric(atac_topic[k,]))

      if (corr<=0){d[k]=0}#positive constraint

    }

    d[is.na(d)]=0

    topic=c(topic,peak_topic_v%*%d) 

  }

  simu_topic_rna=rbind(simu_topic_rna,topic)

}

write.csv(simu_topic_rna,file=paste0('../data/',version,'/simu_topic_rna.csv'),row.names = F,quote = F)

##

#imu_topic_rna=read.csv('../data/rna_atac/v2/simu_topic_rna.csv',header=T)

simu_rna_topic=t(simu_topic_rna)

rownames(simu_rna_topic)=rna_names

colnames(simu_rna_topic)=paste0('topic',1:num_topic)

### analyze rna_topic and simu_rna_topic

##correlation

plot_corr=c()

for (i in 1:num_topic){

  plot_corr=c(plot_corr,cor(rna_topic[,i],simu_rna_topic[,i]))

}

save_path=paste0('../plot/',version)

write.csv(plot_corr,file=paste0(save_path,'/plot_correlation.csv'),quote=F)

corr_color=rep(NA,100)

corr_color[plot_corr>=0]='blue'

corr_color[plot_corr<0]='red'

corr_plot=data.frame(x=1:100,y=plot_corr,color=corr_color)

save_name=paste0(save_path,'/gene_atac_corr.png')

plot_correlation(corr_plot,save_name)

## top features

gene_ids=get_feature_orders(rna_topic,num_topic=100)

simu_gene_ids=get_feature_orders(simu_rna_topic,num_topic=100)

topic_select=c(3,32,30,83,21,79,50)

top_feature_num=5

save_path=paste0('../plot/',version,'/')

plot_top_feature_in_selected_topic('atac',top_feature_num,simu_gene_ids,simu_rna_topic,topic_select,save_path,'select')

plot_top_feature_in_selected_topic('gene',top_feature_num,gene_ids,rna_topic,topic_select,save_path,'select')

## save top5 features as csv

top_rna=c()

top_atac=c()

for (i in topic_select){

  top_rna=cbind(top_rna,gene_ids[1:5,i])

  top_atac=cbind(top_atac,simu_gene_ids[1:5,i])

}

##sample x topic heatmap

version='rna_atac'

is_plot=T;is_diff_topic=F

type='rna_atac'

topic_select=','

is_all=F

args=c(num_topic,version,type,topic_select,is_plot,is_diff_topic,is_all)

system(paste(c('Rscript','./plot_sample_by_topic.R',args),collapse=' '))

#

is_plot=F;is_diff_topic=T

topic_select='30,78,3,72,84,51,55,92,50,73,39,46,33,59,12,21,79,16,98,32,25,26,38,83,88'

args=c(num_topic,version,type,topic_select,is_plot,is_diff_topic,is_all)

system(paste(c('Rscript','./plot_sample_by_topic.R',args),collapse=' '))

## prepare GSEA

############# prepare GSEA files

#rna_topic

#simu_rna_topic=t(simu_topic_rna)

save_path_g=paste0('../data/',version,'/rna/')

save_path_p=paste0('../data/',version,'/atac/')

for (i in 1:num_topic){

  value=cbind(rownames(rna_topic),rna_topic[,i])

  write.table(value,file=paste0(save_path_g,'geneName_topic',i,'.rnk'),quote=F,sep='\t',row.names = F,col.names=F)

  value=cbind(rownames(simu_rna_topic),simu_rna_topic[,i])

  write.table(value,file=paste0(save_path_p,'geneName_topic',i,'.rnk'),quote=F,sep='\t',row.names = F,col.names=F)

}

##### prepare MEME input

## top 1000 peaks

library(httr)

library(jsonlite)

library(xml2)

server <- "http://rest.ensembl.org"

###

topic_select=c(3,32,30,83,21,79,50)

#t=3

for (to in topic_select){

  peak=atac_topic[,to]

  top_peaks=peak_names[order(peak,decreasing=T)[1:100]]

  top_peaks_input=c()

  for (i in top_peaks){

    s=strsplit(i,'-')[[1]]#split

    t=paste0(s[1],':',s[2],'..',s[3],':1?')

    top_peaks_input=c(top_peaks_input,t)

  }

  #get sequences

  seqs=c()

  for (i in top_peaks_input){

    ext=paste0("/sequence/region/human/",i)

    r <- GET(paste(server, ext, sep = ""), content_type("text/x-fasta"))

    stop_for_status(r)

    seqs=c(seqs,content(r))

  }

  out_file=file(paste0("../data/rna_atac/seq/seqs100_topic",to,".fasta"))

  writeLines(seqs,out_file)

  close(out_file)

}