simuations with sim1000g
0
0
Entering edit mode
10 months ago
Eliza ▴ 30

Hi, I'm trying to simulate the 22 chromosomes with sim1000G package in R using the for loop the problem is that for the first run of the loop, the code works well but after it cant produce the genetics map.

This is my code

for (i in c(1:2)) {
  print("#####")
  print(i)
  print("######")
  vcf_file = file.path("chr_22_final2.vcf")
  vcf = readVCF( vcf_file , maxNumberOfVariants = 2000 , min_maf = 0.02 , max_maf = NA )
  startSimulation(vcf, totalNumberOfIndividuals = 1200)
  ids = generateUnrelatedIndividuals(1000)
  genotype = retrieveGenotypes(ids)
  E_vector=genotype[,i+300]
  #E_vector[E_vector==2]=1
  #genotype[genotype==2]=1
  result_matrix=genotype[,-(i+300)]
  data_after_sim_funce=data.gen(sample_size=900, p=1999, 
                                                     n_g_non_zero=200, n_gxe_non_zero=20, 
                                                     family="binomial", mode="strong_hierarchical")
  tune_model = gesso.cv(data_after_sim_funce$G_train, data_after_sim_funce$E_train, data_after_sim_funce$Y_train, 
                      grid_size=20, tolerance=1e-4,
                      parallel=TRUE, nfold=4,
                      normalize=FALSE, normalize_response=FALSE,
                      seed=1)
gxe_coefficients = gesso.coef(tune_model$fit, tune_model$lambda_min)$beta_gxe
g_coefficients = gesso.coef(tune_model$fit, tune_model$lambda_min)$beta_g
print(length(g_coefficients[g_coefficients>0]))
print(length(gxe_coefficients[gxe_coefficients>0]))
num_main=length(g_coefficients[g_coefficients>0]) 
num_iter=length(gxe_coefficients[gxe_coefficients>0])
num_main_large_zero[i]=num_main
num_iter_large_zero[i]=num_iter
}

for index i =1 it works correctly but when the loop moves to i=2 it returns an error : [#####...] Creating SIM object [#####...] Haplodata object created Error in startSimulation(vcf, totalNumberOfIndividuals = 1200) : Error: Genetic map does not match the region being simulated

I can't seem to understand what is happening because the code worked perfectly well for i=1 but crashes on i=2. I would be happy for your opinion :) should i delete after every iteration the "old" SIM object ?

snps simulation sim1000G chromosome • 332 views
ADD COMMENT

Login before adding your answer.

Traffic: 3900 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6