Stats141 HW 7

Remember that homework where I thought myself an awesome programmer?

Check out the solutions:

source("http://www.bioconductor.org/biocLite.R")
biocLite("GEOquery")
Loading GEOquery (needed each time you re-start R)
library(GEOquery)
my.gds <- getGEO('GDS3257')
my.eset <- GDS2eSet(my.gds)
abstract(my.eset)
expressions = exprs(my.eset) # extracts the matrix of gene expressions
phenotypes = pData(my.eset) # extracts a data.frame containing the
annotation = fData(my.eset) # extracts the annotation information for the platform

pvalue=function(x){
    fit=lm(x~individual+tissue+individual:tissue+disease.state+gender, data=phenotypes)
    return(anova(fit)$Pr[5])
}
p.value=apply(expressions, 1, pvalue)
p.order=order(p.value)
FDR=round(p.value[p.order][100]*length(p.value)/100, 3)
rej=T
i=1
while(rej==T){
  if(p.value[p.order][i]>FDR*i/length(p.value)){
    print(i)
    select=i-1
    rej=F
  }else i=i+1
}
anno=annotation$Gene.title[p.order][1:select]
for (i in 1:select) print(anno[i])

## Why is this interesting? What distinguishes genes that are chosen by having a small p-value in this test?
# Genes choosen for small p-value are affected by not only individual, issus, gender
# and disease.state, but also the interaction between tissue and individual
# smoking status. In other words, the effect of the tissue type on the gene expression
# level also depends on whether the person is a smoking.

## Choose an FDR level so that the Benjamini-Hochberg procedure yields approximately 100 discoveries.
# With FDR around 0.23~0.24, approximately 100 genes will be selected.

So short, so...elegant, with so many functions truncated using one powerful new function instead.

I'm so jealous, and feel so inadequate.  But at least my analysis was correct?