R/get_HiCUP.R

In Skhakmardan/MHiC: MHiC

get_HiCUP<-function(reads_file,Digest_file){
  
  compressed<-grepl("\\.gz$", reads_file)
  filetype<-grepl("\\.table$", reads_file)
  if(!filetype){
  filetype<-grepl("\\.txt$", reads_file)}
  ##########################################################
  
  #import digest file
  ranges<-hicup_ranges(Digest_file)
  
  ##########################################################
  
  #import  reads file
  if(compressed){
    
    reads<-read.table(gzfile(reads_file)) 
    
  } 
  else{
  if(filetype){
    reads<-read.table(reads_file) 
  }
 
  else{stop('HiCUP output must be in a table format')}
  
  }
  names(reads) <- c("id", "flag", "chr", "locus")
  reads<-reads[order(reads$id, reads$flag, reads$chr, reads$locus),]
  ##########################################################
  
  
  #in the sam file the two ends of the interactions are consecutive rows
  odd <- reads[seq(1,nrow(reads),2), ]
  even <- reads[seq(2,nrow(reads),2), ]
  names(odd) <- c("id1", "flag1", "chr1", "locus1")
  names(even) <- c("id2", "flag2", "chr2", "locus2")	
  
  levels(odd$chr1)<-fixChromosomeNames(levels(odd$chr1))
  levels(even$chr2)<-fixChromosomeNames(levels(even$chr2))
  
  
  if(nrow(odd) != nrow(even)){stop("odd and even rows of data must have same counts")}
  joined <- cbind(odd, even)
  if(!all(joined$id1==joined$id2)){stop(" each odd and even row of data must have same id ")}
  ##########################################################
  
  
  interactions<-joined[,c("chr1", "locus1","chr2", "locus2")]
  interactingLoci<-map_reads_HiCUP(interactions,ranges)
  ##########################################################
  
  return(interactingLoci)
  
}