josefin-werme / LAVA

51 stars 9 forks source link

All r2.lower, r2.upper, lo.lower and lo.upper are NA #31

Closed zh-zhang1984 closed 2 years ago

zh-zhang1984 commented 2 years ago

Hi I have run the run.bivar successfully but the results seem some errors because all r2.lower, r2.upper, lo.lower and lo.upper are NA. When I extend the locus to 100, the these quantitatives are still all NAs

library(LAVA)
#generate input info 
filename <- list.files("/Users/zhangzhongheng/Documents/2022/GWAS_sepsis/munge",
                       full.names = TRUE)
filename <- filename[!str_detect(filename,".log$")]
filename <- data.frame(
  id = sub(".*munge/","",filename),
  filename = filename
) %>% 
  mutate(id = sub("\\.sumstats.gz","",id))
input.info <- metaDat %>% 
  select(id,phenotype = trait,cases = ncase,controls=ncontrol) %>% 
  filter(phenotype%in%names(datHeatmap)) %>% 
  left_join(
    filename,
    by = c("id"="id")
  )
write_delim(
  input.info,
  file = "/Users/zhangzhongheng/Documents/2022/GWAS_sepsis/LAVAdata/input.info.txt"
)
#generate sample overlap
write.table(
  datHeatmap,
  file = "/Users/zhangzhongheng/Documents/2022/GWAS_sepsis/LAVAdata/sample.overlap.txt",
  row.names= T
)
input = process.input(
  input.info.file="/Users/zhangzhongheng/Documents/2022/GWAS_sepsis/LAVAdata/input.info.txt",           # input info file
  sample.overlap.file="/Users/zhangzhongheng/Documents/2022/GWAS_sepsis/LAVAdata/sample.overlap.txt",   # sample overlap file (can be set to NULL if there is no overlap)
  ref.prefix="/Users/zhangzhongheng/Documents/2022/GWAS_sepsis/LAVAdata/g1000_eur/g1000_eur",                    # reference genotype data prefix
  phenos=input.info$phenotype
  )  
loci <- read.loci(
  "/Users/zhangzhongheng/Documents/2022/GWAS_sepsis/LAVAdata/blocks_s2500_m25_f1_w200.GRCh37_hg19.locfile"
)

dat_localrg <- NULL
for (ii in 1:nrow(loci[1:5,])) {
  locus = process.locus(
    loci[ii,], input,
    phenos = names(datHeatmap)[1:20]
  )
  dat_localrg1 <- run.bivar(locus)
  dat_localrg <- bind_rows(
    dat_localrg,
    dat_localrg1 %>% mutate(Locus = loci$LOC[ii])
    )
}
head(dat_localrg)
            phen1                                      phen2       rho rho.lower
1 Type 2 diabetes Low density lipoprotein cholesterol levels 0.0805565        NA
2 Type 2 diabetes                              triglycerides 1.0000000        NA
3 Type 2 diabetes                        basophil cell count 0.3251890        NA
4 Type 2 diabetes                     white blood cell count        NA        NA
5 Type 2 diabetes                        monocyte cell count 0.5435930        NA
6 Type 2 diabetes                      lymphocyte cell count 1.0000000        NA
  rho.upper       r2 r2.lower r2.upper         p Locus
1        NA 0.000000       NA       NA 0.9260910     1
2        NA 0.405397       NA       NA       NaN     1
3        NA 0.128226       NA       NA 0.5281700     1
4        NA       NA       NA       NA        NA     1
5        NA 0.301576       NA       NA 0.3887600     1
6        NA 1.000000       NA       NA 0.0716352     1

Are there any hints to solve this issue?

cadeleeuw commented 2 years ago

Hi,

Apologies for the slow response. To diagnose this issue, could you send me (c.a.de.leeuw@vu.nl) the locus object returned from process.locus() for the first couple of loci? Then I can run through the code manually with that as input and figure out what's going on. You can just save them to a file using something like save.image(locus, file=paste0("locus", ii, ".rdata")) right after the process.locus() call in the loop above should do the trick.

Best, Christiaan

zh-zhang1984 commented 2 years ago

I can solve this issue, this time some of the rows contains NAs; and most rows do not contain NA.

shenlan17 commented 1 year ago

I can solve this issue, this time some of the rows contains NAs; and most rows do not contain NA.

Hello, I have also encountered this problem recently. I do not know the reason and how to solve it. Could you give me some hints