R/fdrFunctions.R

In SwathXtend: SWATH extended library generation and statistical data analysis

fdr.crit = function(dswat.fdr) {
	dswat.fdr = dswat.fdr[dswat.fdr$Decoy=='FALSE',]; 
	dswat.fdr$nfdr.pass = apply(dswat.fdr[,-c(1:7)] <= 0.01, 1, function(y) sum(na.omit(y)) );

	dswat.fdr
}


getFdrBins = function(mat.fdr, Bins = c(0, .01, .1, .2, .3, .4, .5, .8, 1)) {
	v = 100*sum(as.vector(mat.fdr) >= Bins[1] & 
				as.vector(mat.fdr) <= Bins[2]) /length(as.vector(mat.fdr))
				
	for(bb in 2:8)
		v = c(v, 100*sum(as.vector(mat.fdr) > Bins[bb] & 
			as.vector(mat.fdr) <= Bins[bb+1])/length(as.vector(mat.fdr)) )
	
	v		
}

			
coverage = function(a, b)
{
	100*round(length(intersect(a,b))/length(unique(a)),3)

}


bland.altman = function(v1,v2)
{
	round(abs(v1-v2)/mean(c(v1,v2)),3)
}

cv = function(v)
{
	round(sd(na.omit(v) )/mean(na.omit(v)),4)
}

# sample-by-sample quantification accuracy/reproducbility of two swath data
quantification.accuracy = function(dswat1, dswat2, Sample=NULL, method=c('cor','cv','bland.altman'), 
	cor.method=c('pearson', 'spearman', 'kendall'),
	log=FALSE)
{

	method = match.arg(method)
	cor.method = match.arg(cor.method)
	
	stopifnot(ncol(dswat1) == ncol(dswat2)) 
	
	stopifnot(all(names(dswat1) == names(dswat2)))	
	
	if(is.null(Sample)) Sample = colnames(dswat1)[-c(1:5)]		
	
	nsample = nlevels(factor(Sample))
	
	vcor = rep(NA,nsample)
	rcor = rep(NA,nsample) #randomised test
	
	dswat1.means = swath.means(dswat1, Sample)
	dswat2.means = swath.means(dswat2, Sample)
	
	m = merge(dswat1.means, dswat2.means, by=1, all=F)

	m = m[,-1]
	
	m.1 = m[,1:nsample]
	
	m.2 = m[,-c(1:nsample)]
	
	stopifnot(nrow(m.1) == nrow(m.2) & ncol(m.1) == ncol(m.2))	
	
	if(log) {
		m.1 = log(m.1)
		m.2 = log(m.2)
	}
	
	if(method == 'cor') {
	for(ii in 1:nsample) vcor[ii] = round(cor(m.1[,ii],m.2[,ii], method=cor.method),2)
				set.seed(1)
	for(ii in 1:nsample) rcor[ii] = round(cor(m.1[,ii], m.2[,sample((1:nsample)[-ii],1)], , method=cor.method),2)
	} else if(method == 'cv') {
		for(ii in 1:nsample) vcor[ii] = median(sapply(1:nrow(m.1), function(x) cv(c(m.1[x,ii],m.2[x,ii]))) )
			set.seed(1)
		for(ii in 1:nsample) rcor[ii] = median(sapply(1:nrow(m.1), function(x) cv(c(m.1[,ii], 
			m.2[,sample((1:nsample)[-ii],1)]))) )
	} else { # bland-altman
	for(ii in 1:nsample) vcor[ii] = median(sapply(1:nrow(m.1), function(x) bland.altman(m.1[x,ii],m.2[x,ii])) )
				set.seed(1)
	for(ii in 1:nsample) rcor[ii] = median(sapply(1:nrow(m.1), function(x) bland.altman(m.1[,ii], 
			m.2[,sample((1:nsample)[-ii],1)])))
	}

	list(vcor=vcor, rcor=rcor)
	
}



swath.means = function(dswath, Sample=NULL) 
{

	if(is.null(Sample)) Sample = colnames(dswath)[-c(1:5)]

	mat = dswath[,-c(1:5)]
	
	agg = aggregate(mat, by=list(dswath$Protein), function(x) sum(na.omit(x)))

	log.means = aggregate(t(agg[,-1]), by=list(Sample), function(x) mean(log(na.omit(x))) )

	means = exp(log.means[,-1])
	
	
	dat.means = data.frame(Protein=agg[,1], t(means) )

	colnames(dat.means)[-1] = log.means[,1]
		
	dat.means
}