I am presently working with the PopGenome package in R to calculate the haplotype diversity in a large list of exons in immune genes. All went well for around 550 exons. However, the analysis has stopped at one of the exons, giving the following error message:
Error in tmp[[xx]] : subscript out of bounds Calls: diversity.stats -> diversity.stats Execution halted
The code is as follows:
library(PopGenome)
AFR <- readLines("/home/.../AFR.txt")
AMR <- readLines("/home/.../AMR.txt")
EAS <- readLines("/home/.../EAS.txt")
EUR <- readLines("/home/.../EUR.txt")
SAS <- readLines("/home/.../SAS.txt")
GENOME.class <- readVCF("/data/.../ALL.chr19.phase3_shapeit2_mvncall_integrated_v5.20130502.genotypes.vcf.gz", 10000,"19", 54701102, 54701178, include.unknown=FALSE, approx=FALSE, parallel=FALSE)
GENOME.class.pop <- set.populations(GENOME.class,list(AFR,AMR,EAS,EUR,SAS), diploid=TRUE)
div.stat <- diversity.stats(GENOME.class.pop)
hap.div <- div.stat@hap.diversity.within
outputfile <-write.csv(hap.div, "/home/.../1K_KIR2DL1_1_HaplotypeDiversity.csv", row.names = FALSE)
I am presently working with the PopGenome package in R to calculate the haplotype diversity in a large list of exons in immune genes. All went well for around 550 exons. However, the analysis has stopped at one of the exons, giving the following error message:
Error in tmp[[xx]] : subscript out of bounds Calls: diversity.stats -> diversity.stats Execution halted
The code is as follows:
library(PopGenome)
AFR <- readLines("/home/.../AFR.txt")
AMR <- readLines("/home/.../AMR.txt")
EAS <- readLines("/home/.../EAS.txt")
EUR <- readLines("/home/.../EUR.txt")
SAS <- readLines("/home/.../SAS.txt")
GENOME.class <- readVCF("/data/.../ALL.chr19.phase3_shapeit2_mvncall_integrated_v5.20130502.genotypes.vcf.gz", 10000,"19", 54701102, 54701178, include.unknown=FALSE, approx=FALSE, parallel=FALSE)
GENOME.class.pop <- set.populations(GENOME.class,list(AFR,AMR,EAS,EUR,SAS), diploid=TRUE)
div.stat <- diversity.stats(GENOME.class.pop)
hap.div <- div.stat@hap.diversity.within
outputfile <-write.csv(hap.div, "/home/.../1K_KIR2DL1_1_HaplotypeDiversity.csv", row.names = FALSE)
I tried to do a google search and figure out what is wrong, but I am unable to find out how to solve the problem.
The code for another exon that worked fine is as follows:
library(PopGenome)
AFR <- readLines("/home/.../AFR.txt")
AMR <- readLines("/home/.../AMR.txt")
EAS <- readLines("/home/.../EAS.txt")
EUR <- readLines("/home/.../EUR.txt")
SAS <- readLines("/home/.../SAS.txt")
GENOME.class <- readVCF("/data/.../ALL.chr19.phase3_shapeit2_mvncall_integrated_v5.20130502.genotypes.vcf.gz", 10000,"19", 54701102, 54701178, include.unknown=FALSE, approx=FALSE, parallel=FALSE)
GENOME.class.pop <- set.populations(GENOME.class,list(AFR,AMR,EAS,EUR,SAS), diploid=TRUE)
div.stat <- diversity.stats(GENOME.class.pop)
hap.div <- div.stat@hap.diversity.within
outputfile <-write.csv(hap.div, "/home/.../1K_KIR2DL1_1_HaplotypeDiversity.csv", row.names = FALSE)
Any help would be most appreciated. Thank you!!
