From 70528756fc554001578114ffcf976a98dbb7622a Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Tue, 25 Feb 2020 11:52:41 +0100 Subject: [PATCH] Added computeAllFreqs2 function --- R/admix1.R | 56 +++++++++++------------------------------ R/computeAllFreqs2.R | 23 +++++++++++++++++ man/computeAllFreqs2.Rd | 12 +++++++++ 3 files changed, 50 insertions(+), 41 deletions(-) create mode 100644 R/computeAllFreqs2.R create mode 100644 man/computeAllFreqs2.Rd diff --git a/R/admix1.R b/R/admix1.R index 28ce731..c50d4fe 100644 --- a/R/admix1.R +++ b/R/admix1.R @@ -5,7 +5,7 @@ admix1 <- function(tietue) { if (!is.list(tietue)) { # c(filename, pathname) = uigetfile('*.mat', 'Load mixture result file'); -# if (filename==0 & pathname==0), return; +# if (filename==0 & pathname==0), return; # else # disp('---------------------------------------------------'); # disp(['Reading mixture result from: ',[pathname filename],'...']); @@ -13,7 +13,7 @@ admix1 <- function(tietue) { # pause(0.0001); # h0 = findobj('Tag','filename1_text'); # set(h0,'String',filename); clear h0; - + # struct_array = load([pathname filename]); # if isfield(struct_array,'c') #Matlab versio # c = struct_array.c; @@ -31,7 +31,7 @@ admix1 <- function(tietue) { # disp('Incorrect file format'); # return; # end - + # if isfield(c, 'gene_lengths') && ... # (strcmp(c.mixtureType,'linear_mix') | ... # strcmp(c.mixtureType,'codon_mix')) # if the mixture is from a linkage model @@ -40,7 +40,7 @@ admix1 <- function(tietue) { # linkage_admix(c); # return # end - + # PARTITION = c.PARTITION; COUNTS = c.COUNTS; SUMCOUNTS = c.SUMCOUNTS; # alleleCodes = c.alleleCodes; adjprior = c.adjprior; popnames = c.popnames; # rowsFromInd = c.rowsFromInd; data = c.data; npops = c.npops; noalle = c.noalle; @@ -148,7 +148,7 @@ admix1 <- function(tietue) { # for iterationNum = 1:iterationCount # disp(['Iter: ' num2str(iterationNum)]); # allfreqs = simulateAllFreqs(noalle); # Allele frequencies on this iteration. - + # for ind=to_investigate # #disp(num2str(ind)); # omaFreqs = computePersonalAllFreqs(ind, data, allfreqs, rowsFromInd); @@ -172,7 +172,7 @@ admix1 <- function(tietue) { # PARTITION(ind)=isoimman_indeksi; # end # logml = computeIndLogml(omaFreqs, osuusTaulu); - + # for osuus = [0.5 0.25 0.05 0.01] # [osuusTaulu, logml] = etsiParas(osuus, osuusTaulu, omaFreqs, logml); # end @@ -227,23 +227,23 @@ admix1 <- function(tietue) { # for pop = admix_populaatiot' # for level = 1:n_missing_levels(pop) - + # potential_inds_in_this_pop_and_level = ... # find(PARTITION==pop & missing_level_partition==level &... # likelihood>3); # Potential admix individuals here. - + # if ~isempty(potential_inds_in_this_pop_and_level) - + # #refData = simulateIndividuals(nrefIndsInPop,rowsFromInd,allfreqs); # refData = simulateIndividuals(nrefIndsInPop, rowsFromInd, allfreqs, ... # pop, missing_levels(pop,level)); - + # disp(['Analysing the reference individuals from pop ' num2str(pop) ' (level ' num2str(level) ').']); # refProportions = zeros(nrefIndsInPop,npops); # for iter = 1:iterationCountRef # #disp(['Iter: ' num2str(iter)]); # allfreqs = simulateAllFreqs(noalle); - + # for ind = 1:nrefIndsInPop # omaFreqs = computePersonalAllFreqs(ind, refData, allfreqs, rowsFromInd); # osuusTaulu = zeros(1,npops); @@ -312,7 +312,7 @@ admix1 <- function(tietue) { # end # end -# tulostaAdmixtureTiedot(proportionsIt, uskottavuus, alaRaja, iterationCount); +# tulostaAdmixtureTiedot(proportionsIt, uskottavuus, alaRaja, iterationCount); # viewPartition(proportionsIt, popnames); @@ -335,7 +335,7 @@ admix1 <- function(tietue) { # if (~isstruct(tietue)) -# c.proportionsIt = proportionsIt; +# c.proportionsIt = proportionsIt; # c.pvalue = uskottavuus; # Added by Jing # c.mixtureType = 'admix'; # Jing # c.admixnpops = npops; @@ -355,7 +355,7 @@ admix1 <- function(tietue) { # function [npops] = poistaLiianPienet(npops, rowsFromInd, alaraja) # % Muokkaa tulokset muotoon, jossa outlier yksilöt on -# % poistettu. Tarkalleen ottaen poistaa ne populaatiot, +# % poistettu. Tarkalleen ottaen poistaa ne populaatiot, # % joissa on vähemmän kuin 'alaraja':n verran yksilöit? # global PARTITION; @@ -400,30 +400,4 @@ admix1 <- function(tietue) { # global COUNTS; COUNTS = []; # global SUMCOUNTS; SUMCOUNTS = []; # global PARTITION; PARTITION = []; -# global POP_LOGML; POP_LOGML = []; - -# %-------------------------------------------------------- - - -# function allFreqs = computeAllFreqs2(noalle) -# % Lisää a priori jokaista alleelia -# % joka populaation joka lokukseen j 1/noalle(j) verran. - -# global COUNTS; -# global SUMCOUNTS; - -# max_noalle = size(COUNTS,1); -# nloci = size(COUNTS,2); -# npops = size(COUNTS,3); - -# sumCounts = SUMCOUNTS+ones(size(SUMCOUNTS)); -# sumCounts = reshape(sumCounts', [1, nloci, npops]); -# sumCounts = repmat(sumCounts, [max_noalle, 1 1]); - -# prioriAlleelit = zeros(max_noalle,nloci); -# for j=1:nloci -# prioriAlleelit(1:noalle(j),j) = 1/noalle(j); -# end -# prioriAlleelit = repmat(prioriAlleelit, [1,1,npops]); -# counts = COUNTS + prioriAlleelit; -# allFreqs = counts./sumCounts; \ No newline at end of file +# global POP_LOGML; POP_LOGML = []; \ No newline at end of file diff --git a/R/computeAllFreqs2.R b/R/computeAllFreqs2.R new file mode 100644 index 0000000..290c41c --- /dev/null +++ b/R/computeAllFreqs2.R @@ -0,0 +1,23 @@ +#' @title Compute all freqs - version 2 +#' @description Lisää a priori jokaista alleelia joka populaation joka lokukseen +#' j 1/noalle(j) verran. +computeAllFreqs2 <- function (noalle, COUNTS = matrix(NA, 0, 0), + SUMCOUNTS = sum(COUNTS)) { + + max_noalle <- size(COUNTS, 1) + nloci <- size(COUNTS,2) + npops <- size(COUNTS,3) + + sumCounts <- SUMCOUNTS + ones(size(SUMCOUNTS)) + sumCounts <- reshape(t(sumCounts), c(1, nloci, npops)) + sumCounts <- repmat(sumCounts, c(max_noalle, 1, 1)) + + prioriAlleelit <- zeros(max_noalle, nloci) + for (j in 1:nloci) { + prioriAlleelit[1:noalle[j], j] <- 1 / noalle[j] + } + prioriAlleelit <- repmat(prioriAlleelit, c(1, 1, npops)) + counts <- COUNTS + prioriAlleelit + allFreqs <- counts / sumCounts + return(allFreqs) +} \ No newline at end of file diff --git a/man/computeAllFreqs2.Rd b/man/computeAllFreqs2.Rd new file mode 100644 index 0000000..fcb125b --- /dev/null +++ b/man/computeAllFreqs2.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/computeAllFreqs2.R +\name{computeAllFreqs2} +\alias{computeAllFreqs2} +\title{Compute all freqs - version 2} +\usage{ +computeAllFreqs2(noalle, COUNTS = matrix(NA, 0, 0), SUMCOUNTS = sum(COUNTS)) +} +\description{ +Lisää a priori jokaista alleelia joka populaation joka lokukseen +j 1/noalle(j) verran. +}