R/subsetSegments.r

In igordot/copynumber: Segmentation of single- and multi-track copy number data by penalized least squares regression (with hg38 and mm10).

####################################################################
## Author: Gro Nilsen, Knut Liestøl and Ole Christian Lingjærde.
## Maintainer: Gro Nilsen <gronilse@ifi.uio.no>
## License: Artistic 2.0
## Part of the copynumber package
## Reference: Nilsen and Liestøl et al. (2012), BMC Genomics
####################################################################

#Get a subset of segments for given chromosomes and/or sampleIDs

##Input:
### segments : either a data frame or the name of a file containing the segments found by either pcf, multipcf or aspcf
### chrom : a numeric vector with chromosome number(s) for which data or segments should be selected. If unspecified, all chromosomes in data or segments will be selected
### sample : a numeric vector indicating which sample(s) is to selected. The number(s) should correspond to the sample's place (in order of appearance) in the data 
### sep : the separator of the input files if \code{data} is a file. Default is "\t"
### ... : optional parameters to be passed to \code{read.table} in the case where data are to be read from files

## Output:
### sel.segments : the data frame with segments selected on chromosomes and sampleIds

##Required by:
### none

##Requires:
### is.multiseg
### numericChrom
### pullOutContent


subsetSegments <- function(segments,chrom=NULL,sample=NULL,sep="\t",...){

  #Check if segments is a file:
  isfile <- class(segments)=="character"
  
  #get header and chrom from data
  if(isfile){
    #read segment-file
    segments <- read.table(segments,header=TRUE,sep=sep,as.is=TRUE)
  }else{
    #Make sure segments is a data frame
    segments <- pullOutContent(res=segments,what="segments")
  }
  
  #Check whether we have multiPCF segments:
  multi <- is.multiseg(segments)
  
  if(multi){
    seg.sampleid <- colnames(segments)[-c(1:5)]
    seg.chrom <- segments[,1]
  }else{
    seg.sampleid <- segments[,1]
    seg.chrom <- segments[,2]
  }
  
  #make sure chrom in segments are numeric
  seg.chrom <- numericChrom(seg.chrom)
  
  #Pick out relevant chromosomes:
	keepchrom <- 1:length(seg.chrom)
	if(!is.null(chrom)){
    #make sure selected chrom are numeric:
    chrom <- numericChrom(chrom)
    #Check that these chrom are found in data.chrom
    use <- chrom%in%unique(seg.chrom)
    use.chrom <- chrom[use]
    if(length(use.chrom)==0){
      msg <- "None of the specified chromosomes are found in segments"
      stop(msg,call.=FALSE)
    }else if(length(use.chrom)!=length(chrom)){
      not.use <- paste(chrom[!use],sep="",collapse=",")
      msg <- paste("The following chromosome(s) are not found in segments:",not.use,sep=" ")
      warning(msg,call.=FALSE,immediate.=TRUE)
    }
    keepchrom <- which(seg.chrom%in%use.chrom)
	}
	
	#Select segments for the desired chromosome number:
	sel.segments <- segments[keepchrom,,drop=FALSE]

  #Check that specified sampleIDs are found in segments:
  if(!is.null(sample)){
    #Get sampleID for selected samples:
    id <- as.character(unique(seg.sampleid))
    sampleID <- id[sample]
    if(any(is.na(sampleID))){
      stop("Input in 'sample' is outside the number of samples represented in segments",.call=FALSE)
    }

    #selct segments for this sampleID(s)
    if(multi){
      keepsample <- which(colnames(segments)%in%sampleID)
		  sel.segments <- sel.segments[,c(1:5,keepsample)]
      colnames(sel.segments) <- c(colnames(segments)[1:5],colnames(segments)[keepsample])  
    }else{
      keepsample <- which(sel.segments[,1]%in%sampleID)
      sel.segments <- sel.segments[keepsample,,drop=FALSE]
    }
  }
  
  return(sel.segments)

}