search

josefin-werme / LAVA

54 stars 9 forks source link

r2.lower, r2.upper, lo.lower and lo.upper are NA #48

Closed shenlan17 closed 1 year ago

shenlan17 commented 1 year ago

Hi I run the code as follows:

library(LAVA) input = process.input(input.info.file="/gpfsnew/lab/groupYU/members/liumengge/test/LAVA/info.txt", # input info file sample.overlap.file=NULL, # sample overlap file (can be set to NULL if there is no overlap) ref.prefix="/gpfs/lab/groupYU/members/liumengge/g1000_eur/g1000_eur", # reference genotype data prefix phenos=c("1499","1498","1497")) # subset of phenotypes listed in the input info file that we want to process loci = read.loci("/gpfs/lab/groupYU/members/liumengge/anaconda/anaconda3/lib/R/library/LAVA/blocks_s2500_m25_f1_w200.GRCh37_hg19.locfile") n.loc = nrow(loci) univ.p.thresh = 0.05/2495 print(paste("Starting LAVA analysis for",n.loc,"loci")) u=b=list() t1=proc.time() for (i in 1:50) { locus = process.locus(loci[i,], input)
if (!is.null(locus)) { loc.info = data.frame(locus = locus$id, chr = locus$chr, start = locus$start, stop = locus$stop, n.snps = locus$n.snps, n.pcs = locus$K) loc.out = run.univ.bivar(locus, univ.thresh = univ.p.thresh) u[[i]] = cbind(loc.info, loc.out$univ) if(!is.null(loc.out$bivar)) b[[i]] = cbind(loc.info, loc.out$bivar) } } t2=proc.time() print(paste("time:",t2[3]-t1[3])) write.table(do.call(rbind,u), paste0(out.fname,".univ.lava"), row.names=F,quote=F,col.names=T) write.table(do.call(rbind,b), paste0(out.fname,".bivar.lava"), row.names=F,quote=F,col.names=T)

and in the file of binary analyses results , r2.lower, r2.upper, lo.lower and lo.upper are all NA image I don't know why and how to resolve it, could you give me a hint? Thanks!

cadeleeuw commented 1 year ago

Hi,

The most likely reason for this is that the rho estimates for those loci / phenotype pairs were too far out of bounds, which causes the estimates, CI bounds and p-values to all be set to NA. There should be a notification in the screen output about this for every occurrence.

You can manually change the cut-off value used to determine when estimates are too far out of bounds (the param.lim argument), though unless they are quite close to the default cut-off this is not generally recommended. If correlation estimates are very far out of bounds, this usually means they are just too unstable to be reliably estimated and tested.

Best, Christiaan