2014-04-06

Script to compute bp of a band in a gel electrophoresis picture

mm - milimeter from where sample was put into gele
bp - base pair


You need:
  1. Text file with marker sizes and distances (marker.csv)
  2. Text file with distance of your sample (sample.csv)

# R Script file get_bp_from_mm.R:

getbpfrommm <- function(marker, unknown, unknown.filename){
  # log bp of marker
  markerlog <- marker
  markerlog$bp <- log(markerlog$bp)
 
  # make linear model
  marker.lm <- lm(markerlog$bp ~ markerlog$mm)
 
  b=coef(marker.lm)[1]
  m=coef(marker.lm)[2]
 
  unknown.bp <- c()
 
  for(unknown.mm.i in unknown){
    unknown.bp.i = round(exp(m*unknown.mm.i+b))
    unknown.bp <- c(unknown.bp, unknown.bp.i)}
 
  unknown.bp.mm <- cbind(unknown, unknown.bp)
  colnames(unknown.bp.mm) <- c("mm", "bp")
 
  print(unknown.bp.mm)
 
  write.table(unknown.bp.mm, file=paste(unknown.filename, "with_bp.csv", sep="_"), quote=FALSE, sep="\t", row.names=FALSE, col.names=TRUE)
}

#read command line arguments
args <- commandArgs(trailingOnly = TRUE)
marker <- read.delim(args[1])
unknown <- read.delim(args[2])

getbpfrommm(marker, unknown, args[2])





# Invoke from command line
$ cat marker.csv
mm    bp
132.60    75
116.80    200
106.20    300
95.60    400
89.20    500
72.20    700
62.80    1000
48.30    1500
38.40    2000
19.1    5000
 

$ cat sample.csv
mm
103.4
98.8 


$ Rscript get_bp_from_mm.R marker.csv sample.csv
     mm  bp
1 103.4 272
2  98.8 317


$ cat sample.csv_with_bp.csv
mm    bp
103.4    272
98.8    317

No comments:

Post a Comment