massiR_Vignette.R
In massiR: massiR: MicroArray Sample Sex Identifier

### R code from vignette source 'massiR_Vignette.Rnw'

###################################################
### code chunk number 1: load the example data
###################################################
library(massiR)
data(massi.test.dataset)


###################################################
### code chunk number 2: load the test probe list
###################################################
data(massi.test.probes)


###################################################
### code chunk number 3: extract Y from test dataset calculate CV
###################################################
massi.y.out <- massi_y(massi.test.dataset, massi.test.probes)


###################################################
### code chunk number 4: plot the results from massi.y
###################################################
massi_y_plot(massi.y.out)


###################################################
### code chunk number 5: fig1too
###################################################
massi.y.out <- massi_y(massi.test.dataset, massi.test.probes)


###################################################
### code chunk number 6: fig1
###################################################

  barplot(height=massi.y.out[[2]], names.arg=massi.y.out[[1]], xpd=T,
          cex.names=0.5, las=2, ylab="Probe CV (%)")
  # Get the quantile values from the massi.y output
  quantiles <- massi.y.out[[3]]
  
  # add lines for the 0%, 25%, 50%, and 75% quartiles
  abline(h=quantiles[1:4], col=c("black", "red", "blue", "green"), lwd=2)
  legend("topleft",cex=0.7, title="Threshold (Quantile)",
         col=c("black", "red", "blue", "green"),
         fill= c("black", "red", "blue", "green"),
         legend=c("1 (0%)", "2 (25%)",
                  "3 (50%)", "4 (75%)"))


###################################################
### code chunk number 7: run massi.select
###################################################
massi.select.out <- 
  massi_select(massi.test.dataset, massi.test.probes, threshold=4)


###################################################
### code chunk number 8: head massi select out
###################################################
head(massi.select.out)[,1:5]


###################################################
### code chunk number 9: massiR_Vignette.Rnw:109-110
###################################################
results <- massi_cluster(massi.select.out)


###################################################
### code chunk number 10: massiR_Vignette.Rnw:114-115
###################################################
sample.results <- data.frame(results[[2]])


###################################################
### code chunk number 11: massiR_Vignette.Rnw:118-119
###################################################
head(sample.results)


###################################################
### code chunk number 12: massiR_Vignette.Rnw:129-130 (eval = FALSE)
###################################################
## massi_cluster_plot(massi.select.out, results)


###################################################
### code chunk number 13: fig2
###################################################
ord <- order(rowSums(abs(massi.select.out)),decreasing=T)
  heatmap.2(x=as.matrix(massi.select.out[ord,]), keysize=2, cexRow=0.7,
            key=T, trace="none", dendrogram="row", col=redgreen(75), scale="row")


###################################################
### code chunk number 14: fig3
###################################################
massi.cluster.results <- data.frame(results[[2]])
  massi.cluster.results.sort <- massi.cluster.results[order(massi.cluster.results$sex),] # sort data by sex
  probe.means <- massi.cluster.results.sort$mean_y_probes_value # samples probe mean values
  probe.sd <- massi.cluster.results.sort$y_probes_sd # sample probe sd values
  sample.names <- massi.cluster.results.sort$ID # set x-axis names
  plot.top <- ceiling(max(probe.means+probe.sd*1.1)) # set y-axis upper limit
  plot.bottom <- floor(min(probe.means-probe.sd*1.1)) # set y-axis lower limit
  sample.sex <- massi.cluster.results.sort$sex # set the factor for bar color
  # create the plot
  barCenters <- barplot(probe.means, xpd=F, names.arg=results$ID, cex.names=0.7,
                        ylab="Chr.Y mean probe value +/- SD",
                        xlab="",
                        col=c("red", "green")[as.factor(sample.sex)],
                        las=2, ylim=c(plot.bottom,plot.top))
  segments(barCenters, probe.means-probe.sd, # add the sd bars
           barCenters, probe.means+probe.sd, lwd=0.8)
  legend("topleft", fill=c("red", "green"), title="predicted sex", ## add legend to plot
         legend=c("female", "male"), cex=0.5, )


###################################################
### code chunk number 15: fig4
###################################################
## generate PC plot of clusters
  k.medoids.results <- results[[1]]
  clusplot(t(massi.select.out), k.medoids.results$clustering, color=TRUE, shade=FALSE, main="",cex.txt=0.5,
           labels=2, lines=0)


###################################################
### code chunk number 16: massi.dip
###################################################
dip.result <- massi_dip(massi.select.out)


###################################################
### code chunk number 17: fig5
###################################################
dip.result <- massi_dip(massi.select.out)
plot(dip.result[[3]])


###################################################
### code chunk number 18: fig6
###################################################
dip.result <- massi_dip(massi.select.out)
hist(x=dip.result[[2]], breaks=20)


###################################################
### code chunk number 19: bias_dip
###################################################
male.ids <-
  subset(sample.results$ID,
                   subset=sample.results$sex=="male")

female.ids <-
  subset(sample.results$ID,
         subset=sample.results$sex=="female")


###################################################
### code chunk number 20: massiR_Vignette.Rnw:246-249
###################################################
bias.subset.ids <- c(female.ids[1:18], male.ids[1:2])

bias.subset <- massi.select.out[bias.subset.ids]


###################################################
### code chunk number 21: massiR_Vignette.Rnw:252-253
###################################################
bias.dip <- massi_dip(bias.subset)


###################################################
### code chunk number 22: massiR_Vignette.Rnw:264-265
###################################################
data(massi.eset, massi.test.probes)


###################################################
### code chunk number 23: massiR_Vignette.Rnw:269-274
###################################################
eset.select.out <- 
  massi_select(massi.eset, massi.test.probes)

eset.results <- 
  massi_cluster(eset.select.out)


###################################################
### code chunk number 24: massiR_Vignette.Rnw:278-297
###################################################
# Get the sex for each sample from the massi_cluster results
eset.sample.results <- 
  data.frame(eset.results[[2]])

sexData <- 
  data.frame(eset.sample.results[c("ID", "sex")])

# Extract the order of samples in the ExpressionSet and match with results
eset.names <- 
  colnames(exprs(massi.eset))

# match the sample order in massiR results to the same as the ExpressionSet object
sexData <- sexData[match(eset.names, sexData$ID),]

# create an annotatedDataFrame to add to ExpressionSet
pData <- new("AnnotatedDataFrame", data = sexData)

# add the annotatedDataFrame to the Expressionset as phenoData
phenoData(massi.eset) <- pData


###################################################
### code chunk number 25: massiR_Vignette.Rnw:302-304
###################################################
# check the phenodata is now within the ExpressionSet
phenoData(massi.eset)


###################################################
### code chunk number 26: massiR_Vignette.Rnw:307-310
###################################################
# check that all phenodata id's match expressionSet column names.
# This must return "TRUE"
all(massi.eset$ID == colnames(exprs(massi.eset)))


###################################################
### code chunk number 27: load the included probe lists
###################################################
 data(y.probes)


###################################################
### code chunk number 28: y.probes names
###################################################
  names(y.probes)


###################################################
### code chunk number 29: massiR_Vignette.Rnw:327-328
###################################################
  illumina.v2.probes <- data.frame(y.probes["illumina_humanwg_6_v2"])


###################################################
### code chunk number 30: massiR_Vignette.Rnw:340-351 (eval = FALSE)
###################################################
## library(biomaRt)
## mart <- useMart('ensembl', dataset="hsapiens_gene_ensembl")
## filters <- listFilters(mart)
## attributes <- listAttributes(mart)
## 
## gene.attributes <- 
##  getBM(mart=mart, values=TRUE,
##        filters=c("with_illumina_humanwg_6_v2"),
##        attributes= c("illumina_humanwg_6_v2", "entrezgene",
##                       "chromosome_name", "start_position",
##                       "end_position", "strand"))


###################################################
### code chunk number 31: massiR_Vignette.Rnw:354-356 (eval = FALSE)
###################################################
## unique.probe <- 
##   subset(gene.attributes, subset=!duplicated(gene.attributes[,1]))


###################################################
### code chunk number 32: massiR_Vignette.Rnw:359-361 (eval = FALSE)
###################################################
## y.unique <- 
##   subset(unique.probe, subset=unique.probe$chromosome_name == "Y")


###################################################
### code chunk number 33: massiR_Vignette.Rnw:365-367 (eval = FALSE)
###################################################
## illumina.v2.probes <- 
##   data.frame(row.names=y.unique$illumina_humanwg_6_v2)