Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Error: genind2hierfstat conversion #75

Open
oninlaurel opened this issue Nov 14, 2024 · 3 comments
Open

Error: genind2hierfstat conversion #75

oninlaurel opened this issue Nov 14, 2024 · 3 comments

Comments

@oninlaurel
Copy link

Hi! I have been trying to convert my genind file to hierfstat so I can calculate some population statistics but I keep on getting an error when I run my code. Here's the error:

Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE,  : 
  arguments imply differing number of rows: 189, 193, 190, 196, 195, 194, 191, 192, 187, 188

For reference, I am using an LD-pruned SNP (r2=0.1) dataset obtained using plink. Before LD-pruning, I only included SNPs with F_MISSING=0.1. I am not sure what's happening with the conversion to hierfstat. I would appreciate if you could help me with this matter. Thank you!

library("poppr")
library("ape")
library("magrittr")
library("adegenet")
library("reshape")
library('vcfR')
library('parallel')
library('hierfstat')

vcf<- read.vcfR("Sz_13540_pruned0.1.vcf") #Pruned SNPs (r2=0.1) 98 strains
Sz_grp<-read.delim("Sz_popdata_clade.txt", header = TRUE, sep = "\t", fill = TRUE, strip.white = TRUE, check.names = TRUE)

#check that all sample names are retained and in the right order
setdiff(colnames(vcf@gt)[-1], Sz_grp$name)# character(0) = TRUE - all the names are correct and all the samples are in there
all(colnames(vcf@gt)[-1] == Sz_grp$name) # FALSE - samples are not in the correct order

#put samples in the right order (to match vcf order)
to_be_ordered_strata<- as.data.frame(cbind(Sz_grp$name, Sz_grp$clade))
sort_on<- as.data.frame(colnames(vcf@gt)[-1])
ordered_strata<- to_be_ordered_strata[order(match(to_be_ordered_strata[,1],sort_on[,1])),]

#check that it's in order now
all(colnames(vcf@gt)[-1] == ordered_strata$V1) # TRUE - pop now ordered to match vcf. 

##Convert to genID object
genind_ob<- vcfR2genind(vcf, ploidy = 1)

#set pop as strata
colnames(ordered_strata)<- c("strain", "pop")
strata(genind_ob) #NULL
strata(genind_ob)<-ordered_strata
strata(genind_ob) # now set to clade

#check pop specifically 
head(pop(genind_ob)) #nope
setPop(genind_ob) <- ~pop
head(pop(genind_ob)) #now set to the pop strata

#clean up large vcf to free memory
rm(vcf)

#create hierfstat object
my_genind2 <- genind2hierfstat(genind_ob)
@jgx65
Copy link
Owner

jgx65 commented Nov 14, 2024 •
edited
Loading

Hi @oninlaurel , could you share the result of str(genind_ob) ? Thanks. (note that many functions in hierfstat accept genind objects as arguments

@oninlaurel
Copy link
Author

Hi @jgx65! Thanks for your response. Here's the result of str(genind_ob)

Formal class 'genind' [package "adegenet"] with 11 slots
  ..@ tab      : int [1:98, 1:65552] NA NA 1 1 1 1 1 0 1 1 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:98] "Sz10" "Sz104" "Sz107" "Sz108" ...
  .. .. ..$ : chr [1:65552] "scaffold_1:583:A:G.0/0" "scaffold_1:583:A:G.1/1" "scaffold_1:607:G:A.0/0" "scaffold_1:607:G:A.1/1" ...
  ..@ loc.fac  : Factor w/ 32776 levels "scaffold_1:583:A:G",..: 1 1 2 2 3 3 4 4 5 5 ...
  ..@ loc.n.all: Named int [1:32776] 2 2 2 2 2 2 2 2 2 2 ...
  .. ..- attr(*, "names")= chr [1:32776] "scaffold_1:583:A:G" "scaffold_1:607:G:A" "scaffold_1:624:C:T" "scaffold_1:630:C:T" ...
  ..@ all.names:List of 32776
  .. ..$ scaffold_1:583:A:G     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:607:G:A     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:624:C:T     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:630:C:T     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:643:C:T     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:668:C:T     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:683:G:A     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:691:C:T     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:832:C:T     : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:868:C:T     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:871:G:A     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:881:G:A     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:886:A:G     : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:893:C:T     : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:899:T:C     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:905:G:A     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:912:G:A     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:925:T:C     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:945:G:A     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:990:G:A     : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:1013:A:G    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:1019:G:A    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:1025:C:T    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:1191:T:A    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:2563:T:C    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:2705:G:A    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:2745:G:A    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:2904:G:A    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:3086:G:A    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:3145:G:A    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:3200:C:T    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:3807:G:A    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:4224:G:A    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:4525:T:C    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:4602:G:A    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:5002:A:G    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:5093:G:A    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:5495:A:T    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:5987:T:G    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:6135:G:A    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:6712:A:G    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:6776:A:G    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:7000:A:T    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:7615:A:G    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:7945:G:A    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:7983:A:G    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:8067:G:A    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:8184:T:C    : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:8627:C:G    : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:10324:A:C   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:10829:C:T   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:13510:C:T   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:13967:A:T   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:14286:A:G   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:19881:G:A   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:20565:C:T   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:21617:G:C   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:21746:G:T   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:22043:C:T   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:22878:G:C   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:23801:G:A   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:23931:C:T   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:24876:G:A   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:25583:T:G   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:25830:A:G   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:26174:C:A   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:29320:A:G   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:29876:T:C   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:33982:T:C   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:35973:A:G   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:36324:C:G   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:36780:T:C   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:36912:T:G   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:36928:G:A   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:37331:C:T   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:37722:G:A   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:38487:A:G   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:38976:C:G   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:39204:A:T   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:39575:T:C   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:39640:C:A   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:40316:C:A   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:40990:C:T   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:41314:G:A   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:42796:G:A   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:44282:G:T   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:44926:A:G   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:45775:G:A   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:47421:A:G   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:48637:T:G   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:51487:T:G   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:51762:C:G   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:51777:T:C   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:52397:C:A   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:52433:T:C   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:52511:C:G   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:52528:C:T   : chr [1:2] "0/0" "1/1"
  .. ..$ scaffold_1:52679:T:G   : chr [1:2] "1/1" "0/0"
  .. ..$ scaffold_1:53278:C:T   : chr [1:2] "0/0" "1/1"
  .. .. [list output truncated]
  ..@ ploidy   : int [1:98] 1 1 1 1 1 1 1 1 1 1 ...
  ..@ type     : chr "codom"
  ..@ other    : NULL
  ..@ call     : language adegenet::df2genind(X = t(x), sep = sep, ploidy = 1)
  ..@ pop      : NULL
  ..@ strata   : NULL
  ..@ hierarchy: NULL

@jgx65
Copy link
Owner

jgx65 commented Nov 14, 2024

This seems correct. I see @pop is NULL, didn't you set it? What is best is to post a small example (dataset+commands) reproducing the error. Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@jgx65 @oninlaurel