bcm-uga / pcadapt

Performing highly efficient genome scans for local adaptation with R package pcadapt v4
https://bcm-uga.github.io/pcadapt
37 stars 10 forks source link

BLAS/LAPACK routine 'DSYEVR' gave error code -9 #9

Closed V-Z closed 6 years ago

V-Z commented 6 years ago

Hello, I ran the analysis like this:

DATA.pcadapt.data <- read.pcadapt(input.filename="DATA.vcf.gz", type="vcf", local.env=FALSE, ploidy=2,  pop.sizes=NULL, blocksize=100000)
Meta line 55413 read in.
All meta lines processed.
Character matrix gt created.
Character matrix gt rows: 12559
Character matrix gt cols: 128
skip: 0
nrows: 12559
row_num: 0

Processed variant: 12559
All variants processed
 12559 variants have been processed.
 0 variants have been discarded as they are not SNPs.

minua.pcadapt <- pcadapt(input=DATA.pcadapt.data, K=35, method="mahalanobis", data.type="genotype", min.maf=0.05, ploidy=2, output.filename="DATA.pcadapt", clean.files=TRUE, cover.matrix=NULL)
Reading file DATA.vcf.gz...
Number of SNPs: 33222
Number of individuals: 3
Number of SNPs with minor allele frequency lower than 0.05 ignored: 33222
Error in pcadapt(input = minua.pcadapt.data, K = 35, method = "mahalanobis",  : 
  BLAS/LAPACK routine 'DSYEVR' gave error code -9

This is weird as the genind object greated from the input VCF looks correctly like this:

/// GENIND OBJECT /////////

 // 116 individuals; 12,559 loci; 25,118 alleles; size: 17.8 Mb

 // Basic content
   @tab:  116 x 25118 matrix of allele counts
   @loc.n.all: number of alleles per locus (range: 2-2)
   @loc.fac: locus factor for the 25118 columns of @tab
   @all.names: list of allele names for each locus
   @ploidy: ploidy of each individual  (range: 2-2)
   @type:  codom
   @call: .local(x = x, i = i, j = j, drop = drop)

 // Optional content
   @pop: population of each individual (group size range: 1-10)

So there is something very wrong with pcadapt. Any idea, what is wrong, please?

keurcien commented 6 years ago

Hi,

Thanks for reporting the issue. It seems that the issue comes from the converter vcf2pcadapt (which we are recoding currently due to some inconsistencies). It shouldn't have displayed

Reading file DATA.vcf.gz...

but instead it should have shown

Reading file DATA.pcadapt...

That means that pcadapt is running on the vcf.gz file instead of the converted file. The issue probably comes from the fact that I never considered running pcadapt on vcf.gz, so the output file name is bugged.

So from what I understand, there are two options:

  1. Do you have any files in your directory that ends with .pcadapt (or that contains .pcadapt in its name) that has been created by read.pcadapt, let's say filename.pcadapt ? If yes, run pcadapt on this file:
pcadapt("filename.pcadapt", K = 35, ...)
  1. If no, extract the vcf.gz to a standard vcf file using gunzip, and rerun read.pcadapt on the vcf file.

We are recoding pcadapt so any suggestion (or bug report) is welcome. Tell me if this worked for you.

Cheers,

Keurcien

V-Z commented 6 years ago

Hello, thank You for fast reply. Interestingly, If I run read.pcadapt as in my 1st comment, the R object DATA.pcadapt.data contains:

> DATA.pcadapt.data
[1] "DATA.vcf.gz"

pcadat creates in the working directory file tmp.pcadapt. If I run same read.pcadapt on uncompressed VCF, I see

> DATA.pcadapt.data
[1] "DATA.pcadapt"

but next step fails:

> DATA.pcadapt <- pcadapt(input=DATA.pcadapt.data, ...)
Error in pcadapt(input = DATA.pcadapt.data, K = 35, method = "mahalanobis",  : 
  File DATA.pcadapt does not exist.

In the working directory, there is only file tmp.pcadapt. If I run pcadapt() on that file like pcadapt(input="tmp.pcadapt", ...) it works well.

keurcien commented 6 years ago

Ok thanks for the feedback, sorry for the inconvenience. I'm kinda relieved as it is less a hassle to debug. But still, thanks for pointing that out. So yeah you should use the tmp.pcadapt for the moment, until next version (we expect the new release to come out in January 2018).

keurcien commented 6 years ago

I'm closing this issue since the new version has been released.