R/applyttestPep.R

In SwathXtend: SWATH extended library generation and statistical data analysis

applyttestPep <-
function(peptides, Group, doLogs=TRUE, numerator=levels(as.factor(Group))[1]) {

# columns: protein, peptide, values

Protein = peptides[,1]
peptides = peptides[,-c(1:2)]

Group = as.factor(as.vector(Group))

pval = rep(NA, nrow(peptides))
FC = rep(NA, nrow(peptides))

if (doLogs == TRUE) peptides = log(peptides+1);
if (!(nlevels(Group) == 2)) stop("T test needs two levels only in group.");

mn.1 =  apply(peptides[,Group == numerator], 1, FUN=function(v){mean(na.omit(v))})
mn.2 =  apply(peptides[,Group == setdiff(levels(Group), numerator)], 1, FUN=function(v){mean(na.omit(v))})

# calculate FC by peptide
if (doLogs) {
	FC = exp(mn.1)/exp(mn.2)
} else {
	FC = mn.1/mn.2
}

# calculate pval by protein from peptide ratios
ag.peptide.pval = aggregate(FC, by=list(Protein=Protein), 
		FUN=function(v){ res=NA;tp=try(t.test(log(v))); 
		if (!inherits(tp, "try-error")) res=tp$p.value; res})

		
# calculate FC by protein		
ag.peptide.FC = aggregate(FC, by=list(Protein=Protein), 
		FUN=function(v){ exp(mean(log(na.omit(v))))})		

if ( sum(ag.peptide.FC[,1] != ag.peptide.pval[,1]) > 0) stop("Error in calculating peptide level tests.")

res = data.frame(ag.peptide.pval, ag.peptide.FC[,-1])
colnames(res) = c("Protein", "pval", "FC")

res[match(unique(Protein), res[,1]),]

}