PSEA_RNAmixtures.R
In PSEA: Population-Specific Expression Analysis.

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

###################################################
### code chunk number 1: style-Sweave
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: PSEA_RNAmixtures.Rnw:30-33
###################################################
library(PSEA)
library(GEOquery)
library(affy)


###################################################
### code chunk number 3: PSEA_RNAmixtures.Rnw:38-42 (eval = FALSE)
###################################################
## dataset<-getGEO(GEO="GSE19380",destdir=".")
## information<-pData(phenoData(dataset[["GSE19380_series_matrix.txt.gz"]]))
## sample_IDs<-as.character(information[,"geo_accession"])
## datafiles<-sapply(sample_IDs,function(x){rownames(getGEOSuppFiles(x))})


###################################################
### code chunk number 4: PSEA_RNAmixtures.Rnw:47-52 (eval = FALSE)
###################################################
## dataset<-getGEO(GEO="GSE19380",filename="GSE19380_series_matrix.txt.gz")
## dataset<-list("GSE19380_series_matrix.txt.gz"=dataset)
## information<-pData(phenoData(dataset[["GSE19380_series_matrix.txt.gz"]]))
## sample_IDs<-as.character(information[,"geo_accession"])
## datafiles<-file.path(sample_IDs,paste(sample_IDs,".CEL.gz",sep=""))


###################################################
### code chunk number 5: PSEA_RNAmixtures.Rnw:57-59 (eval = FALSE)
###################################################
## raw_data<-ReadAffy(filenames=datafiles,compress=TRUE)
## expression_GSE19380<-2^exprs(rma(raw_data))


###################################################
### code chunk number 6: PSEA_RNAmixtures.Rnw:64-65
###################################################
data(expression_GSE19380)


###################################################
### code chunk number 7: PSEA_RNAmixtures.Rnw:70-71
###################################################
expression<-expression_GSE19380[1:31042,]


###################################################
### code chunk number 8: PSEA_RNAmixtures.Rnw:76-79
###################################################
neuron_probesets<-list(c("1370058_at","1370059_at"),"1387073_at","1367845_at")
astro_probesets<-list("1372190_at","1386903_at",c("1375120_at","1375183_at","1385923_at"))
oligo_probesets<-list("1398257_at","1368861_a_at",c("1368263_a_at","1370434_a_at","1370500_a_at"))


###################################################
### code chunk number 9: PSEA_RNAmixtures.Rnw:84-88
###################################################
mixedsamples<-c(17:24)
neuron_reference<-marker(expression,neuron_probesets,sampleSubset=mixedsamples,targetMean=100)
astro_reference<-marker(expression,astro_probesets,sampleSubset=mixedsamples,targetMean=100)
oligo_reference<-marker(expression,oligo_probesets,sampleSubset=mixedsamples,targetMean=100)


###################################################
### code chunk number 10: figNeuronalReferenceSignal
###################################################
par(cex=0.7)
plot(neuron_reference,type="l")


###################################################
### code chunk number 11: PSEA_RNAmixtures.Rnw:100-102
###################################################
model1<-lm(expression["1367660_at",]~neuron_reference+astro_reference+
oligo_reference,subset=mixedsamples)


###################################################
### code chunk number 12: fig1367660_at_NAO
###################################################
par(mfrow=c(1,3),cex=0.7)
crplot(model1,"neuron_reference",newplot=FALSE)
crplot(model1,"astro_reference",newplot=FALSE,ylim=c(-250,250))
crplot(model1,"oligo_reference",newplot=FALSE,ylim=c(-250,250))


###################################################
### code chunk number 13: PSEA_RNAmixtures.Rnw:116-117
###################################################
summary(model1)


###################################################
### code chunk number 14: PSEA_RNAmixtures.Rnw:122-123
###################################################
model_matrix<-cbind(intercept=1,neuron_reference,astro_reference,oligo_reference)


###################################################
### code chunk number 15: PSEA_RNAmixtures.Rnw:128-129
###################################################
model_subset<-em_quantvg(c(2,3,4), tnv=3, ng=1)


###################################################
### code chunk number 16: PSEA_RNAmixtures.Rnw:134-135
###################################################
models<-lmfitst(t(expression), model_matrix, model_subset, subset=mixedsamples)


###################################################
### code chunk number 17: PSEA_RNAmixtures.Rnw:140-145
###################################################
regressor_names<-as.character(1:4)
coefficients<-coefmat(models[[2]], regressor_names)
pvalues<-pvalmat(models[[2]], regressor_names)
models_summary<-lapply(models[[2]], summary)
adjusted_R2<-slt(models_summary, 'adj.r.squared')


###################################################
### code chunk number 18: PSEA_RNAmixtures.Rnw:150-153
###################################################
negativecoefficient<-apply(coefficients[,-1]<0 & pvalues[,-1]<0.05,1,function(x){any(x,na.rm=TRUE)})
average_expression<-apply(expression[,mixedsamples], 1, mean)
filter<-!negativecoefficient & (coefficients[,1] / average_expression) < 0.5 & adjusted_R2 > 0.6


###################################################
### code chunk number 19: PSEA_RNAmixtures.Rnw:159-160
###################################################
sum(!negativecoefficient)


###################################################
### code chunk number 20: PSEA_RNAmixtures.Rnw:164-165
###################################################
sum(coefficients[,1] / average_expression < 0.5)


###################################################
### code chunk number 21: PSEA_RNAmixtures.Rnw:169-170
###################################################
sum(adjusted_R2 > 0.6)


###################################################
### code chunk number 22: PSEA_RNAmixtures.Rnw:174-175
###################################################
sum(filter)


###################################################
### code chunk number 23: PSEA_RNAmixtures.Rnw:180-182
###################################################
selectedpsname<-"1370431_at"
selectedps<-which(rownames(expression)==selectedpsname)


###################################################
### code chunk number 24: PSEA_RNAmixtures.Rnw:187-188
###################################################
coefficients[selectedps,2]


###################################################
### code chunk number 25: PSEA_RNAmixtures.Rnw:193-194
###################################################
pvalues[selectedps,2]


###################################################
### code chunk number 26: PSEA_RNAmixtures.Rnw:199-200 (eval = FALSE)
###################################################
## crplot(models[[2]][[selectedps]],"2",ylim=c(0,950))


###################################################
### code chunk number 27: fig1370431
###################################################
par(cex=0.7)
crplot(models[[2]][[selectedps]],"2",ylim=c(0,950),newplot=FALSE)


###################################################
### code chunk number 28: PSEA_RNAmixtures.Rnw:213-214
###################################################
sessionInfo()