-
Notifications
You must be signed in to change notification settings - Fork 51
/
make_error_models.R
65 lines (55 loc) · 2.01 KB
/
make_error_models.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
## create error models to be released wtih polyester
## from the GemSIM/GemErr estimation model
modfolder = '../polyester_paper/error_models'
platforms = c('ill100v4_mate1', 'ill100v4_mate2',
'ill100v4_single', 'ill100v5_mate1', 'ill100v5_mate2',
'ill100v5_single', 'r454ti_single')
for(platform in platforms){
i = which(platforms == platform)
model = read.table(paste(modfolder, platform, sep='/'), header=TRUE)
names(model)[1] = 'refbase'
model$refbase = substr(model$refbase, 4, 4)
eval(parse(text=paste0('model', i, ' <- model')))
save(list=paste0('model',i), file=paste0('data/', platform, '.rda'),
compress='xz')
}
## figures for paper:
## put one example in true manuscript,
## include the others in supplementary data.
library(RSkittleBrewer)
library(usefulstuff)
colrs = RSkittleBrewer('tropical')
getColor = function(nt){
switch(nt, A='black', C=colrs[1], G=colrs[2], T='deeppink', N=colrs[4])
}
nts = c('A','C','G','T','N')
plot_nt = function(model, nt){
d = subset(model, refbase==nt)
wrongnts = nts[-which(nts==nt)]
errCols = paste0('read', wrongnts)
mnum = as.matrix(model[,2:6])
ymax = max(mnum[mnum<0.8])
plot(1:100, as.matrix(d[errCols[1]])[-1], col=makeTransparent(getColor(wrongnts[1])),
pch=19, cex=0.5, xlab='Read Position', ylab='Error Probability',
type='o', ylim=c(0,ymax))
for(i in 2:4){
points(1:100, as.matrix(d[errCols[i]])[-1], pch=19, cex=0.5, type='o',
col=makeTransparent(getColor(wrongnts[i])))
}
title(nt)
legend('topleft', wrongnts, pch=19, col=sapply(wrongnts, getColor), cex=0.7)
}
plot_model = function(model, file){
pdf(file)
par(mfrow=c(2,2))
for(nt in c('A','C','G','T')){
plot_nt(model, nt)
}
dev.off()
}
plot_model(model1, 'illumina4_mate1.pdf')
plot_model(model2, 'illumina4_mate2.pdf')
plot_model(model3, 'illumina4_single.pdf')
plot_model(model4, 'illumina5_mate1.pdf')
plot_model(model5, 'illumina5_mate2.pdf')
plot_model(model6, 'illumina5_single.pdf')