getting a matrix methylation status matrix via QuasR

getCMethMatrix.R

  require(data.table)
  require(QuasR)
 
# arguments:  
# proj: qProject object
# range: GRanges object with ONE!!!! range
# samp: sample.name
#
getCMethMatrix<-function(proj,range,samp){  
  Cs=qMeth(proj, query=range,mode="allC",reportLevel="alignment")
  
  # use data.table to get a 1,0 matrix of methylation profiles
  all.cids=unique(Cs[[samp]]$Cid) # get all possible C locations
  
  # make the data.table object
  dt=data.table(meth=Cs[[samp]]$meth ,aid=Cs[[samp]]$aid ,cid=Cs[[samp]]$Cid) 
  
  # this function converts cids to columns
  myfun2<-function(x,all.cids){
    vec=rep(-1,length(all.cids))
    names(vec)=as.character(all.cids)
    b=as.list((vec))        
    b[ as.character(x$cid)]=as.double(x$meth)
    return(b)
  }
  dtm=dt[,myfun2(.SD,all.cids), by=aid]
  ronames=dtm$aid
  dtm[,aid:=NULL] # remove unwanted row
  CpGm=as.matrix(dtm) 
  CpGm[CpGm == -1]=NA # put NAs
  rownames(CpGm)=ronames
  return(CpGm)
}