gistfile1.r

# Create a vector of sex-specific positive control genes
sex.symbols <- c("XIST","CYorf15A","DDX3Y","KDM5D","RPS4Y1","USP9Y","UTY");
sex.genes <- featureNames(dataset)[match(sex.symbols, fData(dataset)$Symbol)];
names(sex.genes) <- sex.symbols;
sex.genes <- na.omit(sex.genes);

# Draw stripcharts of expression of each sex-specific gene
pdf.filename <- file.path(project.path, sprintf("%s_sex_stripcharts.pdf", project.prefix));
if (!file.exists(pdf.filename)) {
	pdf(pdf.filename, width=11, height=8.5);
	layout(matrix(1:length(sex.genes), nrow=1));
	for (gene in names(sex.genes)) {
		p <- t.test(exprs(dataset)[sex.genes[gene],] ~ dataset$sex, var.equal=TRUE)$p.value;
		fold.change <-
			mean(2^(exprs(dataset)[sex.genes[gene], dataset$sex=="M"])) /
			mean(2^(exprs(dataset)[sex.genes[gene], dataset$sex=="F"]));
		if (fold.change < 1) fold.change <- -1 / fold.change;
		# stripchart is a pain in the butt;
		# to properly assign separate colors to each point, each class must be plotted separately
		# Set up a blank stripchart
		stripchart(
			exprs(dataset)[sex.genes[gene],] ~ dataset$sex,
			method="jitter", vertical=TRUE,
			ylab="log2 (expression)", col=NA,
			xlim=c(0.25,nlevels(dataset$sex)+0.75),
			ylim=range(exprs(dataset)[sex.genes,]),
			main=sprintf("%s\nFC = %0.2f, p = %0.2g", gene, fold.change, p)
		);
		# Then draw points for each sample
		for (i in 1:nlevels(dataset$sex)) {
			j <- which(dataset$sex == levels(dataset$sex)[i]);
			stripchart(
				exprs(dataset)[sex.genes[gene], j],
				method="jitter", vertical=TRUE, add=TRUE, at=i,
				pch=21, col="black", cex=1.5, bg=dataset$color[j]
			);
		}
	}
	dev.off();
}

# Density plot of GSTT1 expression
pdf.filename <- file.path(project.path, sprintf("%s_GSTT1_density.pdf", project.prefix));
if (!file.exists(pdf.filename)) {
	pdf(pdf.filename, width=11, height=8.5);
	plot(density(exprs(dataset)["2952_at",]), xlab="log2 (expression)", main="Distribution of GSTT1 expression");
	dev.off();
}

heatmap.legend <- "Colors scaled by row/gene (blue/white/red=below/at/above average)";