commit bugfix

Author: snaketron

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: IgGeneUsage
 Type: Package
 Title: Differential gene usage in immune repertoires
-Version: 1.17.20
+Version: 1.17.21
 Authors@R: 
     person(given = "Simo",  
            family = "Kitanovski",

R/dgu.R

@@ -18,8 +18,7 @@ DGU <- function(ud,
                   hdi_lvl = hdi_lvl,
                   paired = paired)
 
-  udr <- ud
-  ud <- get_usage(u = udr)
+  ud <- get_usage(u = ud)
 
   # setup control list
   control_list <- list(adapt_delta = adapt_delta,

R/utils_dgu.R

@@ -1,7 +1,7 @@
 get_contrast_map <- function(ud) {
   # condition map
-  condition_map <- data.frame(condition_name = ud$condition_names,
-                              condition_id = ud$condition_id)
+  condition_map <- data.frame(condition_name = ud$condition_name_of_sample,
+                              condition_id = ud$condition_id_of_sample)
   condition_map <- condition_map[duplicated(condition_map)==FALSE,]
   rownames(condition_map) <- condition_map$condition_id
 

R/utils_gu.R

@@ -37,7 +37,7 @@ get_condition_prop_dgu <- function(glm, hdi_lvl, ud) {
   gu_summary$condition <- NA
   for(i in 1:ud$N_condition) {
     gu_summary$condition[which(gu_summary$condition_id == i)]<-
-      ud$condition_names[which(ud$condition_id == i)[1]]
+      ud$condition_name_of_sample[which(ud$condition_id_of_sample == i)[1]]
   }
 
   # remove unused vars
@@ -71,7 +71,7 @@ get_condition_prop_gu <- function(glm, hdi_lvl, ud) {
 
   gu_summary$gene_id <- as.numeric(par)
   gu_summary$gene_name <- ud$gene_names[gu_summary$gene_id]
-  gu_summary$condition <- ud$condition_names[1]
+  gu_summary$condition <- ud$condition_name_of_sample[1]
 
   # remove unused vars
   gu_summary$gene_id <- NULL

R/utils_usage.R

@@ -2,6 +2,36 @@
 # Description
 # Parse input data for GU analysis
 get_usage <- function(u) {
+  
+  get_condition_of_sample <- function(sample_ids, 
+                                      m) {
+    # condition at sample
+    condition_names <- character(length = length(sample_ids))
+    for(i in 1:length(sample_ids)) {
+      condition_names[i] <- m$condition[m$sample_id == sample_ids[i]]
+    }
+    condition_ids <- as.numeric(as.factor(condition_names))
+    
+    return(list(condition_ids = condition_ids,
+                condition_names = condition_names))
+  }
+  
+  get_condition_of_individual <- function(individual_ids, 
+                                          individual_names, 
+                                          m) {
+    # condition at individual
+    condition_names <- character(length = max(individual_ids))
+    for(i in 1:max(individual_ids)) {
+      q <- individual_names[individual_ids==i][1]
+      condition_names[i] <- m$condition[m$individual_id == q][1]
+    }
+    condition_ids <- as.numeric(as.factor(condition_names))
+    
+    return(list(condition_ids = condition_ids,
+                condition_names = condition_names))
+  }
+  
+  
   u$individual_org_name <- u$individual_id
   if("replicate" %in% colnames(u)) {
 
@@ -89,32 +119,22 @@ get_usage <- function(u) {
   tr <- table(replicate_ids)
   has_balanced_replicates <- ifelse(test = all(tr==tr[1]), yes=TRUE, no=FALSE)
 
-  if(has_replicates) {
-    # condition at individual
-    condition_names <- character(length = max(individual_ids))
-    for(i in 1:max(individual_ids)) {
-      q <- individual_names[individual_ids==i][1]
-      condition_names[i] <- m$condition[m$individual_id == q][1]
-    }
-    condition_ids <- as.numeric(as.factor(condition_names))
-  } else {
-    # condition at sample
-    condition_names <- character(length = length(sample_ids))
-    for(i in 1:length(sample_ids)) {
-      condition_names[i] <- m$condition[m$sample_id == sample_ids[i]]
-    }
-    condition_ids <- as.numeric(as.factor(condition_names))
-  }
+  cos <- get_condition_of_sample(sample_ids = sample_ids, m = m) 
+  coi <- get_condition_of_individual(individual_ids = individual_ids,
+                                     individual_names = individual_names,
+                                     m = m)
 
   return(list(Y = Y, 
               N = as.numeric(N), 
               N_sample = ncol(Y), 
               N_gene = nrow(Y),
               gene_names = gene_names,
               sample_names = sample_ids, 
-              condition_id = condition_ids,
-              condition_names = condition_names,
-              N_condition = max(condition_ids),
+              condition_id_of_sample = cos$condition_ids,
+              condition_name_of_sample = cos$condition_names,
+              condition_id_of_individual = coi$condition_ids,
+              condition_name_of_individual = coi$condition_names,
+              N_condition = max(cos$condition_ids),
               individual_id = individual_ids,
               individual_names = individual_names,
               individual_org_names = individual_org_names,
@@ -125,10 +145,12 @@ get_usage <- function(u) {
               N_replicate = max(replicate_ids),
               proc_ud = u,
               has_replicates = has_replicates,
-              has_conditions = max(condition_ids)>1,
+              has_conditions = max(cos$condition_ids)>1,
               has_balanced_replicates = has_balanced_replicates))
 }
 
+
+
 # Description:
 # get the appropriate model
 get_model <- function(has_replicates, 
@@ -173,10 +195,8 @@ get_model <- function(has_replicates,
       if(paired) {
         model <- stanmodels$dgu_paired_rep
         pars <- c("phi", "kappa", "alpha", 
-                  "sigma_condition", "sigma_individual", "sigma_alpha",
-                  "sigma_alpha_rep", "sigma_beta_rep",
-                  "alpha_sample", "beta_sample", 
-                  "alpha_individual", "beta_individual", 
+                  "sigma_condition", "sigma_individual", "sigma_beta_rep",
+                  "mu_individual", "beta_sample", "beta_individual", 
                   "beta_condition", 
                   "Yhat_rep", "Yhat_rep_prop", "Yhat_condition_prop", 
                   "log_lik", "dgu", "dgu_prob", "theta")

inst/stan/dgu.stan

@@ -28,19 +28,16 @@ functions {
 }
 
 data {
-  int<lower=0> N_sample;                   // number of repertoires 
-  int<lower=0> N_gene;                     // gene
-  int<lower=0> N_individual;               // number of individuals
-  int<lower=0> N_condition;                // number of conditions
-  array [N_sample] int N;                  // number of tries (repertoire size)
-  array [N_gene, N_sample] int Y;          // number of heads for each coin
-  array [N_sample] int condition_id;       // id of conditions
-  array [N_sample] int individual_id;      // id of replicate
+  int<lower=0> N_sample;                        // # samples
+  int<lower=0> N_gene;                          // # genes
+  int<lower=0> N_condition;                     // # conditions
+  array [N_sample] int N;                       // rep size
+  array [N_gene, N_sample] int Y;               // gene usage
+  array [N_sample] int condition_id_of_sample;  // sample condition id
 }
 
 transformed data {
-  // convert int N -> real N fo convenient division
-  // in generated quantities block
+  // convert int to real N for convenient division
   array [N_sample] real Nr;
   Nr = N;
 }
@@ -68,7 +65,7 @@ transformed parameters {
   }
 
   for(i in 1:N_sample) {
-    beta_individual[i]  = beta_condition[condition_id[i]] + sigma_individual[condition_id[i]] * z_beta_individual[i];
+    beta_individual[i]  = beta_condition[condition_id_of_sample[i]] + sigma_individual[condition_id_of_sample[i]] * z_beta_individual[i];
     theta[i] = inv_logit(alpha + beta_individual[i]);
   }
 }

inst/stan/dgu_paired.stan

@@ -28,14 +28,14 @@ functions {
 }
 
 data {
-  int<lower=0> N_sample;                   // number of repertoires 
-  int<lower=0> N_gene;                     // gene
-  int<lower=0> N_individual;               // number of individuals
-  int<lower=0> N_condition;                // number of conditions
-  array [N_sample] int N;                  // number of tries
-  array [N_gene, N_sample] int Y;          // number of heads for each coin
-  array [N_sample] int condition_id;       // id of conditions
-  array [N_sample] int individual_id;      // id of replicate
+  int<lower=0> N_sample;                        // # samples 
+  int<lower=0> N_gene;                          // # genes
+  int<lower=0> N_individual;                    // # individuals
+  int<lower=0> N_condition;                     // # conditions
+  array [N_sample] int N;                       // rep size
+  array [N_gene, N_sample] int Y;               // gene usage
+  array [N_sample] int condition_id_of_sample;  // condition of sample
+  array [N_sample] int individual_id;           // individual id
 }
 
 transformed data {
@@ -70,7 +70,7 @@ transformed parameters {
   }
 
   for(i in 1:N_sample) {
-    theta[i] = inv_logit(alpha + beta_individual[individual_id[i]] + beta_condition[condition_id[i]]);
+    theta[i] = inv_logit(alpha + beta_individual[individual_id[i]] + beta_condition[condition_id_of_sample[i]]);
   }
 }
 

inst/stan/dgu_paired_rep.stan

@@ -28,14 +28,14 @@ functions {
 }
 
 data {
-  int<lower=0> N_sample;                   // number of repertoires 
-  int<lower=0> N_gene;                     // gene
-  int<lower=0> N_individual;               // number of individuals
-  int<lower=0> N_condition;                // number of conditions
-  array [N_sample] int N;                  // number of tries
-  array [N_gene, N_sample] int Y;          // number of heads for each coin
-  array [N_sample] int condition_id;       // id of conditions
-  array [N_sample] int individual_id;      // id of replicate
+  int<lower=0> N_sample;                        // # samples 
+  int<lower=0> N_gene;                          // # genes
+  int<lower=0> N_individual;                    // # individuals
+  int<lower=0> N_condition;                     // # conditions
+  array [N_sample] int N;                       // rep size
+  array [N_gene, N_sample] int Y;               // gene usage
+  array [N_sample] int condition_id_of_sample;  // condition of sample
+  array [N_sample] int individual_id;           // replicate id
 }
 
 transformed data {
@@ -73,7 +73,7 @@ transformed parameters {
 
   for(i in 1:N_sample) {
     beta_sample[i]  = beta_individual[individual_id[i]] + sigma_beta_rep * z_beta_sample[i];
-    theta[i] = inv_logit(alpha + beta_sample[i] + beta_condition[condition_id[i]]);
+    theta[i] = inv_logit(alpha + beta_sample[i] + beta_condition[condition_id_of_sample[i]]);
   }
 }
 

inst/stan/dgu_rep.stan

@@ -28,14 +28,14 @@ functions {
 }
 
 data {
-  int<lower=0> N_sample;                   // number of repertoires
-  int<lower=0> N_gene;                     // gene
-  int<lower=0> N_individual;               // number of individuals
-  int<lower=0> N_condition;                // number of conditions
-  array [N_sample] int N;                   // number of tries (repertoire size)
-  array [N_gene, N_sample] int Y;           // number of heads for each coin
-  array [N_individual] int condition_id;    // id of conditions
-  array [N_sample] int individual_id;       // id of replicate
+  int<lower=0> N_sample;                                // # samples
+  int<lower=0> N_gene;                                  // # genes
+  int<lower=0> N_individual;                            // # individuals
+  int<lower=0> N_condition;                             // # conditions
+  array [N_sample] int N;                               // rep size
+  array [N_gene, N_sample] int Y;                       // gene usage size
+  array [N_individual] int condition_id_of_individual;  // id of conditions
+  array [N_sample] int individual_id;                   // id of replicate
 }
 
 transformed data {
@@ -70,7 +70,7 @@ transformed parameters {
   }
 
   for(i in 1:N_individual) {
-    beta_individual[i] = beta_condition[condition_id[i]] + sigma_individual[condition_id[i]] * z_beta_individual[i];
+    beta_individual[i] = beta_condition[condition_id_of_individual[i]] + sigma_individual[condition_id_of_individual[i]] * z_beta_individual[i];
   }
 
   for(i in 1:N_sample) {

inst/stan/gu.stan

@@ -30,8 +30,6 @@ functions {
 data {
   int<lower=0> N_sample;                   // number of repertoires
   int<lower=0> N_gene;                     // gene
-  int<lower=0> N_individual;               // number of individuals
-  int<lower=0> N_condition;                // number of conditions
   array [N_sample] int N;                  // number of tries (repertoire size)
   array [N_gene, N_sample] int Y;          // number of heads for each coin
 }

inst/stan/gu_rep.stan

@@ -31,7 +31,6 @@ data {
   int<lower=0> N_sample;                    // number of repertoires
   int<lower=0> N_gene;                      // gene
   int<lower=0> N_individual;                // number of individuals
-  int<lower=0> N_condition;                 // number of conditions
   array [N_sample] int N;                   // number of tries (repertoire size)
   array [N_gene, N_sample] int Y;           // number of heads for each coin
   array [N_sample] int individual_id;       // id of replicate