Student 1

Student 1 wrote this answer to Question 1:

library(seqinr)
sequences <- read.fasta("NC_000913.fna")

genes_gc_content <- function(genome) {
  positions <- seq(from=1, to=length(genome)-window_size, by= window_size)
  nA<-0
  nG<-0
  nC<-0
  nT<-0
  
  for(i in 1:length(genome)) {
    positions[i]<- genome_gc_content(genome[[i]])
    if(genome[i]=="a") {
      na <- na+1
    }else if(genome[i]=="t") {
      nt <- nt+1
    }else if(genome[i]=="g") {
      ng <- ng+1
    }else if(genome[i]=="c") {
      nc <- nc+1
    }
  }
  answer<- (ng+nc)/(ng+nc+nt+na)
  return(positions)
}

The first two lines are right, but they are not necessary. You can assume that the user of the function will do that work. Remember that there are two scenarios: first when you define the function, and second when someone else uses the function.

The line positions <- seq(from=1, to=length(genome)-window_size, by= window_size) is for Question 2 and is useless here.

Then we have the line positions[i]<- genome_gc_content(genome[[i]]). For one side, if we use genome[[i]] it means that genome is a list. But on the next lines we use genome as a vector. So, what is it?

Now, read again the text of the question. The function must be called window_gc_content, not genes_gc_content as here. Also, it takes three inputs: sequence, position, and size.

There is a detail in the if() commands. Here we test only for lower case letters, but as we saw on classes, sometimes they are upper-case letters.

Student 2

Student 2 wrote something like this:

gc_content<- function(coli_seq){
  A<- sum(coli_seq=="a"|| "A")
  T<- sum(coli_seq=="t"|| "T")
  G<- sum(coli_seq=="g"|| "G")
  C<- sum(coli_seq=="c"|| "C")
  ans<- (G+C)/(G+C+A+T)
  return(ans)
}

This does not work, because we are mixing apples and oranges. The part coli_seq=="a" is a comparison between a vector of letters and a letter; the result is a logic vector. Then the part "A" is a single letter, not a vector, not a logic value.

As a human I can understand that the student wants to say ‘if the sequence is “a” or “A”’, but the computer is dumb and does not understand. We need to say ‘if the sequence is “a” or the sequence is “A”’.

Use help(``||``) to se the manual Moreover, the operand || may not be the correct one, since it returns a single value. Better use |, as indicated in the help page.

window_gc_content <- function (genome, position, size){
  genome <- seq(from=250000, length=100)
  position <- 250000
  size <- 100
  nA <- 0
  nC <- 0
  nG <- 0
  nT <- 0
  for(i in 1:length(seq(from=250000, length=100))) {
    if(genome[i]=="A" | genome[i]=="a"){
      nA <- nA+1
    } else  if(genome[i]=="C" | genome[i]=="c"){
      nC <- nC+1
    } else  if(genome[i]=="G" | genome[i]=="g"){
      nG <- nG+1
    } else  if(genome[i]=="T" | genome[i]=="t"){
      nT <- nT+1
    }
  }
  answer <- (nG+nC)/(nA+nT+nG+nC)
 
  return(answer)
}

#'title: "Homework 5"
#'subtitle: "Computing in Molecular Biology and Genetics II"
#'author: "Put your name here"
#'number: STUDENT_NUMBER
#'date: "March 11, 2020"

window_gc_content <- function(sequence, position, size) {
  sequen <- rep("a", 1000)
  sequen <- genome[[1]][seq(from=position, length=size)]
  for(i in 1: length(sequen)) {
    sequen [i] <- gc_content(sequen[i])
    answer <- gc_content(sequne)
  }
  return(sequne)
}

Student 5

At the last minute I received this:

window_gc_content <- function(sequence, position, size)
  for(i in 1:length(genome)) {
    if(genome[i]=="A"| genome[i]=="a") {
      nA <- nA+1
    } else if (genome[i]=="C" | genome[i]=="c"){
      nC <- nC+1
    } else if (genome[i]=="G" | genome[i]=="g") {
      nG <- nG+1
    } else if (genome[i]=="T" | genome[i]== "t") {
      nT <- nT+1
    }
  }
  answer <- (nG+nC)/(nA+nT+nG+nC)

  answer <- rep("a", 1000)
  answer <- seq(from=sequences[[1]][250000], length=1000)
  return(answer)
}

The first line is nice, but it does not have the { symbol. The second one uses genome, which is not an input nor something we created. It appeared from nowhere… it is magic. In other words, it is wrong, since magic does not exist.

Again, the three inputs are never used in the code. They should be used.

Finally, answer is assigned three times. The first time it gets a number, which is what it should be. The second time we destroy the good answer and we replace it with a vector 1000 letters "a". Finally we destroy it again and we assign a vector of numbers, that do not depend on the GC content.

Moreover, the question says that sequences is a vector of character, not a list. So we cannot use sequences[[1]], which is only valid for lists.

Finally, the function seq() can only take numeric inputs. In other words, when we say seq(from=x, length=y), then x and y must be numbers. But sequences[[1]][250000] is not a number.

By the way, the numbers 1000 and 250000 should not be in your answer. They are examples. Your function should not depend on them. Always think on the general case, not the particular one.

Conclusion

• Be sure to understand the question. If you do not understand, ask in the forum. Explain what do you understand and what you do not understand.
• This is LEGO. Identify all the pieces and understand how they connect.
  • Write the name of each variable, its structure (vector, list, data frame, etc.) and its type (numeric, logic, character, etc.).
• You can only use the inputs and variables that you create. If you did not create it, do not change it.
• The output should change if the input change. Check that each input is used somewhere.
• Understand what you have and what you want to have. It is like a biochemistry process.
• Sometimes it is useful to work backwards. Start with what you want to have, and decompose it in simpler terms. For example, to get the GC-skew, you need to know nG and nC. The you just need to find nG and nC.
• Look at previous examples and recycle them, by using the old functions inside the new function.
  • If you cannot recycle, then you can adapt the old code for the new case.
• It is always wise to return to our very first class: How to Solve It.