R/get_Hic_pro.R

In Skhakmardan/MHiC: MHiC

get_Hic_pro<-function(reads_file,Digest_file = NULL){

  filetype<-grepl("\\.matrix$", reads_file)
  if(!filetype){filetype<-grepl("\\.table$", reads_file)}
  else if(!filetype){
  filetype<-grepl("\\.txt$", reads_file)}
  ##########################################################
  #import digest file

  ranges<-read.table(Digest_file, header = FALSE)

  ranges$chr<-ranges$V1
  ranges$start<-ranges$V2
  ranges$end<-ranges$V3
  ranges$id<-ranges$V4
  ranges<-ranges[,c("chr","start","end","id")]
  levels(ranges$chr) <- fixChromosomeNames(levels(ranges$chr))
  ##########################################################
  #import  reads file
  if(filetype){
    reads<-read.table(reads_file)
  }
  else{stop('Hic_pro output must be in a matrix format')}

  colnames(reads)<-c("locus1","locus2","interaction_counts")
  ##########################################################
  # map reads
  fragments1<- reads$locus1
  fragments2<- reads$locus2

  loci1 <- ranges[fragments1,]
  loci2 <- ranges[fragments2,]
  rm(fragments1,fragments2)
  ##########################################################
  interactions<-data.frame(loci1$chr,loci1$start,loci2$chr,loci2$start
                           ,reads$interaction_counts)

  names(interactions)<-c("chr1", "locus1_start", "chr2", "locus2_start"
                         ,"frequencies")
  return(interactions)

}