jgx65 / hierfstat

the hierfstat package
24 stars 14 forks source link

pairwise.WCfst Error - non-conformable arrays #66

Closed silknets closed 1 year ago

silknets commented 1 year ago

Hello all,

I'm having an issue similar to closed post here, but the solutions there aren't sufficient to resolve my issue. I'm trying to generate 2x pairwise FST matrices, one with original microsatellite genotype data (A), and another where that same genotype data has been recoded so that no populations share alleles (B). I've been able to generate the FST matrices for other datasets (original and recoded versions), but for the current dataset, I have only been able to generate the matrix for the original dataset. For some reason, the recoded dataset throws an error no matter what I try! I have use genind, fasta (.dat), and a dataframe (organized like data(gtrunchier)) for inputs, all with the same error:

recode_pair <- pairwise.WCfst(recode) Error in 2 nal p - mho : non-conformable arrays

Any suggestions for how to resolve this error? Help would be much appreciated! I've attached zipped .dat files for both the original and recoded data.

Sam S. gafa.zip gafa_recode.zip

jgx65 commented 1 year ago

Hi Sam,

There is indeed an issue, it has to do with correctly splitting the alleles. You choose an encoding with 3 digits per allele, but none of your alleles are larger than 99. If you add 100100 to all genotypes, things work:

 d1<-read.fstat("gafa.dat")
d2<-read.fstat("gafa_out.dat")
d3<-d2
d3[,-1]<-d3[,-1]+100100 #the fix
genet.dist(d1,method="WC84") #works
genet.dist(d2,method="WC84") #does not work
genet.dist(d3,method="WC84") #works

I'll work on fixing this, thanks for reporting

jgx65 commented 1 year ago

@silknets ,

The issue should be resolved (see code)

> round(genet.dist(d2,method="WC84"),digits=3)
       1     2     3     4     5     6     7     8     9
2  0.712                                                
3  0.669 0.677                                          
4  0.804 0.812 0.792                                    
5  0.632 0.637 0.578 0.737                              
6  0.647 0.653 0.591 0.763 0.558                        
7  0.691 0.697 0.652 0.788 0.617 0.631                  
8  0.646 0.653 0.590 0.763 0.558 0.568 0.631            
9  0.728 0.737 0.695 0.838 0.651 0.669 0.713 0.668      
10 0.635 0.641 0.583 0.740 0.554 0.563 0.621 0.563 0.654
> all.equal(genet.dist(d2,method="WC84"),genet.dist(d3,method="WC84"))
[1] TRUE

If you are happy with the fix, could you please close the issue? Thanks.

silknets commented 1 year ago

I re-installed the package from your GitHub and I was able to run the pairwise.WCfst(recode) script successfully!

Much appreciated,

Sam Silknetter

On Mon, Mar 27, 2023 at 4:13 PM Jerome Goudet @.***> wrote:

@silknets https://github.com/silknets ,

The issue should be resolved (see code https://github.com/jgx65/hierfstat/commit/c8d5ea935c40b6a9a12efb17de169aba4c2c9a47) . Let me know if this is the case. Cheers

— Reply to this email directly, view it on GitHub https://github.com/jgx65/hierfstat/issues/66#issuecomment-1485803431, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALH7V5F37X3IWILHRF64IUDW6HYGTANCNFSM6AAAAAAWEAXWWE . You are receiving this because you were mentioned.Message ID: @.***>