R/get_peak_blocks_cor.R

In jaspershen/MSannotator: xMSannotator: R package for network based annotation of metabolomics data.

get_peak_blocks_cor <-
function(dataA, max.rt.diff=10, alignment.tool="apLCMS", clust.method="correlation",cor.method="pearson",corthresh=0.7,corpvaluethresh=0.05,time_step=3,mult.test.cor=FALSE)
{	
        diff_mz_num=1
        unique_mz={}
        ppm_v={}
        rt_v={}
        cnames=colnames(dataA)
        dataA<-as.data.frame(dataA)
	
		#dataA<-as.data.frame(dataA,2,as.numeric)
         if(alignment.tool=="apLCMS")
        {
              sample.col.start=2
        }
        else
        {
              if(alignment.tool=="XCMS")
              {
                    sample.col.start=2
                    cnames[1]="mz"

                    cnames[2]="time"
                    colnames(dataA)=cnames
                    
              }
              else
              {
              		 sample.col.start=3
                    cnames[1]="mz"

                    cnames[2]="time"
                    colnames(dataA)=cnames
                    #stop(paste("Invalid value for alignment.tool. Please use either \"apLCMS\" or \"XCMS\"", sep=""))
              }
        }

	
	if(clust.method=="correlation"){
				 cl<-makeCluster(10)

        #clusterExport(cl, "getCorchild")

        #goodfeats[mnum,-c(1:2)]),data_mt

        #system.time(pearson_res<-parRapply(cl,dataA[,-c(1:sample.col.start)],getCorchild,t(dataA[,-c(1:sample.col.start)]),cor.method))


        #stopCluster(cl)
		 pearson_resmat<-{}
		 
		 
		 #return(list(cormat=cormat,complete_pearsonpvalue_mat=complete_pearsonpvalue_mat,complete_pearsonqvalue_mat=complete_pearsonqvalue_mat))
	
			pearson_Res_all<-{}

			cormat<-WGCNA::cor(t(dataA[,-c(1:sample.col.start)]), use = 'p')
		    cormat<-as.data.frame(cormat)
        cormat<-as.matrix(cormat)
        
           
      		cor_list<-{}
		
		cormat<-cor2pcor(cormat)

		cor_groups<-new("list")
		
		#for(cm in 1:dim(cormat)[1]){
			
			cor_groups<-lapply(1:dim(cormat)[1],function(cm){
			corind_list=new("list")
			cor_ind<-which(cormat[,cm]>=corthresh)
			corind_list[[diff_mz_num]]=cor_ind #dataA[getbind_same,]
            diff_mz_num=diff_mz_num+1
			return(corind_list)
			})
		
		cor_group_list<-{}
		#
	  for(m in 1:length(cor_groups)){
	  
	   if((m%in%cor_group_list)==FALSE){
		for(n in (m+1):length(cor_groups)){
		
			if(n>length(cor_groups)){
				break;
			}
			com1<-intersect(cor_groups[[m]][[1]],cor_groups[[n]][[1]])
			if(length(com1)>0){
				
				cor_groups[[m]][[1]]<-c(cor_groups[[m]][[1]],cor_groups[[n]][[1]])
				cor_groups[[m]][[1]]<-unique(cor_groups[[m]][[1]])
				cor_group_list<-c(cor_group_list,n)
				cor_groups[[n]][[1]]<-c(0)
			}
		}
			cor_groups[[m]][[1]]<-unique(cor_groups[[m]][[1]])
		
		}
	  }
	  
	  if(length(cor_group_list)>0){
		cor_groups<-cor_groups[-cor_group_list]
		}
	 
	  diffmat<-{}
	cor_groups<-unique(cor_groups)
	levelAnum<-1
	 for(m in 1:length(cor_groups)){
	 	if((m%in%cor_group_list)==FALSE){
	 	
	 	
	 	cur_group_ind<-unique(cor_groups[[m]][[1]])
	 	
	 	cur_group<-dataA[cur_group_ind,]
	 	cur_group<-cbind(levelAnum,cur_group)
	 	
	 	diffmat<-rbind(diffmat,cur_group)
	 	levelAnum<-levelAnum+1
	 	}
	 	}
	

	#cor.test(as.numeric(cur_group[1,-c(1:10)]),as.numeric(cur_group[10,-c(1:10)]))

	}else{
		
		if(clust.method=="wgcna"){
			
				levelAnum<-mlabels$colors
				diffmat<-cbind(levelAnum,dataA)
				diffmat<-diffmat[order(diffmat$levelAnum),]
				
				for(m1 in mlabels){
					subdata<-diffmat[which(mlabels==m1),]
					
				}
				
		}
		
	}
	#b1<-hist(diffmat$time,breaks=seq(0,max(diffmat$time),max.rt.diff))
	
	
	diffmat<-as.data.frame(diffmat)
	
	
	mod_list<-diffmat$levelAnum
	
	t1<-table(mod_list)
mod_names<-names(t1)
mod_names<-as.numeric(mod_names)


 
time_mult_fact<-1


 
 diffmatC<-{}

#t1<-seq(from=1,to=max(dataA$time),by=time_step/2)
 #t1<-hist(subdata$time,breaks=seq(from=1,to=max(subdata$time),by=10/2))
 
# sigma<-min(sd(subdata$time),(quantile(subdata$time)[4]-quantile(subdata$time)[2])/1.34)
 #s1<-min(sqrt(var(x)), h)
 
 #bw1<-1.06 * sigma * length(subdata$time)^(-1/5)
 
 d1<-density(dataA$time,bw="nrd",from=min(dataA$time),to=(10+max(dataA$time)))
 
 
 time_step<-abs(d1$x[1]-d1$x[time_step+1])
 
 t1<-d1$x
 
 time_step<-1*time_step

diffmatB<-{}
for(i in 1:length(mod_names)){

#groupB_res<-sapply(1:length(mod_names),function(i){
	
	groupA_num<-mod_names[i]
	
	subdata<-dataA[which(m2$colors==groupA_num),]
	
 subdata<-subdata[order(subdata$time),]

 groupB<-group_by_rt(subdata,time_step,max.rt.diff=max.rt.diff,groupnum=groupA_num)
	 
	  diffmatB<-rbind(diffmatB,groupB)
	  
	  }
	  
	 return(diffmatB)	
	
       
                               												
			        #return(diffmat)
}