diff --git a/NAMESPACE b/NAMESPACE index afc1b39..816fc08 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,6 +8,7 @@ export(computeIndLogml) export(computePersonalAllFreqs) export(computeRows) export(etsiParas) +export(inputdlg) export(isfield) export(laskeMuutokset4) export(learn_simple_partition) diff --git a/R/admix1.R b/R/admix1.R index e678b49..9cf5bb1 100644 --- a/R/admix1.R +++ b/R/admix1.R @@ -1,49 +1,64 @@ #' @title Admixture analysis -#' @param tietue record -#' @details If the record == -1, the mixture results file is loaded. Otherwise, will the required variables be retrieved from the record fields? +#' @param tietue a named record list +#' @details If the record == -1, the mixture results file is loaded. Otherwise, +#' will the required variables be retrieved from the record fields? +#' `tietue`should contain the following elements: PARTITION, COUNTS, SUMCOUNTS, +#' alleleCodes, adjprior, popnames, rowsFromInd, data, npops, noalle #' @export -admix1 <- function(tietue) { +admix1 <- function(tietue, PARTITION = matrix(NA, 0, 0), + COUNTS = matrix(NA, 0, 0), SUMCOUNTS = NA) { if (!is.list(tietue)) { -# c(filename, pathname) = uigetfile('*.mat', 'Load mixture result file'); -# if (filename==0 & pathname==0), return; -# else -# disp('---------------------------------------------------'); -# disp(['Reading mixture result from: ',[pathname filename],'...']); -# end -# pause(0.0001); -# h0 = findobj('Tag','filename1_text'); -# set(h0,'String',filename); clear h0; + message('Load mixture result file. These are the files in this directory:') + print(list.files()) + pathname_filename <- file.choose() + if (!file.exists(pathname_filename)) { + stop( + "File ", pathname_filename, + " does not exist. Check spelling and location." + ) + } else { + cat('---------------------------------------------------\n'); + message('Reading mixture result from: ', pathname_filename, '...') + } + sys.sleep(0.0001) #ASK: what for? -# struct_array = load([pathname filename]); -# if isfield(struct_array,'c') #Matlab versio -# c = struct_array.c; -# if ~isfield(c,'PARTITION') | ~isfield(c,'rowsFromInd') -# disp('Incorrect file format'); -# return -# end -# elseif isfield(struct_array,'PARTITION') #Mideva versio -# c = struct_array; -# if ~isfield(c,'rowsFromInd') -# disp('Incorrect file format'); -# return -# end -# else -# disp('Incorrect file format'); -# return; -# end + # ASK: what is this supposed to do? What do graphic obj have to do here? + # h0 = findobj('Tag','filename1_text'); + # set(h0,'String',filename); clear h0; -# if isfield(c, 'gene_lengths') && ... -# (strcmp(c.mixtureType,'linear_mix') | ... -# strcmp(c.mixtureType,'codon_mix')) # if the mixture is from a linkage model -# # Redirect the call to the linkage admixture function. -# c.data = noIndex(c.data,c.noalle); # call function noindex to remove the index column -# linkage_admix(c); -# return -# end + struct_array <- load(pathname_filename) + if (isfield(struct_array, 'c')) { #Matlab versio + c <- struct_array$c + if (!isfield(c, 'PARTITION') | !isfield(c,'rowsFromInd')) { + stop('Incorrect file format') + } + } else if (isfield(struct_array, 'PARTITION')) { #Mideva versio + c <- struct_array + if (!isfield(c,'rowsFromInd')) stop('Incorrect file format') + } else { + stop('Incorrect file format') + } -# 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; + if (isfield(c, 'gene_lengths') & + strcmp(c$mixtureType, 'linear_mix') | + strcmp(c$mixtureType, 'codon_mix')) { # if the mixture is from a linkage model + # Redirect the call to the linkage admixture function. + # call function noindex to remove the index column + c$data <- noIndex(c$data, c$noalle) + # linkage_admix(c) # ASK: translate this function to R or drop? + # return + stop("linkage_admix not implemented") + } + 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 } else { PARTITION <- tietue$PARTITION COUNTS <- tietue$COUNTS @@ -57,297 +72,334 @@ admix1 <- function(tietue) { noalle <- tietue$noalle } -# answers = inputdlg({['Input the minimum size of a population that will'... -# ' be taken into account when admixture is estimated.']},... -# 'Input minimum population size',[1],... -# {'5'}); -# if isempty(answers) return; end -# alaRaja = str2num(answers{1,1}); -# [npops] = poistaLiianPienet(npops, rowsFromInd, alaRaja); + answers <- inputdlg( + prompt = paste( + "Input the minimum size of a population that will", + "be taken into account when admixture is estimated." + ), + definput = 5 + ) + alaRaja <- as.num(answers) + npops <- poistaLiianPienet(npops, rowsFromInd, alaRaja) -# nloci = size(COUNTS,2); -# ninds = size(data,1)/rowsFromInd; + nloci <- size(COUNTS, 2) + ninds <- size(data, 1) / rowsFromInd -# answers = inputdlg({['Input number of iterations']},'Input',[1],{'50'}); -# if isempty(answers) return; end -# iterationCount = str2num(answers{1,1}); + answers <- inputdlg('Input number of iterations', 50) + if (isempty(answers)) return() + iterationCount <- as.numeric(answers[1, 1]) # maybe [[]]? -# answers = inputdlg({['Input number of reference individuals from each population']},'Input',[1],{'50'}); -# if isempty(answers) nrefIndsInPop = 50; -# else nrefIndsInPop = str2num(answers{1,1}); -# end + answers <- inputdlg( + prompt = 'Input number of reference individuals from each population', + definput = 50 + ) + if (isempty(answers)) { + nrefIndsInPop <- 50 + } else { + nrefIndsInPop <- as.numeric(answers[1, 1]) + } -# answers = inputdlg({['Input number of iterations for reference individuals']},'Input',[1],{'10'}); -# if isempty(answers) return; end -# iterationCountRef = str2num(answers{1,1}); + answers <- inputdlg( + prompt = 'Input number of iterations for reference individuals', + definput = 10 + ) + if (isempty(answers)) return() + iterationCountRef <- as.numeric(answers[1, 1]) + + # First calculate log-likelihood ratio for all individuals: + likelihood <- zeros(ninds, 1) + allfreqs <- computeAllFreqs2(noalle) + for (ind in 1:ninds) { + omaFreqs <- computePersonalAllFreqs(ind, data, allfreqs, rowsFromInd) + osuusTaulu <- zeros(1, npops) + if (PARTITION[ind] == 0) { + # Yksil?on outlier + } else if (PARTITION[ind] != 0) { + if (PARTITION[ind] > 0) { + osuusTaulu(PARTITION[ind]) <- 1 + } else { + # Yksilöt, joita ei ole sijoitettu mihinkään koriin. + arvot <- zeros(1, npops) + for (q in 1:npops) { + osuusTaulu <- zeros(1, npops) + osuusTaulu[q] <- 1 + arvot[q] <- computeIndLogml(omaFreqs, osuusTaulu) + } + iso_arvo <- max(arvot) + isoimman_indeksi <- match(max(arvot), arvot) + osuusTaulu <- zeros(1, npops) + osuusTaulu[isoimman_indeksi] <- 1 + PARTITION[ind] <- isoimman_indeksi + } + logml <- computeIndLogml(omaFreqs, osuusTaulu) + logmlAlku <- logml + for (osuus in c(0.5, 0.25, 0.05, 0.01)) { + etsiResult <- etsiParas(osuus, osuusTaulu, omaFreqs, logml) + osuusTaulu <- etsiResult[1] + logml <- etsiResult[2] + } + logmlLoppu <- logml + likelihood[ind] <- logmlLoppu - logmlAlku + } + } + + # Analyze further only individuals who have log-likelihood ratio larger than 3: + to_investigate <- t(find(likelihood > 3)) + cat('Possibly admixed individuals:\n') + for (i in 1:length(to_investigate)) { + cat(as.character(to_investigate[i])) + } + cat(' ') + cat('Populations for possibly admixed individuals:\n') + admix_populaatiot <- unique(PARTITION[to_investigate]) + for (i in 1:length(admix_populaatiot)) { + cat(as.character(admix_populaatiot[i])) + } + + # THUS, there are two types of individuals, who will not be analyzed with + # simulated allele frequencies: those who belonged to a mini-population + # which was removed, and those who have log-likelihood ratio less than 3. + # The value in the PARTITION for the first kind of individuals is 0. The + # second kind of individuals can be identified, because they do not + # belong to "to_investigate" array. When the results are presented, the + # first kind of individuals are omitted completely, while the second kind + # of individuals are completely put to the population, where they ended up + # in the mixture analysis. These second type of individuals will have a + # unit p-value. -# # First calculate log-likelihood ratio for all individuals: -# likelihood = zeros(ninds,1); -# allfreqs = computeAllFreqs2(noalle); -# for ind = 1:ninds -# omaFreqs = computePersonalAllFreqs(ind, data, allfreqs, rowsFromInd); -# osuusTaulu = zeros(1,npops); -# if PARTITION(ind)==0 -# # Yksil?on outlier -# elseif PARTITION(ind)~=0 -# if PARTITION(ind)>0 -# osuusTaulu(PARTITION(ind)) = 1; -# else -# # Yksilöt, joita ei ole sijoitettu mihinkään koriin. -# arvot = zeros(1,npops); -# for q=1:npops -# osuusTaulu = zeros(1,npops); -# osuusTaulu(q) = 1; -# arvot(q) = computeIndLogml(omaFreqs, osuusTaulu); -# end -# [iso_arvo, isoimman_indeksi] = max(arvot); -# osuusTaulu = zeros(1,npops); -# osuusTaulu(isoimman_indeksi) = 1; -# PARTITION(ind)=isoimman_indeksi; -# end -# logml = computeIndLogml(omaFreqs, osuusTaulu); -# logmlAlku = logml; -# for osuus = [0.5 0.25 0.05 0.01] -# [osuusTaulu, logml] = etsiParas(osuus, osuusTaulu, omaFreqs, logml); -# end -# logmlLoppu = logml; -# likelihood(ind) = logmlLoppu-logmlAlku; -# end -# end + # Simulate allele frequencies a given number of times and save the average + # result to "proportionsIt" array. -# # Analyze further only individuals who have log-likelihood ratio larger than 3: -# to_investigate = (find(likelihood>3))'; -# disp('Possibly admixed individuals: '); -# for i = 1:length(to_investigate) -# disp(num2str(to_investigate(i))); -# end -# disp(' '); -# disp('Populations for possibly admixed individuals: '); -# admix_populaatiot = unique(PARTITION(to_investigate)); -# for i = 1:length(admix_populaatiot) -# disp(num2str(admix_populaatiot(i))); -# end + proportionsIt <- zeros(ninds, npops) + for (iterationNum in 1:iterationCount) { + cat('Iter:', as.character(iterationNum)) + allfreqs <- simulateAllFreqs(noalle) # Allele frequencies on this iteration. -# # THUS, there are two types of individuals, who will not be analyzed with -# # simulated allele frequencies: those who belonged to a mini-population -# # which was removed, and those who have log-likelihood ratio less than 3. -# # The value in the PARTITION for the first kind of individuals is 0. The -# # second kind of individuals can be identified, because they do not -# # belong to "to_investigate" array. When the results are presented, the -# # first kind of individuals are omitted completely, while the second kind -# # of individuals are completely put to the population, where they ended up -# # in the mixture analysis. These second type of individuals will have a -# # unit p-value. + for (ind in to_investigate) { + #disp(num2str(ind)); + omaFreqs <- computePersonalAllFreqs( + ind, data, allfreqs, rowsFromInd + ) + osuusTaulu = zeros(1, npops) + if (PARTITION[ind] == 0) { + # Yksil?on outlier + } else if (PARTITION[ind] != 0) { + if (PARTITION[ind] > 0) { + osuusTaulu(PARTITION[ind]) <- 1 + } else { + # Yksilöt, joita ei ole sijoitettu mihinkään koriin. + arvot <- zeros(1, npops) + for (q in 1:npops) { + osuusTaulu <- zeros(1, npops) + osuusTaulu[q] <- 1 + arvot[q] <- computeIndLogml(omaFreqs, osuusTaulu) + } + iso_arvo <- max(arvot) + isoimman_indeksi <- match(max(arvot), arvot) + osuusTaulu <- zeros(1, npops) + osuusTaulu[isoimman_indeksi] <- 1 + PARTITION[ind] <- isoimman_indeksi + } + logml <- computeIndLogml(omaFreqs, osuusTaulu) + + for (osuus in c(0.5, 0.25, 0.05, 0.01)) { + etsiResult <- etsiParas(osuus, osuusTaulu, omaFreqs, logml) + osuusTaulu <- etsiResult[1] + logml <- etsiResult[2] + } + } + proportionsIt[ind, ] <- proportionsIt[ind, ] * (iterationNum - 1) + + osuusTaulu + proportionsIt[ind, ] <- proportionsIt[ind, ] / iterationNum + } + } + + #disp(['Creating ' num2str(nrefIndsInPop) ' reference individuals from ']); + #disp('each population.'); + + #allfreqs = simulateAllFreqs(noalle); # Simuloidaan alleelifrekvenssisetti + allfreqs <- computeAllFreqs2(noalle); # Koitetaan tällaista. -# # Simulate allele frequencies a given number of times and save the average -# # result to "proportionsIt" array. + # Initialize the data structures, which are required in taking the missing + # data into account: + n_missing_levels <- zeros(npops, 1) # number of different levels of "missingness" in each pop (max 3). + missing_levels <- zeros(npops, 3) # the mean values for different levels. + missing_level_partition <- zeros(ninds, 1) # level of each individual (one of the levels of its population). + for (i in 1:npops) { + inds <- find(PARTITION == i) + # Proportions of non-missing data for the individuals: + non_missing_data <- zeros(length(inds), 1) + for (j in 1:length(inds)) { + ind <- inds[j] + non_missing_data[j] <- length( + find(data[(ind - 1) * rowsFromInd + 1:ind * rowsFromInd, ] > 0) + ) / (rowsFromInd * nloci) + } + if (all(non_missing_data > 0.9)) { + n_missing_levels[i] <- 1 + missing_levels[i, 1] <- mean(non_missing_data) + missing_level_partition[inds] <- 1 + } else { + # TODO: fix syntax + # [ordered, ordering] = sort(non_missing_data); + ordered <- ordering <- sort(non_missing_data) + #part = learn_simple_partition(ordered, 0.05); + part <- learn_partition_modified(ordered) + aux <- sortrows(cbind(part, ordering), 2) + part = aux[, 1] + missing_level_partition[inds]<- part + n_levels <- length(unique(part)) + n_missing_levels[i] <- n_levels + for (j in 1:n_levels) { + missing_levels[i, j] <- mean(non_missing_data[find(part == j)]) + } + } + } -# proportionsIt = zeros(ninds,npops); -# for iterationNum = 1:iterationCount -# disp(['Iter: ' num2str(iterationNum)]); -# allfreqs = simulateAllFreqs(noalle); # Allele frequencies on this iteration. + # Create and analyse reference individuals for populations + # with potentially admixed individuals: + refTaulu <- zeros(npops, 100, 3) + for (pop in t(admix_populaatiot)) { -# for ind=to_investigate -# #disp(num2str(ind)); -# omaFreqs = computePersonalAllFreqs(ind, data, allfreqs, rowsFromInd); -# osuusTaulu = zeros(1,npops); -# if PARTITION(ind)==0 -# # Yksil?on outlier -# elseif PARTITION(ind)~=0 -# if PARTITION(ind)>0 -# osuusTaulu(PARTITION(ind)) = 1; -# else -# # Yksilöt, joita ei ole sijoitettu mihinkään koriin. -# arvot = zeros(1,npops); -# for q=1:npops -# osuusTaulu = zeros(1,npops); -# osuusTaulu(q) = 1; -# arvot(q) = computeIndLogml(omaFreqs, osuusTaulu); -# end -# [iso_arvo, isoimman_indeksi] = max(arvot); -# osuusTaulu = zeros(1,npops); -# osuusTaulu(isoimman_indeksi) = 1; -# PARTITION(ind)=isoimman_indeksi; -# end -# logml = computeIndLogml(omaFreqs, osuusTaulu); + for (level in 1:n_missing_levels[pop]) { -# for osuus = [0.5 0.25 0.05 0.01] -# [osuusTaulu, logml] = etsiParas(osuus, osuusTaulu, omaFreqs, logml); -# end -# end -# proportionsIt(ind,:) = proportionsIt(ind,:).*(iterationNum-1) + osuusTaulu; -# proportionsIt(ind,:) = proportionsIt(ind,:)./iterationNum; -# end -# end + potential_inds_in_this_pop_and_level <- + find( + PARTITION == pop & missing_level_partition == level & + likelihood > 3 + ) # Potential admix individuals here. -# #disp(['Creating ' num2str(nrefIndsInPop) ' reference individuals from ']); -# #disp('each population.'); + if (!isempty(potential_inds_in_this_pop_and_level)) { -# #allfreqs = simulateAllFreqs(noalle); # Simuloidaan alleelifrekvenssisetti -# allfreqs = computeAllFreqs2(noalle); # Koitetaan tällaista. + #refData = simulateIndividuals(nrefIndsInPop,rowsFromInd,allfreqs); + refData <- simulateIndividuals( + nrefIndsInPop, rowsFromInd, allfreqs, pop, + missing_levels[pop, level] + ) + + cat( + 'Analysing the reference individuals from pop', pop, + '(level', level, ').' + ) + refProportions <- zeros(nrefIndsInPop, npops) + for (iter in 1:iterationCountRef) { + #disp(['Iter: ' num2str(iter)]); + allfreqs <- simulateAllFreqs(noalle) + + for (ind in 1:nrefIndsInPop) { + omaFreqs <- computePersonalAllFreqs( + ind, refData, allfreqs, rowsFromInd + ) + osuusTaulu <- zeros(1, npops) + osuusTaulu[pop] <- 1 + logml <- computeIndLogml(omaFreqs, osuusTaulu) + for (osuus in c(0.5, 0.25, 0.05, 0.01)) { + etsiResult <- etsiParas( + osuus, osuusTaulu, omaFreqs, logml + ) + osuusTaulu <- etsiResult[1] + logml <- etsiResult[2] + } + refProportions[ind, ] <- + refProportions[ind, ] * (iter - 1) + osuusTaulu + refProportions[ind, ] <- refProportions[ind, ] / iter + } + } + for (ind in 1:nrefIndsInPop) { + omanOsuus <- refProportions[ind, pop] + if (round(omanOsuus * 100) == 0) { + omanOsuus <- 0.01 + } + if (abs(omanOsuus) < 1e-5) { + omanOsuus <- 0.01 + } + refTaulu[pop, round(omanOsuus*100), level] <- + refTaulu[pop, round(omanOsuus*100),level] + 1 + } + } + } + } + + # Rounding of the results: + proportionsIt <- proportionsIt * 100 + proportionsIt <- round(proportionsIt) + proportionsIt <- proportionsIt / 100 + for (ind in 1:ninds) { + if (!any(to_investigate == ind)) { + if (PARTITION[ind] > 0) { + proportionsIt[ind, PARTITION[ind]] <- 1 + } + } else { + # In case of a rounding error, the sum is made equal to unity by + # fixing the largest value. + if ((PARTITION[ind] > 0) & (sum(proportionsIt[ind, ]) != 1)) { + isoin <- max(proportionsIt[ind, ]) + indeksi <- match(isoin, max(proportionsIt[ind, ])) + erotus <- sum(proportionsIt[ind, ]) - 1 + proportionsIt[ind, indeksi] <- isoin - erotus + } + } + } + + # Calculate p-value for each individual: + uskottavuus <- zeros(ninds, 1) + for (ind in 1:ninds) { + pop <- PARTITION[ind] + if (pop == 0) { # Individual is outlier + uskottavuus[ind] <- 1 + } else if (isempty(find(to_investigate == ind))) { + # Individual had log-likelihood ratio<3 + uskottavuus[ind] <- 1 + } else { + omanOsuus <- proportionsIt[ind, pop] + if (abs(omanOsuus) < 1e-5) { + omanOsuus <- 0.01 + } + if (round(omanOsuus*100)==0) { + omanOsuus <- 0.01 + } + level <- missing_level_partition[ind] + refPienempia <- sum(refTaulu[pop, 1:round(100*omanOsuus), level]) + uskottavuus[ind] <- refPienempia / nrefIndsInPop + } + } + + # ASK: Remove? are these plotting functions? + tulostaAdmixtureTiedot(proportionsIt, uskottavuus, alaRaja, iterationCount) + viewPartition(proportionsIt, popnames) + + talle = inputdlg('Do you want to save the admixture results? [Y/n]', 'y') + if (talle %in% c('y', 'Y', 'yes', 'Yes')) { + #waitALittle; + filename <- inputdlg( + 'Save results as (file name):', 'admixture_results.rda' + ) -# # Initialize the data structures, which are required in taking the missing -# # data into account: -# n_missing_levels = zeros(npops,1); # number of different levels of "missingness" in each pop (max 3). -# missing_levels = zeros(npops,3); # the mean values for different levels. -# missing_level_partition = zeros(ninds,1); # level of each individual (one of the levels of its population). -# for i=1:npops -# inds = find(PARTITION==i); -# # Proportions of non-missing data for the individuals: -# non_missing_data = zeros(length(inds),1); -# for j = 1:length(inds) -# ind = inds(j); -# non_missing_data(j) = length(find(data((ind-1)*rowsFromInd+1:ind*rowsFromInd,:)>0)) ./ (rowsFromInd*nloci); -# end -# if all(non_missing_data>0.9) -# n_missing_levels(i) = 1; -# missing_levels(i,1) = mean(non_missing_data); -# missing_level_partition(inds) = 1; -# else -# [ordered, ordering] = sort(non_missing_data); -# #part = learn_simple_partition(ordered, 0.05); -# part = learn_partition_modified(ordered); -# aux = sortrows([part ordering],2); -# part = aux(:,1); -# missing_level_partition(inds) = part; -# n_levels = length(unique(part)); -# n_missing_levels(i) = n_levels; -# for j=1:n_levels -# missing_levels(i,j) = mean(non_missing_data(find(part==j))); -# end -# end -# end + if (filename == 0) { + # Cancel was pressed + return() + } else { # copy 'baps4_output.baps' into the text file with the same name. + if (file.exists('baps4_output.baps')) { + file.copy('baps4_output.baps', paste0(filename, '.txt')) + file.remove('baps4_output.baps') + } + } -# # Create and analyse reference individuals for populations -# # with potentially admixed individuals: -# refTaulu = zeros(npops,100,3); -# 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); -# osuusTaulu(pop)=1; -# logml = computeIndLogml(omaFreqs, osuusTaulu); -# for osuus = [0.5 0.25 0.05 0.01] -# [osuusTaulu, logml] = etsiParas(osuus, osuusTaulu, omaFreqs, logml); -# end -# refProportions(ind,:) = refProportions(ind,:).*(iter-1) + osuusTaulu; -# refProportions(ind,:) = refProportions(ind,:)./iter; -# end -# end -# for ind = 1:nrefIndsInPop -# omanOsuus = refProportions(ind,pop); -# if round(omanOsuus*100)==0 -# omanOsuus = 0.01; -# end -# if abs(omanOsuus)<1e-5 -# omanOsuus = 0.01; -# end -# refTaulu(pop, round(omanOsuus*100),level) = refTaulu(pop, round(omanOsuus*100),level)+1; -# end -# end -# end -# end - -# # Rounding of the results: -# proportionsIt = proportionsIt.*100; proportionsIt = round(proportionsIt); -# proportionsIt = proportionsIt./100; -# for ind = 1:ninds -# if ~any(to_investigate==ind) -# if PARTITION(ind)>0 -# proportionsIt(ind,PARTITION(ind))=1; -# end -# else -# # In case of a rounding error, the sum is made equal to unity by -# # fixing the largest value. -# if (PARTITION(ind)>0) & (sum(proportionsIt(ind,:)) ~= 1) -# [isoin,indeksi] = max(proportionsIt(ind,:)); -# erotus = sum(proportionsIt(ind,:))-1; -# proportionsIt(ind,indeksi) = isoin-erotus; -# end -# end -# end - -# # Calculate p-value for each individual: -# uskottavuus = zeros(ninds,1); -# for ind = 1:ninds -# pop = PARTITION(ind); -# if pop==0 # Individual is outlier -# uskottavuus(ind)=1; -# elseif isempty(find(to_investigate==ind)) -# # Individual had log-likelihood ratio<3 -# uskottavuus(ind)=1; -# else -# omanOsuus = proportionsIt(ind,pop); -# if abs(omanOsuus)<1e-5 -# omanOsuus = 0.01; -# end -# if round(omanOsuus*100)==0 -# omanOsuus = 0.01; -# end -# level = missing_level_partition(ind); -# refPienempia = sum(refTaulu(pop, 1:round(100*omanOsuus), level)); -# uskottavuus(ind) = refPienempia / nrefIndsInPop; -# end -# end - -# tulostaAdmixtureTiedot(proportionsIt, uskottavuus, alaRaja, iterationCount); - -# viewPartition(proportionsIt, popnames); - -# talle = questdlg(['Do you want to save the admixture results?'], ... -# 'Save results?','Yes','No','Yes'); -# if isequal(talle,'Yes') -# #waitALittle; -# [filename, pathname] = uiputfile('*.mat','Save results as'); - - -# if (filename == 0) & (pathname == 0) -# # Cancel was pressed -# return -# else # copy 'baps4_output.baps' into the text file with the same name. -# if exist('baps4_output.baps','file') -# copyfile('baps4_output.baps',[pathname filename '.txt']) -# delete('baps4_output.baps') -# end -# end - - -# if (~isstruct(tietue)) -# c.proportionsIt = proportionsIt; -# c.pvalue = uskottavuus; # Added by Jing -# c.mixtureType = 'admix'; # Jing -# c.admixnpops = npops; -# # save([pathname filename], 'c'); -# save([pathname filename], 'c', '-v7.3'); # added by Lu Cheng, 08.06.2012 -# else -# tietue.proportionsIt = proportionsIt; -# tietue.pvalue = uskottavuus; # Added by Jing -# tietue.mixtureType = 'admix'; -# tietue.admixnpops = npops; -# # save([pathname filename], 'tietue'); -# save([pathname filename], 'tietue', '-v7.3'); # added by Lu Cheng, 08.06.2012 -# end -# end + if (!isstruct(tietue)) { + c$proportionsIt <- proportionsIt + c$pvalue <- uskottavuus # Added by Jing + c$mixtureType <- 'admix' # Jing + c$admixnpops <- npops; + save(c, file=filename) + } else { + tietue$proportionsIt <- proportionsIt + tietue$pvalue <- uskottavuus; # Added by Jing + tietue$mixtureType <- 'admix' + tietue$admixnpops <- npops + save(tietue, file=filename) + } + } } \ No newline at end of file diff --git a/R/find.R b/R/find.R new file mode 100644 index 0000000..be417dd --- /dev/null +++ b/R/find.R @@ -0,0 +1,9 @@ +#' @title Find indices and values of nonzero elements +#' @description Emulates behavior of `find` +find <- function(x) { + if (is.logical(x)) { + return(which(x)) + } else { + return(which(x > 0)) + } +} \ No newline at end of file diff --git a/R/inputdlg.R b/R/inputdlg.R new file mode 100644 index 0000000..c803ddc --- /dev/null +++ b/R/inputdlg.R @@ -0,0 +1,17 @@ +#' @title Gather user input +#' @description Replicates the functionality of the homonymous function in Matlab (sans dialog box) +#' @param prompt Text field with user instructions +#' @param dim number of dimensions in the answwers +#' @param definput default value of the input +#' @export +inputdlg <- function(prompt, definput=NULL, dims=1) { + if (!is.null(definput)) { + prompt <- append(prompt, paste0(" (default: ", definput, ")")) + } + input_chr <- readline(paste0(prompt, ": ")) + if (input_chr == "") input_chr <- definput + input_chr_or_num <- tryCatch( + as.numeric(input_chr), warning = function(w) input_chr + ) + return(input_chr_or_num) +} \ No newline at end of file diff --git a/R/sortrows.R b/R/sortrows.R new file mode 100644 index 0000000..28c4d25 --- /dev/null +++ b/R/sortrows.R @@ -0,0 +1,15 @@ +#' @title Sort rows of matrix or table +#' @description Emulates the behavior of the `sortrows` function on Matlab +#' @param A matrix +#' @param column ordering column +sortrows <- function(A, column = 1) { + if (length(column) == 1) { + new_row_order <- order(A[, column]) + } else if (length(column) == 2) { + new_row_order <- order(A[, column[1]], A[, column[2]]) + } else { + stop("Not yet implemented for 2+ tie-breakers") + } + A_reordered <- A[new_row_order, ] + return(A_reordered) +} \ No newline at end of file diff --git a/R/viewPartition.R b/R/viewPartition.R new file mode 100644 index 0000000..80a3b73 --- /dev/null +++ b/R/viewPartition.R @@ -0,0 +1,121 @@ +viewPartition <- function(osuudet, popnames, COUNTS = matrix(0, 0, 0)) { + + npops <- size(COUNTS, 3) + nind <- size(osuudet,1) + +# TODO: translate if necessary. Remove if this function won't be used +# disp(['Number of populations: ' num2str(npops)]); +# if npops>30 +# disp(' '); +# disp('Figure can be drawn only if the number of populations'); +# disp('is less or equal to 30.'); +# disp(' '); +# return; +# end + + +# varit = givecolors(npops); +# korkeinviiva = 1.05; +# pieninarvo = -korkeinviiva; + + +# h0 = figure; +# set(h0, 'NumberTitle', 'off'); %image_figure; %Muutettu +# tiedot.popnames = popnames; +# tiedot.info = osuudet; +# set(h0,'UserData',tiedot); + +# set(gca, 'Xlim', [-.5 ,nind+.5], 'YLim', [pieninarvo ,korkeinviiva], ... +# 'XTick', [], 'XTickLabel', [], 'YTick', [], 'YTickLabel', []); + +# for i=1:nind + +# if any(osuudet(i,:)>0) +# cumOsuudet = cumsum(osuudet(i,:)); + +# % Pylv��n piirt�minen +# for j=1:npops +# if j==1 +# if cumOsuudet(1)>0 +# h0 =patch([i-1, i, i, i-1], [0, 0, cumOsuudet(1), cumOsuudet(1)], varit(j,:)); +# set(h0,'EdgeColor','none'); % Midevaa varten kommentoitava! +# end +# else +# if (cumOsuudet(j)>cumOsuudet(j-1)) +# h0 = patch([i-1, i, i, i-1], [cumOsuudet(j-1), cumOsuudet(j-1), ... +# cumOsuudet(j), cumOsuudet(j)], varit(j,:)); +# set(h0,'EdgeColor','none'); % Midevaa varten kommentoitava! +# end +# end +# end +# end +# end + + + +# if ~isempty(popnames) +# npops = size(popnames,1); +# for i=1:npops +# firstInd = popnames{i,2}; +# line([firstInd-1, firstInd-1], [0,1], 'Color', 'k'); %Populaatioiden rajat + +# if i 36) stop('Maximum number of colors 36') + colors <- matrix( + data = c( + 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, + 0.4, 0, 0, 0, 0.4, 0, 0, 0, 0.4, 0.4, 0.4, 0, 0.4, 0, + 0.4, 0, 0.4, 0.4, 0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2, 0.2, + 0.2, 0, 0.2, 0, 0.2, 0, 0.2, 0.2, 0.8, 0, 0, 0, 0.8, 0, + 0, 0, 0.8, 0.8, 0.8, 0, 0.8, 0, 0.8, 0, 0.8, 0.8, + 0.6, 0, 0, 0, 0.6, 0, 0, 0, 0.6, 0.6, 0.6, 0, 0.6, 0, + 0.6, 0, 0.6, 0.6, 0.6, 0.2, 0.4, 0.2, 0.4, 0.8, 0.8, + 0.4, 0.2, 0, 0.6, 0.2, 0.2, 0.8, 0.6, 0.5, 0.2, 0.1, + 0.6, 0.3, 0.1 + ), + ncol = 3, + byrow = TRUE + ) + colors = colors[1:n, ] + # red; green; blue; yellow + # RGB format: [red green blue] + return(colors) +} \ No newline at end of file diff --git a/R/zeros_ones.R b/R/zeros_ones.R index 2bc6f00..c4487bf 100644 --- a/R/zeros_ones.R +++ b/R/zeros_ones.R @@ -8,8 +8,15 @@ #' @note Actually works for any `x`, but there's no need to bother imposing #' validation controls here. zeros_or_ones <- function(n, x) { + # Expanding n to length 2 if necessary if (length(n) == 1) n <- c(n, n) - return(matrix(x, n[1], n[2])) + + # Returning a matrix or an array + if (length(n) == 2) { + return(matrix(x, n[1], n[2])) + } else { + return(array(x, dim=n)) + } } #' @title Matrix of zeros @@ -17,8 +24,14 @@ zeros_or_ones <- function(n, x) { #' the `zeros()` function on Matlab #' @param n1 number of rows #' @param n2 number of columns -zeros <- function(n1, n2 = n1) { - return(zeros_or_ones(c(n1, n2), 0)) +#' @param ... extra dimensions +zeros <- function(n1, n2 = n1, ...) { + if (length(n1) == 1) { + n <- c(n1, n2, ...) + } else { + n <- n1 + } + return(zeros_or_ones(n, 0)) } #' @title Matrix of ones @@ -26,6 +39,12 @@ zeros <- function(n1, n2 = n1) { #' the `ones()` function on Matlab #' @param n1 number of rows #' @param n2 number of columns -ones <- function(n1, n2 = n1) { - return(zeros_or_ones(c(n1, n2), 1)) +#' @param ... extra dimensions +ones <- function(n1, n2 = n1, ...) { + if (length(n1) == 1) { + n <- c(n1, n2, ...) + } else { + n <- n1 + } + return(zeros_or_ones(n, 1)) } \ No newline at end of file diff --git a/man/find.Rd b/man/find.Rd new file mode 100644 index 0000000..be4de89 --- /dev/null +++ b/man/find.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/find.R +\name{find} +\alias{find} +\title{Find indices and values of nonzero elements} +\usage{ +find(x) +} +\description{ +Emulates behavior of `find` +} diff --git a/man/inputdlg.Rd b/man/inputdlg.Rd new file mode 100644 index 0000000..35096d6 --- /dev/null +++ b/man/inputdlg.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/inputdlg.R +\name{inputdlg} +\alias{inputdlg} +\title{Gather user input} +\usage{ +inputdlg(prompt, definput = NULL, dims = 1) +} +\arguments{ +\item{prompt}{Text field with user instructions} + +\item{definput}{default value of the input} + +\item{dim}{number of dimensions in the answwers} +} +\description{ +Replicates the functionality of the homonymous function in Matlab (sans dialog box) +} diff --git a/man/ones.Rd b/man/ones.Rd index e739460..b110cf8 100644 --- a/man/ones.Rd +++ b/man/ones.Rd @@ -4,12 +4,14 @@ \alias{ones} \title{Matrix of ones} \usage{ -ones(n1, n2 = n1) +ones(n1, n2 = n1, ...) } \arguments{ \item{n1}{number of rows} \item{n2}{number of columns} + +\item{...}{extra dimensions} } \description{ wrapper of `zeros_or_ones()` that replicates the behavior of diff --git a/man/sortrows.Rd b/man/sortrows.Rd new file mode 100644 index 0000000..9680899 --- /dev/null +++ b/man/sortrows.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sortrows.R +\name{sortrows} +\alias{sortrows} +\title{Sort rows of matrix or table} +\usage{ +sortrows(A, column = 1) +} +\arguments{ +\item{A}{matrix} + +\item{column}{ordering column} +} +\description{ +Emulates the behavior of the `sortrows` function on Matlab +} diff --git a/man/zeros.Rd b/man/zeros.Rd index 039b95b..a5fa748 100644 --- a/man/zeros.Rd +++ b/man/zeros.Rd @@ -4,12 +4,14 @@ \alias{zeros} \title{Matrix of zeros} \usage{ -zeros(n1, n2 = n1) +zeros(n1, n2 = n1, ...) } \arguments{ \item{n1}{number of rows} \item{n2}{number of columns} + +\item{...}{extra dimensions} } \description{ wrapper of `zeros_or_ones()` that replicates the behavior of diff --git a/tests/testthat/test-convertedBaseFunctions.R b/tests/testthat/test-convertedBaseFunctions.R index 357ce70..1800e15 100644 --- a/tests/testthat/test-convertedBaseFunctions.R +++ b/tests/testthat/test-convertedBaseFunctions.R @@ -39,9 +39,11 @@ test_that("zeros and ones work as expected", { expect_equal(zeros(2), matrix(0, 2, 2)) expect_equal(zeros(2, 1), matrix(0, 2, 1)) expect_equal(zeros(1, 10), matrix(0, 1, 10)) + expect_equal(zeros(3, 2, 4), array(0, c(3, 2, 4))) expect_equal(ones(8), matrix(1, 8, 8)) expect_equal(ones(5, 2), matrix(1, 5, 2)) expect_equal(ones(2, 100), matrix(1, 2, 100)) + expect_equal(ones(3, 2, 4, 2), array(1, c(3, 2, 4, 2))) }) test_that("times works as expected", { @@ -143,4 +145,19 @@ test_that("isempty works as expected", { expect_false(isempty(cat1)) expect_true(isempty(cat2)) expect_false(isempty(str1)) +}) + +test_that("find works as expected", { + X <- matrix(c(1, 0, 2, 0, 1, 1, 0, 0, 4), 3, byrow=TRUE) + Y <- seq(1, 19, 2) + expect_equal(find(X), c(1, 5, 7, 8, 9)) + expect_equal(find(!X), c(2, 3, 4, 6)) + expect_equal(find(Y == 13), 7) +}) + +test_that("sortrows works as expected", { + mx <- matrix(c(3, 2, 2, 1, 1, 10, 0, pi), 4) + expect_equal(sortrows(mx), matrix(c(1, 2, 2, 3, pi, 10, 0, 1), 4)) + expect_equal(sortrows(mx, 2), matrix(c(2, 3, 1, 2, 0, 1, pi, 10), 4)) + expect_equal(sortrows(mx, 1:2), mx[order(mx[, 1], mx[, 2]), ]) }) \ No newline at end of file