bp - base pair
You need:
- Text file with marker sizes and distances (marker.csv)
- 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