R/rareMETALS.single.correctRefAlt.R

In dajiangliu/rareGWAMA: Meta-Analysis of Gene-level Rare Variant Association Test in Whole Genome Datasets

#' Single variant meta-analysis
#'
#' @param score.stat.file files of score statistics
#' @param cov.file covariance matrix files
#' @param range tabix range of variants to be analyzed
#' @param refaltList A list of reference alternative and positions of variants to be analyzed; Each variant in the dataset will be match against ref/alt alleles specified in refaltList; Only variants with matched ref and alt alleles can be included; 
#' @param alternative alternative hypothesis to be specified
#' @param callrate.cutoff Cutoffs of call rate, lower than which will NOT be analyzed (labelled as missing)
#' @param hwe.cutoff Cutoffs of HWE p-values
#' @param refaltList A list that contains pos, ref and alt for calibrating variant sites to have the right ref and alt alleles;
#' @return a list consisting of results
#' @export
rareMETALS.single.correctRefAlt <- function(score.stat.file,cov.file,range,refaltList,alternative=c('two.sided','greater','less'),callrate.cutoff=0,hwe.cutoff=0)
{
    ix.gold <- 1;
    extra.par <- list(ix.gold=ix.gold,QC.par=list(callrate.cutoff=callrate.cutoff,hwe.cutoff=hwe.cutoff));
    capture.output(raw.data.all <- rvmeta.readDataByRange( score.stat.file, cov.file, range));
    if(length(raw.data.all)==0)
        return(list(list(p.value=NA,
                         skip=1,
                         statistic=NA,
                         no.var=0,
                         no.sample=NA)));
    if(length(alternative)>1) alternative <- alternative[1];
    ix.gold <- extra.par$ix.gold;
    if(length(extra.par$ix.gold)==0) ix.gold <- 1;
    raw.data <- raw.data.all[[1]];
    raw.data.ori <- raw.data;
    if(length(extra.par$QC.par)>0) 
        raw.data <- QC(raw.data,extra.par$QC.par,cov=0);
    log.mat <- raw.data$log.mat;
    ix.pop <- 1:length(raw.data$ref);
    
    U.stat <- 0;V.stat.sq <- 0;V.stat.sq <- 0;
    p.value <- 0;statistic <- 0;
    direction.by.study <- "+";
    beta1.est <- 0;
    beta1.sd <- 0;
    hsq.est <- 0;
    no.sample <- 0;
    maf.vec <- 0;
    no.sample.pop <- 0;
    no.sample.var <- 0;
    ustat <- 0;
    vstat <- 0;
    nref <- 0;nalt <- 0;nhet <- 0;
    no.sample.mat <- matrix(0,nrow=length(ix.pop),ncol=length(raw.data$nSample[[1]]));
    ix.var <- 1;nearby <- character(0);
    if(length(raw.data$ref[[ix.gold]])==0)
        return(list(p.value=NA,
                    ref=NA,
                    alt=NA,
                    raw.data=raw.data.ori,
                    raw.data.clean=raw.data,
                    statistic=NA,
                    direction.by.study=NA,
                    anno=NA,
                    maf=NA,
                    no.sample=NA,
                    beta1.est=NA,
                    beta1.sd=NA,
                    hsq.est=NA,
                    nearby=NA,
                    pos=NA));
    ix.match <- match(raw.data$pos,refaltList$pos);
    refaltList <- list(pos=refaltList$pos[ix.match],ref=refaltList$ref[ix.match],alt=refaltList$alt[ix.match]);
    ref.gold <- refaltList$ref;
    alt.gold <- refaltList$alt;
    maf.sd.vec <- 0;maf.maxdiff.vec <- 0;ix.maf.maxdiff.vec <- 0;
    maf.pop.ori <- 0;
    QC.by.study <- "";
    for(ix.var in 1:length(raw.data$ref[[ix.gold]]))
    {
        if(ix.var/1000==as.integer(ix.var/1000))
        {
        }
        direction.by.study.var <- rep("^",length(ix.pop));
        U.stat <- 0;V.stat.sq <- 0;maf.pop <- 0;
        nref.var <- 0;nhet.var <- 0;nalt.var <- 0;
        ix.include <- rep(0,length(ix.pop));
        for(ii in 1:length(ix.pop))
        {
            
            if(!is.na(raw.data.ori$ref[[ii]][ix.var]) | !is.na(raw.data.ori$alt[[ii]][ix.var]))
                if(rm.na(raw.data$af[[ii]][ix.var])==0 | rm.na(raw.data$af[[ii]][ix.var])==1)
                {
                    U.stat <- U.stat+rm.na(raw.data$ustat[[ii]][ix.var]);
                    V.stat.sq <- V.stat.sq+(rm.na(raw.data$vstat[[ii]][ix.var]))^2;
                    nref.var <- nref.var+rm.na(raw.data$nref[[ii]][ix.var]);
                    nalt.var <- nalt.var+rm.na(raw.data$nalt[[ii]][ix.var]);
                    nhet.var <- nhet.var+rm.na(raw.data$nhet[[ii]][ix.var]);
                    no.sample.mat[ii,ix.var] <- raw.data.ori$nalt[[ii]][ix.var]+raw.data.ori$nref[[ii]][ix.var]+raw.data.ori$nhet[[ii]][ix.var]
                    if(is.na(no.sample.mat[ii,ix.var])) no.sample.mat[ii,ix.var] <- raw.data.ori$nSample[[ii]][ix.var];
                    ix.include[ii] <- 1;
                    log.mat[ix.var,ii] <- paste(log.mat[ix.var,ii],"0",sep="",collapse="");
                }           
        }
        if(!is.na(raw.data$ref[[ii]][ix.var]) & !is.na(raw.data$alt[[ii]][ix.var]))
        {
            flip.ref.alt <- (raw.data$ref[[ii]][ix.var]==refaltList$alt[ix.var] & raw.data$alt[[ii]][ix.var]==refaltList$ref[ix.var]);
            match.ref.alt <- (refaltList$ref[ix.var]==(raw.data$ref[[ii]][ix.var]) & refaltList$alt[ix.var]==(raw.data$alt[[ii]][ix.var]));
            mono <- (((raw.data$ref[[ii]][ix.var]==".") | (raw.data$ref[[ii]][ix.var]==0) | (raw.data$alt[[ii]][ix.var]==".") | (raw.data$alt[[ii]][ix.var]==0) | (rm.na(raw.data$af[[ii]][ix.var])==0)| rm.na(raw.data$af[[ii]][ix.var])==1) & !is.na(raw.data.ori$nSample[[ii]][ix.var]))
            if(flip.ref.alt)
            {
                raw.data$af[[ii]][ix.var] <- 1-raw.data$af[[ii]][ix.var];
                raw.data$ustat[[ii]][ix.var] <- (-1)*(raw.data$ustat[[ii]][ix.var]);
                U.stat <- U.stat+rm.na(raw.data$ustat[[ii]][ix.var]);
                nref.var <- nref.var+rm.na(raw.data$nalt[[ii]][ix.var]);
                nalt.var <- nalt.var+rm.na(raw.data$nref[[ii]][ix.var]);
                nhet.var <- nhet.var+rm.na(raw.data$nhet[[ii]][ix.var]);
                V.stat.sq <- V.stat.sq+(rm.na(raw.data$vstat[[ii]][ix.var]))^2;                    
                no.sample.mat[ii,ix.var] <- raw.data$nalt[[ii]][ix.var]+raw.data$nref[[ii]][ix.var]+raw.data$nhet[[ii]][ix.var]
                if(is.na(no.sample.mat[ii,ix.var])) no.sample.mat[ii,ix.var] <- raw.data$nSample[[ii]][ix.var];
                ix.include[ii] <- 1;
                log.mat[ix.var,ii] <-  paste(log.mat[ix.var,ii],"F",sep="",collapse="");
            }
            if(match.ref.alt)
            {
                U.stat <- U.stat+rm.na(raw.data$ustat[[ii]][ix.var]);
                V.stat.sq <- V.stat.sq+(rm.na(raw.data$vstat[[ii]][ix.var]))^2;
                nref.var <- nref.var+rm.na(raw.data$nref[[ii]][ix.var]);
                nalt.var <- nalt.var+rm.na(raw.data$nalt[[ii]][ix.var]);
                nhet.var <- nhet.var+rm.na(raw.data$nhet[[ii]][ix.var]);
                no.sample.mat[ii,ix.var] <- raw.data$nalt[[ii]][ix.var]+raw.data$nref[[ii]][ix.var]+raw.data$nhet[[ii]][ix.var]
                if(is.na(no.sample.mat[ii,ix.var])) no.sample.mat[ii,ix.var] <- raw.data$nSample[[ii]][ix.var];
                ix.include[ii] <- 1;
                log.mat[ix.var,ii] <-  paste(log.mat[ix.var,ii],"M",sep="",collapse="");
                
            }
            if(mono)
            {
                ix.include[ii] <- 1;
            }
            if(!match.ref.alt & !flip.ref.alt & !mono)
            {
                ix.include[ii] <- 1;
                log.mat[ix.var,ii] <- paste(log.mat[ix.var,ii],"O",sep="",collapse="");
                raw.data$nSample[[ii]][ix.var] <- NA;
                raw.data$af[[ii]][ix.var] <- NA;
                raw.data$ac[[ii]][ix.var] <- NA;
                raw.data$nref[[ii]][ix.var] <- NA;
                raw.data$nhet[[ii]][ix.var] <- NA;
                raw.data$nalt[[ii]][ix.var] <- NA;
            }
            if(is.na(raw.data$ustat[[ii]][ix.var])) direction.by.study.var[ii] <- "^";
            if(!is.na(raw.data$ustat[[ii]][ix.var]))
            {
                if(raw.data$ustat[[ii]][ix.var]>0) direction.by.study.var[ii] <- "+";
                if(raw.data$ustat[[ii]][ix.var]<0) direction.by.study.var[ii] <- "-";
                if(raw.data$ustat[[ii]][ix.var]==0) direction.by.study.var[ii] <- "?";
            }
        }
        maf.pop[ii] <- ((raw.data$af[[ii]])[ix.var]);
        maf.pop.ori[ii] <- (raw.data$af[[ii]])[ix.var];
        
        maf.vec[ix.var] <- sum(maf.pop*no.sample.mat[,ix.var],na.rm=TRUE)/sum(no.sample.mat[,ix.var],na.rm=TRUE);
        maf.sd.vec[ix.var] <- sqrt(sum((no.sample.mat[,ix.var])*((maf.pop.ori-maf.vec[ix.var])^2),na.rm=TRUE)/sum(no.sample.mat[,ix.var],na.rm=TRUE));
        maf.maxdiff.tmp <- max(abs(maf.pop.ori-maf.vec[ix.var]),na.rm=TRUE);
        maf.maxdiff.vec[ix.var] <- 0;ix.maf.maxdiff.vec[ix.var] <- sample(1:length(maf.vec),1);
        if(maf.maxdiff.tmp>0)
        {
            maf.maxdiff.vec[ix.var] <- maf.maxdiff.tmp;
            ix.maf.maxdiff.vec[ix.var] <- which.max(abs(maf.pop.ori-maf.vec[ix.var]));
        }
        ix.0 <- which(maf.pop==0 | maf.pop==1);
        if(length(ix.0)>0)
        {
            maf.sd.vec[ix.var] <- sqrt(sum((no.sample.mat[-ix.0,ix.var])*((maf.pop[-ix.0]-maf.vec[ix.var])^2),na.rm=TRUE)/sum(no.sample.mat[,ix.var],na.rm=TRUE));
            maf.maxdiff.vec[ix.var] <- 0;ix.maf.maxdiff.vec[ix.var] <- ix.0[1];
            if(max(abs(maf.pop[-ix.0]-maf.vec[ix.var]),na.rm=TRUE)>0)
            {
                maf.maxdiff.vec[ix.var] <- max(abs(maf.pop[-ix.0]-maf.vec[ix.var]),na.rm=TRUE);
                ix.rm0 <- (1:length(maf.pop))[-ix.0];
                ix.maf.maxdiff.vec[ix.var] <- ix.rm0[which.max(abs(maf.pop[-ix.0]-maf.vec[ix.var]))];
            }
            
        }
        direction.by.study[ix.var] <- paste(direction.by.study.var,sep='',collapse='');
        QC.by.study[ix.var] <- paste(log.mat[ix.var,],sep=",",collapse=",");
        beta1.est[ix.var] <- U.stat/V.stat.sq;
        beta1.sd[ix.var] <- sqrt(1/V.stat.sq);
        hsq.est[ix.var] <- (beta1.est[ix.var])*(beta1.est[ix.var])*V.stat.sq/sum(no.sample.mat[,ix.var],na.rm=TRUE);
        no.sample.var[ix.var] <- sum(no.sample.mat[,ix.var],na.rm=TRUE);
        nearby[ix.var] <- paste(get.gene.inWindow(raw.data$pos[ix.var],1e6),sep=",",collapse=",");
        if(alternative=='two.sided'){
            statistic[ix.var] <- U.stat^2/V.stat.sq;
            ustat[ix.var] <- U.stat;
            vstat[ix.var] <- sqrt(V.stat.sq);
            nref[ix.var] <- nref.var;
            nalt[ix.var] <- nalt.var;
            nhet[ix.var] <- nhet.var;
            ix.inf <- which(V.stat.sq==0);
            if(length(ix.inf)>0) statistic[ix.var][ix.inf] <- NA;
            p.value[ix.var] <- pchisq(statistic[ix.var],df=1,lower.tail=FALSE)
        }
        if(alternative=='greater'){
            statistic[ix.var] <- (U.stat/sqrt(V.stat.sq));
            p.value[ix.var] <- pnorm(statistic[ix.var],lower.tail=FALSE)
        }
        if(alternative=='less'){
            statistic[ix.var] <- (U.stat/sqrt(V.stat.sq));
            p.value[ix.var] <- pnorm(statistic[ix.var],lower.tail=TRUE)
        }
    }
    
    ref.out <- list(ref.gold);
    alt.out <- list(alt.gold);
    nSample.out <- list(no.sample.var);
    af.out <- list(maf.vec);
    ac.out <- list(as.integer(2*maf.vec*no.sample.var));
    callrate <- list(rep(1,length(ref.gold)));
    nref.out <- list(nref);
    nalt.out <- list(nalt);
    nhet.out <- list(nhet);
    ustat.out <- list(ustat);
    vstat.out <- list(vstat);
    callrate.out <- list(rep(1,length(maf.vec)));
    hwe.out <- list(rep(1,length(maf.vec)));
    effect.out <- list(beta1.est);
    pVal.out <- list(p.value);
    pos.out <- raw.data$pos;
    anno.out <- raw.data$anno;
    cov.out <- list(matrix(nrow=0,ncol=0));
    integratedData <- list();
    integratedData[[1]] <- list(ref=ref.out,
                                alt=alt.out,
                                nSample=nSample.out,
                                af=af.out,
                                ac=ac.out,
                                callrate=callrate.out,
                                hwe=hwe.out,
                                nref=nref.out,
                                nalt=nalt.out,
                                ustat=ustat.out,
                                vstat=vstat.out,
                                nhet=nhet.out,
                                effect=effect.out,
                                pVal=pVal.out,
                                cov=cov.out,
                                pos=pos.out,
                                covZZ=list(matrix(nrow=0,ncol=0)),
                                covXZ=list(matrix(ncol=0,nrow=length(ref.out[[1]]))),
                                hweCase=list(rep(NA,length(ref.out[[1]]))),
                                hweCtrl=list(rep(NA,length(ref.out[[1]]))),
                                afCtrl=list(rep(NA,length(ref.out[[1]]))),
                                afCase=list(rep(NA,length(ref.out[[1]]))),
                                anno=anno.out);
    
    
    return(list(p.value=p.value,
                ref=ref.gold,
                alt=alt.gold,
                integratedData=integratedData,
                raw.data=raw.data.ori,
                clean.data=raw.data,
                statistic=statistic,
                direction.by.study=direction.by.study,
                anno=raw.data$anno,
                maf=maf.vec,
                maf.maxdiff.vec=maf.maxdiff.vec,
                ix.maf.maxdiff.vec=ix.maf.maxdiff.vec,
                maf.sd.vec=maf.sd.vec,
                no.sample.mat=no.sample.mat,
                no.sample=no.sample.var,
                beta1.est=beta1.est,
                beta1.sd=beta1.sd,
                QC.by.study=QC.by.study,
                hsq.est=hsq.est,
                nearby=nearby,
                log.mat=log.mat,
                pos=raw.data$pos));
}