R/fitAll.R

In GPA: GPA (Genetic analysis incorporating Pleiotropy and Annotation)

# fit GPA for all possible combination

fitAll <- function( pmat,
  maxIter=2000, stopping="relative", epsStopLL=1e-10,
  parallel=FALSE, nCore=8 ) {

  # consider all possible cases

  combs <- t(combn( seq_len(ncol(pmat)), 2 ))
  comblist <- split( combs, seq_len(nrow(combs)) )


  # fit GPA

  message("-------------------------------------------------------------")
  message( "Fitting GPA for all possible combination of GWAS datasets...", sep="" )
  message("-------------------------------------------------------------")

  if ( parallel==TRUE ) {
    fit.list.GPA <- mclapply( comblist, function(cc) {
      message( "\tGWAS pair: ",cc[1]," and ",cc[2], sep="" )
      fit.GPA.cc <- GPA( pmat[ , c(cc[1],cc[2]) ], NULL, pleiotropyH0=FALSE,
        maxIter=maxIter, stopping=stopping, epsStopLL=epsStopLL, verbose=0 )
      return(fit.GPA.cc)
    }, mc.cores=nCore )
  } else {
    fit.list.GPA <- lapply( comblist, function(cc) {
      message( "\tGWAS pair: ",cc[1]," and ",cc[2], sep="" )
      fit.GPA.cc <- GPA( pmat[ , c(cc[1],cc[2]) ], NULL, pleiotropyH0=FALSE,
        maxIter=maxIter, stopping=stopping, epsStopLL=epsStopLL, verbose=0 )
      return(fit.GPA.cc)
    })
  }

  message("")
  message("Done!")
  message("")


  # fit GPA under H0

  message("-------------------------------------------------------------")
  message( "Fitting GPA for all possible combination under H0...", sep="" )
  message("-------------------------------------------------------------")

  if ( parallel==TRUE ) {
    fit.list.GPA.H0 <- mclapply( comblist, function(cc) {
      message( "Fit GPA for the GWAS ",cc[1]," and ",cc[2], sep="" )
      fit.GPA.H0.cc <- GPA( pmat[ , c(cc[1],cc[2]) ], NULL, pleiotropyH0=TRUE,
        maxIter=maxIter, stopping=stopping, epsStopLL=epsStopLL, verbose=0 )
      return(fit.GPA.H0.cc)
    }, mc.cores=nCore )
  } else {
    fit.list.GPA.H0 <- lapply( comblist, function(cc) {
      message( "Fit GPA for the GWAS ",cc[1]," and ",cc[2], sep="" )
      fit.GPA.H0.cc <- GPA( pmat[ , c(cc[1],cc[2]) ], NULL, pleiotropyH0=TRUE,
        maxIter=maxIter, stopping=stopping, epsStopLL=epsStopLL, verbose=0 )
      return(fit.GPA.H0.cc)
    })
  }

  message("")
  message("Done!")
  message("")


  # calculate pleiotropy test p-values

  message("-------------------------------------------------------------")
  message( "Calculating pleiotropy test p-values...", sep="" )
  message("-------------------------------------------------------------")

  plist <- lapply( seq_len(length(fit.list.GPA)), function(i) {
  	ll.H1 <- get_fit(fit.list.GPA[[i]])$loglik[ length(get_fit(fit.list.GPA[[i]])$loglik) ]
  	ll.H0 <- get_fit(fit.list.GPA.H0[[i]])$loglik[ length(get_fit(fit.list.GPA.H0[[i]])$loglik) ]

  	LRT <- -2 * ( ll.H0 - ll.H1 )
  	pval <- pchisq( LRT, 1, lower.tail=FALSE )
    return(pval)
  } )

  pvec <- unlist(plist)
  pTestPval <- matrix( 0, ncol=ncol(pmat), nrow=ncol(pmat) )
  rownames(pTestPval) <- colnames(pmat)
  colnames(pTestPval) <- colnames(pmat)

  for ( i in seq_len(nrow(combs)) ) {
    pTestPval[ combs[i,1], combs[i,2] ] <- pvec[i]
    pTestPval[ combs[i,2], combs[i,1] ] <- pvec[i]
  }

  message("Done!")

  return(list( pmat=pmat, combs=combs, combList=comblist,
    pTestPval=pTestPval, fitGPA=fit.list.GPA, fitH0=fit.list.GPA.H0 ))
}