hellenthal-group-UCL / fastGLOBETROTTER

faster haplotype-based admixture inference and dating software
Other
8 stars 1 forks source link

Error in "findpeak()" function - potentially derived from error when defining "wsize" #2

Open 00izquierdoguillem opened 1 year ago

00izquierdoguillem commented 1 year ago

Hi,

I recently ran fastGLOBETROTTER according to the guidelines and examples provided in this github repository ("FastGLOBETROTTER_Instructions.pdf" and "FastGLOBETROTTER_Tutorial.pdf"). When using my own data, depending on the input files used I encountered the following error:

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 0 (non-NA) cases Calls: findpeak -> lm -> lm.fit

Eventually, I traced down this error to the chunk of code defining findpeak() that starts around line 160.

 findpeak  = function  (ourmix,means,bin.width)
 {
   dom = which(ourmix ==max(ourmix))
   new = means[length(ourmix)*(dom-1)+dom,]
   ##if (dim(means)[2] > 100) {wsize = 100  ## !!!!! pongsakorn orig !!!!!!!!
   ##} else  {wsize = round(dim(means)[2]/4)}  ## !!!!! pongsakorn orig !!!!!!!!
   wsize = round(dim(means)[2]/2) ## !!!!! garrett change (now can remove up to half the original -- fitted -- x-axis) !!!!!!!!
 ...
     while (j < wsize && ((j+i)<(dim(means)[2]-10)))  ## !!!!! garrett change (now must keep the last 10 gridpoints at least) !!!!!!!!
     {
       temp = lm(new[j:(j+i)]~c(j:(j+i)))
       ...
      }
 }

In my case, the new vector defined at the start of this function only had non-NA values before a certain element (in my case [,75]), yet wsize was set at roughly 140. This meant that temp could only be defined until the 75th element of new before lm() generated the aforementioned error by pulling the NA value found in the 76th element of new.

I ended up working around this issue by using part of the initial code by Pongsakorn found in the commented line 163: I modified the definition of wsize so that the denominator was 4 instead of 2 in line 164:

 wsize = round(dim(means)[2]/4) 

This solved the issue for me, but I would be grateful if you could take a look at it to see whether this correction was wrong.

Thanks in advance,

Guillem