-
Notifications
You must be signed in to change notification settings - Fork 29
/
tradis_essentiality.R
executable file
·78 lines (57 loc) · 3.55 KB
/
tradis_essentiality.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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#!/usr/bin/env Rscript
# PODNAME: tradis_essentiality.R
# ABSTRACT: tradis_essentiality.R
library("MASS")
options(warn=-1)
options(width=80)
args <- commandArgs(trailingOnly = TRUE)
input = args[1]
if( is.na(input) ){
cat(paste("Usage: tradis_essentiality.R data.tab\n\n"))
cat(strwrap("Produces calls of gene essentiality using an adaptation of the method described in Langridge et al. Genome Research 2009 and Barquist et al. NAR 2013. A loess curve is fit to the distribution of insertion indices, and used to identify the minima between the 'essential' and 'non-essential' distributions. These distributions are then used to fit gamma distributions, which are then used to calculate log-odds ratios, which are used to determine an insertion-index threshold for gene essentiality. Note that this analysis requires a saturated mutant library, and is not suitable for the analysis of data sets with low insertion density. The script produces a number of diagnostic plots which can be used to verify that this condition has been met.\n"), fill=80)
q(status=1)
}
STM_baseline <- read.table(input, sep="\t",header=TRUE,stringsAsFactors=F, quote="\"")
ii <- STM_baseline$ins_index
#identify second maxima
h <- hist(ii, breaks=200,plot=FALSE)
maxindex <- which.max(h$density[10:length(h$density)])
maxval <- h$mids[maxindex+3]
# print pdf of loess curve and later on, histogram
pdf(paste(input, "QC_and_changepoint_plots", "pdf", sep = "."))
#find inter-mode minima with loess
nG <- length(STM_baseline$read_count)
r <- floor(maxval *2000)
I = ii < r / 2000
h1 = hist(ii[I],breaks=(0:r/2000))
lo <- loess(h1$density ~ c(1:length(h1$density))) #loess smothing over density
plot(h1$density, main="Density")
lines(predict(lo),col='red',lwd=2)
m = h1$mids[which.min(predict(lo))]
I1 = ((ii < m)&(ii >= 0))
h = hist(ii, breaks="FD",plot=FALSE)
I2 = ((ii >= m)&(ii < h$mids[max(which(h$counts>5))]))
f1 = (sum(I1) + sum(ii == 0))/nG
f2 = (sum(I2))/nG
d1 = fitdistr(ii[I1], "exponential")
d2 = fitdistr(ii[I2], "gamma") #fit curves
# print pdf of histogram
#pdf("Loess_and_changepoint_estimation.pdf")
#plots
hist(ii,breaks="FD", xlim=c(0,max(ii)), freq=FALSE,xlab="Insertion index", main="Gamma fits")
lines(0:200/500, f1*dgamma(0:200/500, 1, d1$estimate[1])) # was [2]
lines(0:200/500, f2*dgamma(0:200/500, d2$estimate[1], d2$estimate[2]))
# print changepoint
#calculate log-odds ratios to choose thresholds
lower <- max(which(log((pgamma(1:1000/10000, d2$e[1],d2$e[2])*(1-pgamma(1:1000/10000, 1,d1$e[1], lower.tail=FALSE)))/(pgamma(1:1000/10000, 1,d1$e[1], lower.tail=FALSE)*(1-pgamma(1:1000/10000, d2$e[1],d2$e[2]))) , base=2) < -2))
upper <- min(which(log((pgamma(1:1000/10000, d2$e[1],d2$e[2])*(1-pgamma(1:1000/10000, 1,d1$e[1], lower.tail=FALSE)))/(pgamma(1:1000/10000, 1,d1$e[1], lower.tail=FALSE)*(1-pgamma(1:1000/10000, d2$e[1],d2$e[2]))) , base=2) > 2))
essen <- lower/10000
ambig <- upper/10000
lines(c(lower/10000, lower/10000), c(0,20), col="red")
lines(c(upper/10000, upper/10000), c(0,20), col="red")
mtext(paste(essen, ":", "Essential changepoint"), side=3, adj=1, padj=2)
mtext(paste(ambig, ":", "Ambiguous changepoint"), side=3, adj=1, padj=3.75)
dev.off()
write.csv(STM_baseline, file=paste(input, "all", "csv", sep="."), row.names = FALSE, col.names= TRUE, quote=FALSE)
write.csv(STM_baseline[STM_baseline$ins_index < essen,], file=paste(input, "essen", "csv", sep="."), row.names = FALSE, col.names= TRUE, quote=FALSE)
write.csv(STM_baseline[STM_baseline$ins_index >= essen & STM_baseline$ins_index < ambig,], file=paste(input, "ambig", "csv", sep="."), row.names = FALSE, col.names= TRUE, quote=FALSE)