R/pcfPlain.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
####################################################################


##Required by:


##Requires:
### findNN
### getMad
### pcf(not the main function, but the rest of the help functions found in the same document)
### handleMissing
### pullOutContent


## Main function for pcf-analysis to be called by the user

pcfPlain <- function(pos.data,kmin=5,gamma=40,normalize=TRUE,fast=TRUE,digits=4,return.est=FALSE,verbose=TRUE){
  
  
  #Input could come from winsorize and thus be a list; check and possibly retrieve data frame wins.data
  pos.data <- pullOutContent(pos.data,what="wins.data")
    
  #Make sure all data columns are numeric:
  if(any(!sapply(pos.data,is.numeric))){
    stop("All input in pos.data must be numeric",call.=FALSE)
  }
  
  #Check data input:
  stopifnot(ncol(pos.data)>=2)  #something is missing in input data
  
  #Extract information from pos.data:
  position <- pos.data[,1]
  cn.data <- as.matrix(pos.data[,-1])
  nSample <- ncol(pos.data)-1
  sampleid <- colnames(pos.data)[-1]  
  nProbe <- length(position)
	
  
  #save user's gamma
	gamma0 <- gamma     
	
	sd <- rep(1,nSample) #sd is only used if normalize=TRUE, and then these values are replaced by MAD-sd
	#If number of probes in entire data set is less than 100K, the MAD sd-estimate is calculated using all obs for every sample
  #Only required if normalize=T
  if(nProbe<100000 && normalize){
    #calculate MAD-sd for each sample:
    for(j in 1:nSample){
      sample.data <- pos.data[,j+1]
      sd[j] <- getMad(sample.data[!is.na(sample.data)],k=25)   #Take out missing values before calculating mad
		}
  }#endif
  
  
	#Initialize
	pcf.names <- c("pos",sampleid)
	seg.names <- c("sampleID","start.pos","end.pos","n.probes","mean")
	
  segments <- data.frame(matrix(nrow=0,ncol=5))
  colnames(segments) <- seg.names
  if(return.est){
    pcf.est <- matrix(nrow=0,ncol=nSample)
	}
  
  		
  #Run PCF separately for each sample:
  for(i in 1:nSample){
    if(return.est){
      #Initialize:
      yhat <- rep(NA,length(nProbe))
    }
    
    sample.data <- cn.data[,i]
    
		#Remove probes with missing obs; Only run pcf on non-missing values
		obs <- !is.na(sample.data)
		obs.data <- sample.data[obs]
	
	  
	  if(length(obs.data)>0){  ##Make sure there are observations for this sample! If not, estimates are left NA as well
		
			#If number of probes in entire data set is >= 100K, the MAD sd-estimate is calculated using obs in this arm for this sample.
      #Only required if normalize=T 
      if(nProbe>=100000 && normalize){
        sd[i] <- getMad(obs.data,k=25)   
      }
    
      #Scale gamma by variance if normalize is TRUE 
      use.gamma <- gamma0
      if(normalize){
        use.gamma <- gamma0*(sd[i])^2
      }
      
  		#Must check that sd!=0 and sd!!=NA -> use.gamma=0/NA. If not, simply calculate mean of observations
			if(use.gamma==0 || is.na(use.gamma)){  
			  if(return.est){
          res <- list(Lengde=length(obs.data),sta=1,mean=mean(obs.data),nIntervals=1,yhat=rep(mean(obs.data)))
        }else{
          res <- list(Lengde=length(obs.data),sta=1,mean=mean(obs.data),nIntervals=1)
        }
        
			}else{
  			# Compute piecewise constant fit
  			#run fast approximate PCF if fast=TRUE and number of probes>400, or exact PCF otherwise
  			if(!fast || length(obs.data)<400){
  				#Exact PCF:	
  			 	res <- exactPcf(y=obs.data,kmin=kmin,gamma=use.gamma,yest=return.est)	 			
  			}else{
          #Run fast PCF:
          res <- selectFastPcf(x=obs.data,kmin=kmin,gamma=use.gamma,yest=return.est)
  			}#endif
  			
  			
			}#endif
			
			#Retrieve segment info from results:
			seg.start <- res$sta
		  seg.stop <- c(seg.start[-1]-1,length(obs.data))
		  seg.npos <- res$Lengde
			seg.mean <- res$mean
			nSeg <- res$nIntervals 
      if(return.est){
        yhat[obs] <- res$yhat
      } 
      #Find genomic positions for start and stop of each segment:
		  pos.start <- position[seg.start]
		  pos.stop <- position[seg.stop]
			
      #Handle missing values:
			if(any(!obs)){
			  #first find nearest non-missing neighbour for missing probes:
        nn <- findNN(pos=position,obs=obs)
     
        #Include probes with missing values in segments where their nearest neighbour probes are located
        new.res <- handleMissing(nn=nn,pos=position,obs=obs,pos.start=pos.start,pos.stop=pos.stop,seg.npos=seg.npos)  
        pos.start <- new.res$pos.start
        pos.stop <- new.res$pos.stop
        seg.npos <- new.res$seg.npos
       
        if(return.est){
          yhat[!obs] <- yhat[nn]
        }
			}
			
		}else{
		  warning(paste("pcf is not run for sample ",i," because all observations are missing. NA is returned.",sep=""),immediate.=TRUE,call.=FALSE)
		  seg.start <- 1
		  seg.stop <- nProbe
	    pos.start <- position[seg.start]
	    pos.stop <- position[seg.stop]
	    nSeg <- 1
	    seg.mean <- NA
	    seg.npos <- nProbe
	  }  
		
		
		#Round:
		seg.mean <- round(seg.mean,digits=digits)
		
    #Add results for this sample to results for other samples in data frame:
    seg <- data.frame(rep(sampleid[i],nSeg),pos.start,pos.stop,seg.npos,seg.mean,stringsAsFactors=FALSE)
		colnames(seg) <- seg.names
		segments <- rbind(segments,seg)

    if(return.est){
      #Rounding:
		  yhat <- round(yhat,digits=digits)
      pcf.est <- cbind(pcf.est,yhat)
    }
  }#endfor
  
  
    	
  if(verbose){
    cat(paste("pcf finished for sample ",i,sep=""),"\n")
  }

  
	if(return.est){
    pcf.est <- data.frame(position,pcf.est,stringsAsFactors=FALSE)
	  colnames(pcf.est) <- pcf.names
		return(list(estimates=pcf.est,segments=segments))
	}else{
	 return(segments)
	}
	
}#endfunction