diff --git a/src/base/vcf.rs b/src/base/vcf.rs index d3df1e3..eda1619 100644 --- a/src/base/vcf.rs +++ b/src/base/vcf.rs @@ -141,7 +141,8 @@ impl Filter for VcfLine { while j < m { q = 0.0; for i in 0..n { - q += allele_frequencies.matrix[(i, j)] * filter_stats.pool_sizes[i]; + q += allele_frequencies.matrix[(i, j)] * filter_stats.pool_sizes[i] + / filter_stats.pool_sizes.iter().sum::(); // We've made sure the pool_sizes sum up to one in phen.rs } if (q < filter_stats.min_allele_frequency) diff --git a/src/imputation/adaptive_ld_knn_imputation.rs b/src/imputation/adaptive_ld_knn_imputation.rs index 0f7a49f..06f2608 100644 --- a/src/imputation/adaptive_ld_knn_imputation.rs +++ b/src/imputation/adaptive_ld_knn_imputation.rs @@ -319,7 +319,7 @@ impl GenotypesAndPhenotypes { } // Impute across pools with missing data } // Impute if we have missing data } // Iterate across alleles across loci within the window - // println!("window_freqs:\n{:?}", &window_freqs); + // println!("window_freqs:\n{:?}", &window_freqs); }, // Parallel processing across windows ); // } diff --git a/src/main.rs b/src/main.rs index 5b473ea..3e58d75 100644 --- a/src/main.rs +++ b/src/main.rs @@ -8,8 +8,8 @@ mod base; mod gp; mod gwas; mod imputation; -mod simulation; mod popgen; +mod simulation; mod tables; use base::{ChunkyReadAnalyseWrite, CrossValidation, LoadAll, Parse, SaveCsv}; use gp::{ @@ -18,8 +18,8 @@ use gp::{ }; use gwas::*; use imputation::*; -use simulation::*; use popgen::*; +use simulation::*; #[derive(Parser, Debug)] #[clap( diff --git a/src/simulation/simulate_genotypes.rs b/src/simulation/simulate_genotypes.rs index 1a775e4..a6d4234 100644 --- a/src/simulation/simulate_genotypes.rs +++ b/src/simulation/simulate_genotypes.rs @@ -1,23 +1,29 @@ use crate::base::*; use ndarray::prelude::*; -use std::io; -use statrs::distribution::{Uniform, MultivariateNormal, Continuous}; use rand::distributions::Distribution; use rand::{seq::IteratorRandom, thread_rng}; +use statrs::distribution::{Continuous, MultivariateNormal, Uniform}; +use std::io; // Simulate genotypes with some linkage disequilibrium -pub fn simulate_genotypes(n: usize, p: usize, n_chr: usize, max_bp: usize, r2_50_perc_bp: usize) -> io::Result> { +pub fn simulate_genotypes( + n: usize, + p: usize, + n_chr: usize, + max_bp: usize, + r2_50_perc_bp: usize, +) -> io::Result> { // Define approximately equally sized chromosomes and coordinates of the characterised loci or markers or SNPs - let mut chromosome_sizes = vec![p/n_chr; n_chr]; + let mut chromosome_sizes = vec![p / n_chr; n_chr]; let p_ = chromosome_sizes.iter().fold(0, |sum, &x| sum + x); // If the we have less or more than the required number of loci we add or subtract loci from the last chromosome if p_ < p { - chromosome_sizes[n_chr-1] += p - p_; + chromosome_sizes[n_chr - 1] += p - p_; } else if p_ > p { - chromosome_sizes[n_chr-1] -= p_ - p; + chromosome_sizes[n_chr - 1] -= p_ - p; } - let max_bp_per_chr = max_bp/n_chr; + let max_bp_per_chr = max_bp / n_chr; println!("Chromosome sizes: {:?}", chromosome_sizes); // Sample uniformly distributed loci per chromosome let dist_unif = Uniform::new(0.0, max_bp_per_chr as f64).unwrap(); @@ -27,7 +33,10 @@ pub fn simulate_genotypes(n: usize, p: usize, n_chr: usize, max_bp: usize, r2_50 let mut chromosomes: Array1 = Array1::from_elem(p, "".to_owned()); let mut positions: Array1 = Array1::from_elem(p, 0); let mut alleles: Array1 = Array1::from_elem(p, "".to_owned()); - let atcgd = vec!["A", "T", "C", "G", "D"].iter().map(|&x| x.to_owned()).collect::>(); + let atcgd = vec!["A", "T", "C", "G", "D"] + .iter() + .map(|&x| x.to_owned()) + .collect::>(); for i in 0..n_chr { // Total number of sites we want to sample in the chromosome let s = chromosome_sizes[i]; @@ -50,8 +59,7 @@ pub fn simulate_genotypes(n: usize, p: usize, n_chr: usize, max_bp: usize, r2_50 let mvn = MultivariateNormal::new(vec![0., 0.], vec![1., 0., 0., 1.]).unwrap(); - - return Ok(Array2::from_elem((1, 1), f64::NAN)) + return Ok(Array2::from_elem((1, 1), f64::NAN)); } ////////////////////////////////////////////////////////////////////////////////////////////////////////////