# Program for testing if an SNP represents a synthetic association
# Written in R, http://cran.r-project.org/
#
# 21 Sept 2010
# Contact: Fumihiko Takeuchi, http://www.fumihiko.takeuchi.name
# Reference: Detection of common single nucleotide polymorphisms synthesizing quantitative trait association of rarer causal variants (submitted)

# Input:
# q: a vector representing quantitative trait value for each individual
# x: a vector representing genotype for each individual
# Testing synthetic association:
# Load the functions below, and try synthetic.test(q, x)

skewness.test = function(q, x) {
	xlist = setdiff(x, NA);
	chilist = rep(NA, length(xlist));
	for (i in 1:length(xlist)) {
		xi = xlist[i];
		qxi = q[x==xi];
		qxi = qxi[!is.na(qxi)];
		chilist[i] = length(qxi) / 6 *
			(sum((qxi - mean(qxi))^3)/length(qxi))^2 / 
			(sum((qxi - mean(qxi))^2)/length(qxi))^3
	}
	pchisq(sum(chilist), df=length(chilist), lower.tail=F)
}
bartlett.skewness.test = function(q, x) {
	pbi = bartlett.test(q, x)$p.value;
	psi = skewness.test(q, x);
	pchisq(-2 * log(pbi * psi), df=4, lower.tail=F) 
}
synthetic.test = function(q, x) {
	#rank-based inverse normal transformation
	c = 3/8;
	qinr = (rank(q, ties.method="average") - c)/(sum(!is.na(q)) - 2 * c + 1);
	qinr[is.na(q)] = NA;
	qinr = qnorm(qinr);
	print(paste("Heteroscedasticity p-value: ",
		bartlett.test(qinr, x)$p.value));
	print(paste("Skewness p-value: ",
		skewness.test(qinr, x)));
	print(paste("Combined p-value: ",
		bartlett.skewness.test(qinr, x)));
}