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
|
maPlot <- function(x,y, logAbundance=NULL, logFC=NULL, normalize=FALSE, plot.it=TRUE, smearWidth=1, col=NULL, allCol="black", lowCol="orange", deCol="red", de.tags=NULL, smooth.scatter=FALSE, lowess=FALSE, ...)
# Low-level function for creating an MA-plot for DGE data.
# Created by Mark Robinson. Last modified by Davis McCarthy, 19 November 2010.
# Edits by Gordon Smyth 20 March 2011.
{
if( !is.null(logAbundance) && !is.null(logFC) ) {
A <- logAbundance
M <- logFC
w <- v <- rep(FALSE, length(A))
w <- A < -25+log2(1e6) # logCPM instead of logConc
if( any(w) ) {
shift <- max(abs(M[w])) - max(abs(M[!w]))
A[w] <- min(A[!w]) - runif(sum(w),min=0,max=smearWidth)
M[w] <- sign(M[w]) * (abs(M[w]) - shift)
}
} else {
if(normalize) {
x <- x/sum(x)
y <- y/sum(y)
}
A <- (log2(x)+log2(y))/2
M <- log2(y) - log2(x)
w <- x==min(x) | y==min(y)
if( any(w) ) {
A[w] <- min(A[!w]) - runif(sum(w),min=0,max=smearWidth)
M[w] <- log2(y[w]+min(y[!w])) - log2(x[w]+min(x[!w]))
}
}
qs <- quantile(M, c(0.05,0.95))
range <- qs[2]-qs[1]
v <- (M < (median(M) - 5*range)) | (M > (median(M) + 5*range))
if( any(v) ) {
M[v] <- sign(M[v]) * (max(abs(M[!v])) + 0.5*range)
}
if(plot.it) {
if( is.null(col) ) {
col <- rep(allCol, length(A))
if( any(w) | any(v) )
col[w | v] <- lowCol
}
if(smooth.scatter) {
smoothScatter(A, M, col=col, ...)
grid()
if( any(w) | any(v) )
points(A[w | v], M[w | v], col=lowCol, ...)
}
else
plot(A,M,col=col,...)
points(A[de.tags],M[de.tags],col=deCol,...)
if(lowess) {
keep <- A > min(A[!(w |v)]) + 1
low <- lowess(A[keep],M[keep], f=1/4)
lines(low,col="red",lwd=4)
}
}
invisible(list(A=A,M=M,w=w,v=v))
}
|
