tshmak/lassosum

!any(is.na(cor)) is not TRUE

maltose1117 opened this issue · 7 comments

Hi, below is my code.

library(lassosum)
library(data.table)
library(parallel)
library(fdrtool)
cl <- makeCluster(7, type="FORK") 

#the data is downloaded from https://blog.nus.edu.sg/agen/summary-statistics/t2d-2020/ T2D_ALL_PRIMARY dataset
ss <- fread("/path/to/SpracklenCN_prePMID_T2D_ALL_Primary.txt") 

ss$P[ss$P < .Machine$double.xmin] <- .Machine$double.xmin

test.bfile <- "bi.merged"
cor <- p2cor(p = ss$P, n = ss$Neff, sign=ss$Beta)
sum(is.na(cor))

ld <- fread("/path/to/Berisa.ASN.hg19.bed")
out <- lassosum.pipeline(cor=cor, chr=ss$Chr, pos = ss$Pos, A1 = ss$EAF, A2 = ss$NEA, 
                         test.bfile = "/path/to/bi.merged",
                         LDblocks = ld, cluster = cl,exclude.ambiguous = F)

v <- pseudovalidate(out)
save(v,file="/path/to/ASN.DM.Rdata")

The output is:

[1] 0
Coordinating summary stats with reference panel...
Splitting genome by LD blocks ...
Running lassosum ...
s =  0.2
Error in lassosum(cor = cor2, bfile = ref.bfile, shrink = s, extract = ref.extract,  :
  !any(is.na(cor)) is not TRUE
Calls: lassosum.pipeline -> lapply -> FUN -> lassosum -> stopifnot
Execution halted

I managed to use the same "bi.merged" bfiles to construct PRS for other phenotypes but not this one. Sorry that I may not be able to provide the bfiles I used. But if you really need them to diagnose the problem, maybe I could revise the .fam data a bit and send them to you.

Hope my question is clear now. Thank you for your time.

A1 = ss$EAF seems odd, as "EAF" is often "European Allele Frequency", which is not the alleles.
You should probably verify what you're passing to A1 and A2.

A1 = ss$EAF seems odd, as "EAF" is often "European Allele Frequency", which is not the alleles.
You should probably verify what you're passing to A1 and A2.

Thanks for your comment. According to the summary statistics website (https://blog.nus.edu.sg/agen/summary-statistics/t2d-2020/), EAF refers to effect allele frequency, and I think this is what I shall pass to 'A1'.

A1 is the alternative allele for each variant, it should be a character vector mainly composed of "A", "T", "C" and "G".
This "EA" in your case I think, not "EAF".

A1 is the alternative allele for each variant, it should be a character vector mainly composed of "A", "T", "C" and "G".

Seems I made a silly mistake...You are right, I will use

A1 = ss$EA

and see the result. Thanks!

Now the output becomes:

[1] 0
Coordinating summary stats with reference panel...
Splitting genome by LD blocks ...
Running lassosum ...
s =  0.2
s =  0.5
s =  0.9
Running lassosum with s=1...
Calculating polygenic scores ...

And another error happens:

  trying to use CRAN without setting a mirror
Calls: install.packages -> startsWith -> contrib.url
Execution halted

Shall I open a new issue to discuss it? Thanks for your comments in advance.

This is weird. lassosum does not call install.packages. Are you sure the error does not come from your own code?

It's my fault. I forgot to # the install.packages("fdrtool") in the pseudovalidation part. Thank you both!