mda14jx_week4-copy.ipynb

Kernel: R (R-Project)

options(jupyter.plot_mimetypes ='image/png')

Exercise 1

library(affy)

Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from ‘package:stats’:

    IQR, mad, xtabs

The following objects are masked from ‘package:base’:

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

?mas5

mas5() executes the command for MAS5.0 expression measure. This function converts an instance of AffyBatch into an instance of ExpressionSet using the implementation of Affymetrix's MAS5.0 expression measure. Using mas5(object, normalize = TRUE, sc = 500, analysis = "absolute", ...)

?rma

rma() executes the command for Robust Multi-Array Average expression measure. This function conversts an AffyBatch object into an ExpressionSet object using the RMA expression measure. Using rma(object, subset=NULL, verbose=TRUE, destructive=TRUE, normalize=TRUE, background=TRUE, bgversion=2, ...)

Exercise 2

setwd("~/Autumn2016/Week4/data_wk4")

getwd()

[1] "/projects/ddda6a8e-2bca-47f5-b1d6-79b2c48d0e30/Autumn2016/Week4/data_wk4"

Exercise 3

load("Dilution.RDA")

show(Dilution)

AffyBatch (storageMode: list)
assayData: 640 features, 640 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: 20A 20B 10A 10B
  varLabels: liver sn19 scanner
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2 

The data shows the size of the array is 640x640 (35221kb), cdf maps each gene that is in the array (12625 genes), and there are 4 samples.

Exercise 4

rma(Dilution)
eset_rma<-rma(Dilution)
show(eset_rma)

Background correcting
Normalizing
Calculating Expression

ExpressionSet (storageMode: lockedEnvironment)
assayData: 12625 features, 4 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: 20A 20B 10A 10B
  varLabels: liver sn19 scanner
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2 

Background correcting
Normalizing
Calculating Expression
ExpressionSet (storageMode: lockedEnvironment)
assayData: 12625 features, 4 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: 20A 20B 10A 10B
  varLabels: liver sn19 scanner
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2 

eset_mas5<-mas5(Dilution)

background correction: mas 
PM/MM correction : mas 
expression values: mas 
background correcting...done.
12625 ids to be processed
|                    |
|####################|

e_rma<-exprs(eset_rma)
show(e_rma)

e_mas5<-exprs(eset_mas5)
show(e_mas5)

par(mfrow=c(1,2))
d1<-density(e_rma)
plot(d1)
d2<-density(e_mas5)
plot(d2)

par(mfrow=c(1,2))
boxplot(e_rma)
boxplot(e_mas5)

log2e_mas5<-log2(e_mas5)
show(log2e_mas5)

par(mfrow=c(1,2))
log2d<-density(log2e_mas5)
plot(log2d)
boxplot(log2e_mas5)

Exercise 6

library(puma)

Error in library(puma): there is no package called ‘puma’
Traceback:
1. library(puma)
2. stop(txt, domain = NA)

setwd("~/Autumn2016/Week4/data_wk4")

load("Dilution.RDA")
show(Dilution)

eset_mmgmos<-mmgmos(Dilution)

e_puma<-exprs(eset_mmgmos)
show(e_puma)

par(mfrow=c(1,2))
pumad<-density(e_puma)
plot(pumad)
boxplot(e_puma)

e_puma[e_puma < 0]<- 0

par(mfrow=c(1,2))
pumad<-density(e_puma)
plot(pumad)
boxplot(e_puma)

setwd("~/Autumn2016/Week4/data_wk4")

load("affybatch.estrogen.RDA")
show(affybatch.estrogen)

pData(affybatch.estrogen) <- data.frame(
"estrogen"=c("absent","absent","present","present"
,"absent","absent","present","present")
, "time.h"=c("10","10","10","10","48","48","48","48")
, row.names=rownames(pData(affybatch.estrogen)))

eset_estrogen_mmgmos<-mmgmos(affybatch.estrogen)

e_estrogen_puma<-exprs(eset_estrogen_mmgmos)
show(e_estrogen_puma)