Jean,

Thank you for your suggestion it has really helped a lot. I ended up changing the expression system in my model, since I realised that if the "genome" has extensive duplication, i.e. same values, which my organism of interest is know to evolve, then because of the ux <- unique(c(x, y))part of the function, duplicate genes within the genome may mutate in the same way, to test this I changed the gene values to a much narrower range between 1 and 10, rather than 10000. What I found was all the 7's would mutate negatively by 1, but all the 5's would mutate positively by 1 etc. I also realised with duplication that's extensive, if it's value comes up in the expression data, then which of those duplicated genes is the one that's expressed is not clear, which complicates things.

I've changed it so as the expressed data is the same length as the genes data, but consists of a 1 and 0 for on and off. The function decided how many effectors to express between 1 and the total number of genes, samples them, makes a data series the same length as the effector genes, with all values set to 0 or off. Then gets the indexes of the genes sampled from the genome to be switched on, and then switches the 0's at the corresponding indexes in the expression data to 1.


Init_Express <- function(x, min.exp=1, max.exp=length(x)) {
  Number_Expressed <- round(runif(1, min=min.exp, max=max.exp))
  Expressed <- sample(x,Number_Expressed)
  E <- rep(0, length(x))
  pos<-match(Expressed, x)
  E[pos]<-1
  return(E)
}

Expression.State <- lapply(seq(Effectors), function(i) Init_Express(Effectors[[i]]))

Based on how you managed to mutate genes using operations without loops, I thought of how to do this with expression states:

Change.Expression <- function(x, prob.express=0.5, prob.repress=0.5) {
  old<-x
  low.old <- old[which(old==0)]
  low.new <- low.old + sample(c(0, 1), size=length(low.old), replace=TRUE, prob=c(1-prob.express, 0+prob.express))
  high.old <- old[which(old==1)]
  high.new <- high.old - sample(c(0, 1), size=length(high.old), replace=TRUE, prob=c(1-prob.repress, 0+prob.repress))
  all.new <- old
  all.new[pmatch(low.old, old, dup=F)] <- low.new
  all.new[pmatch(high.old, old, dup=F)] <- high.new
  return(all.new)
}

Expression.State <- lapply(seq(Expression.State), function(i) Change.Expression(Expression.State[[i]]))

I've heard a lot of talk about the benefits or drawbacks of loops vs. the apply family of functions, be it speed or tidyness of code. I've been offered use of a high performance cluster in my dept. (not sure how R does on those) but I was still concerned about nested loops because of the number of individuals and then genes per individual, I mostly went with loops because of the indexing issue - controlling going through individuals, and then all the genes in turn. Changing the nested part of the loop in the manner you've demonstrated, eliminates that and will be a big help. It's helped me see how I can use lapply with functions that do things in a more vectorised manner and a more R like manner, I was afraid of the eventuality of having to write some processed in C++ to do them quicker.

Best Wishes,

Ben W.

UEA (ENV) and The Sainsbury Laboratory.






________________________________
From: Jean V Adams [jvadams@usgs.gov]
Sent: 27 November 2012 22:01
To: Benjamin Ward (ENV)
Cc: r-help@r-project.org
Subject: Re: [R] Stuck trying to modify a function

Ben,

You can use the sample() function to randomly add -1, 0, or 1 to each observation, and control for the probability of mutation at the same time.  Then you can use the match() function to make sure that any mutations in X are carried through to Y in the same way.  I wrote the function to do each list element separately.  So a gene in X[[1] and Y[[1]] will be mutated in the same way, but the same gene in X[[2]] and Y[[2]] may be mutated in a different way.  Not sure if that is what you want.

Mutate <- function(x, y, prob.mutate=0.9) {
        ux <- unique(c(x, y))
        new.ux <- ux + sample(c(-1, 0, 1), size=length(ux), replace=TRUE, prob=c(prob.mutate/2, 1-prob.mutate, prob.mutate/2))
        new.x <- new.ux[match(x, ux)]
        new.y <- new.ux[match(y, ux)]
        list(xm=new.x, ym=new.y)
        }

Effectors <- lapply(seq(X), function(i) Mutate(X[[i]], Y[[i]]))

Jean



"Benjamin Ward (ENV)" <B.Ward@uea.ac.uk> wrote on 11/27/2012 10:45:23 AM:
>
> Hi,
>
> I have the following data:
>
> Path_Number <- 5
> ID.Path <- c(1:Path_Number) # Make vector of ID's.
> No_of_X <- sample(50:550, length(ID.Path), replace=TRUE) #
> X <- split(sample(1:10000, sum(No_of_X), replace=TRUE), rep(ID.Path,No_of_X))
> Y <- lapply(X,function(x) sample(x, round(runif(1, min=10, max=50))))
>
> X and Y are both lists, and I've made the following function to work
> on that data as part of a simulation I'm building:
>
> Mutate<-function(x){
>   l<-0
>   for(i in x){
>     l2<-0
>     l<-l+1
>     for(i in x[[l]]){
>       l2<-l2+1
>       if(runif(1) < 0.9) ifelse(runif(1) <0.5, x[[l]][l2] <- x[[l]]
> [l2]+1, x[[l]][l2] <- x[[l]][l2]-1)
>     }
>   }
>   return(x)
> }
>
> I call this with Effectors<-Mutate(X)
> The function is designed to alter the values of each element in X by
> either + or - 1 (50:50 chance wether + or -). However Y, elements of
> which are a subset of the corresponding elements of X, need to be
> consistent i.e. if a value in X is changed, and that value is part
> of the Y subset, then the value in Y also needs to be changed.
> however, since Y is a smaller subset it will not be indexed the
> same. My idea was to include in the function an if statement that
> checks if Y contains the value to be changed, removes it, and then
> after the value in X is changed, put the new value in Y. I attemptedthis with:
>
>
> Mutate<-function(x,y){
>   l<-0
>   for(i in x){
>     l2<-0
>     l<-l+1
>     for(i in x[[l]]){
>       l2<-l2+1
>       if(runif(1) < 0.9){
>         if(x[[l]][l2] %in% y[[l]] == TRUE){
>           y[[l]]<-[which(y[[l]]!=x[[l]][l2])]
>           if(runif(1) <0.5){
>             x[[l]][l2] <- x[[l]][l2]+1
>             y[[l]]<-append(x[[l]][l2])
>           }else{
>             x[[l]][l2] <- x[[l]][l2]-1
>             y[[l]]<-append(x[[l]][l2])
>           }
>         }
>         ifelse(runif(1) <0.5, x[[l]][l2] <- x[[l]][l2]+1, x[[l]][l2]
> <- x[[l]][l2]-1)
>       }
>     }
>   }
>   return(list(x,y))
> }
>
> Bit of an eyesore so I've put the altered stuff in bold. I've
> basically taken what the ifelse statement does in the first
> function, (which is still there and run if Y does not contain the X
> value being altered) and broken it down into an if and an else
> segment with multiple operations in curly braces to accommodate the
> extra actions needed to alter Y as well as X. This was all I could
> think of to keep changes between the two "in sync", however this
> does not work when I try to load the function into workspace:
>
> Error: unexpected '}' in "}"
>
> I hope someone can point out what it is I've done that isn't
> working, or a better way to do this.
>
> Best Wishes,
>
> Ben W.
>
> UEA (ENV) & The Sainsbury Laboratory.

	[[alternative HTML version deleted]]


