R tips

If you are looking for R tips for beginners and advanced users there is nothing better than this:
http://pj.freefaculty.org/R/Rtips.html

or in .pdf version:
http://pj.freefaculty.org/R/Rtips.pdf

Topics include:
- exchange data between R and other programs (Excel, etc)
- create variable names on the fly
- return multiple values from a function
- scatterplot: smooth a line connecting points
- test for Normality
- multiple analysis of variance
- character encoding
- R environment in side scripts

and many more. Enjoy.

10 papers that every bioinformatician should read

This list of top 10 papers in bioinformatics does not only consist of purely bioinformatic papers, but shows how gentle and simple bioinformatics can and should be:

1. Dudley JT, Butte AJ (2009) A Quick Guide for Developing Effective Bioinformatics Programming Skills. PLoS Computational Biology 5(12)
Link
"Consequently it is no surprise that many successful bioinformatics apps are written by biologists who lack formal computer science training, as they undoubtedly put scientific utility ahead of architectural elegance and completeness."

2. Kim TK, Hemberg M, Gray JM, Costa AM, Bear DM, Wu J, Harmin DA, Laptewicz M, Barbara-Haley K, Kuersten S, Markenscoff-Papadimitriou E, Kuhl D, Bito H, Worley PF, Kreiman G, Greenberg ME (2010) Widespread transcription at neuronal activity-regulated enhancers. Nature 465(7295)
Link
Simple and brilliant use of ChIP-seq and RNA-seq techniques.

3. Goecks J, Nekrutenko A, Taylor J, and The Galaxy Team (2010) Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome Biology 11(8)
Link
Bioinformatics for the masses.

4. Ameur A, Zaghlool A, Halvardson J, Wetterbom A, Gyllensten U, Cavelier L, Feuk L (2011) Total RNA sequencing reveals nascent transcription and widespread co-transcriptional splicing in the human brain. Nature Structural & Molecular Biology 18
Link
Great example of getting additional features from RNA-seq data.

5. Smedley D, Haider S, Ballester B, Holland R, London D, Thorisson G, Kasprzyk A (2009) BioMart – biological queries made easy. BMC Genomics 10
Link
Data mart for biologists.

6. Krzywinski M, Schein J, Birol I, Jones S, Marra M (2008) Circos - an information aesthetic for comparative genomics. Genome Informatics conference
Link
Beauty of data visualization.

7. Majewski J, Ott J (2004) Distribution and characterization of regulatory elements in the human genome. Genome Research 12
Link
One of the first successful attempts to analyze genomic features at genome scale.

8. Kent WJ, Zweig AS,  Barber G, Hinrichs AS, Karolchik D (2009) BigWig and BigBed: enabling browsing of large distributed datasets. Bioinformatics 26(17)
Link
Standard formats for storing NGS-data. Fast, simple.


The list will continue...

Plotting ROC curves in ggplot2

Default ROC curves in R are disgusting. ggplot2 comes to the rescue.

First, let's write some data generating function that will be useful for ROC:

getExampleDataForROC <- function(len = 100, distort = 1) {
  half <- round(len/2)
  ### Group ###
  group <- vector()
  group[1:half] <- "Disease"
  group[(half + 1):len] <- "Control"
  group <- as.factor(group)
  ### Score ###
  score <- vector()
  score[1:half] <- 1
  score[(half + 1):len] <- -1
  for (i in 1:len)
    score[i] <- score[i] + rnorm(1, 0, distort)
  order <- order(score)
  data.frame(pos = score[order], group = group[order])
}

Let's write ROC plotting functions:

getPointsForROC <- function(len = 100, distort = 1) {
  basal <- getExampleDataForROC(len, distort)
  tp <- vector(); tn <-vector(); fp <-vector(); fn <- vector()
  tpr <- vector(); fpr <- vector()
  acc <- vector(); spc <- vector()
  len <- dim(basal)[1]
  for(i in 1:len-1) {
    tp[i] <- sum(basal[(i+1):len,2] == "Disease")
    tn[i] <- sum(basal[1:i,2] == "Control")
    fp[i] <- sum(basal[(i+1):len,2] == "Control")
    fn[i] <- sum(basal[1:i,2] == "Disease")
    tpr[i] <- tp[i] / (tp[i] + fn[i])
    fpr[i] <- fp[i] / (fp[i] + tn[i])
    acc[i] <- (tp[i] + tn[i]) / ((tp[i] + fn[i]) +  (fp[i] + tn[i]))
    spc[i] <- 1 - fpr[i]
  }
  points <- (cbind(fpr,tpr))[(len-1):1,]
  points <- rbind(points, c(1,1))
  return(points)
} 


plotROC <- function(len = 100) {
  points.part1 <- getPointsForROC(len, 1)
  desc.part1 <- rep("Distort = 1",len)
  points.part2 <- getPointsForROC(len, 1.5)
  desc.part2 <- rep("Distort = 1.5",len)
  points.part3 <- getPointsForROC(len, 2)
  desc.part3 <- rep("Distort = 2",len)
  points.part4 <- getPointsForROC(len, 2.5)
  desc.part4 <- rep("Distort = 2.5",len)

  points <- rbind(points.part1, points.part2, points.part3, points.part4)
  desc <- c(desc.part1, desc.part2, desc.part3, desc.part4)
  data <- data.frame(TPR = points[,2], FPR = points[,1], GeneSet = desc)
  colors = brewer.pal(5, "YlOrRd")[2:5]
  qplot(FPR, TPR, data = data, geom = "blank", main = "ROC curve", xlab = "False Positive Rate (1-Specificity)", ylab = "True Positive Rate (Sensitivity)" ) + geom_line(aes(x = FPR, y = TPR, data = data, colour = GeneSet), size = 2, alpha = 0.7) + scale_colour_manual(values=colors)
}

Finally, let's generate plot:

plotROC(1000)

Plotting PCA results in ggplot2

Default PCA plots in R are disgusting. ggplot2 comes to the rescue.

First, let's write some data generating functions that will be useful for PCA:

getVector <- function(n = 20, dims = 2) {
  out <- vector()
  breaks <- round((n/dims)*(1:dims))
  start <- 1  
  for(i in 1:dims) {
    num <- rnorm(1, 10, 1)
    for(j in start:(breaks[i])) {
      out[j] <- num + rnorm(1, 0, 0.2)
    }
    start <- breaks[i] + 1
  } 
  out
}

getData <- function(rows = 10, cols = 20, groups = 2, dims = 2, distort = 0.5) {
  out <- matrix(data = 10, nrow = rows, ncol = cols)
  ### make groups ###
  breaks <- round((rows/groups)*(1:groups))
  start <- 1  
  for(i in 1:groups) {
      vector <- getVector(n = cols, dim = 2)
    for(j in start:(breaks[i])) {
      out[j,] <- vector
    }
    start <- breaks[i] + 1
  }
  ### make error ###
  for(i in 1:rows) {
    out[i,] <- out[i,] + rnorm(n = cols, mean = 0, sd = distort)
  }
  ### return ###
  rownames(out) <- paste("Row ", 1:rows, sep = "")
  colnames(out) <- paste("Col ", 1:cols, sep = "")
  out
}

getGroup <- function(rows = 10, groups = 2) {
  out <- vector()
  breaks <- round((rows/groups)*(1:groups))
  start <- 1  
  for(i in 1:groups) {
    for(j in start:(breaks[i])) {
      out[j] <- i
    }
    start <- breaks[i] + 1
  }
  out <- as.factor(out)
  out
}

Let's write PCA plotting function:

makePcaPlot <- function(x = getData(), group = getGroup(), title = "") {
  require(ggplot2)
  require(RColorBrewer)
  data <- x
  data <- t(apply(data, 1, scale))
  rownames(data) <- rownames(x)
  colnames(data) <- colnames(x)
  mydata <- t(data)  
  groupFactor <- group
  mydata.pca <- prcomp(mydata, retx=TRUE, center=TRUE, scale.=TRUE)  

  percent <- round((((mydata.pca$sdev)^2 / sum(mydata.pca$sdev^2))*100)[1:2])
  loadings <- mydata.pca$rotation
  rownames(loadings) <- colnames(mydata)
  scores <- mydata.pca$x
  cex = 3
  group <- as.character(groupFactor)
  what <- groupFactor == "1"; group[what] <- "One"
  what <- groupFactor == "2"; group[what] <- "Two"
  group <- factor(group)
  group <- factor(group, levels=c("Two", "One"), ordered=TRUE)
  colors = sample(1:11, size = length(levels(group)))


  PCA1 <- scores[,1]
  PCA2 <- scores[,2]
  base_size = 10

  qplot(PCA2, PCA1, geom="blank", main = title, xlab = paste("PCA2 (", percent[2], "%)", sep = ""), ylab = paste("PCA1 (", percent[1], "%)", sep = "")) + geom_point(aes(colour = group), size = 6, alpha = 3/4) + scale_colour_manual(values=brewer.pal(11, "RdYlGn")[11:1][colors]) + opts(
      axis.text.x = theme_text(size = base_size * 1.3 , lineheight = 0.9, colour = "grey50", hjust = 1, angle = 90),
      axis.text.y =       theme_text(size = base_size * 1.3, lineheight = 0.9, colour = "grey50", hjust = 1)
      )
}

Finally, let's generate plot:

makePcaPlot(getData(30,4,2,distort = 0.7), getGroup(30,2))