Setting up aroma.affymetrix for analysis

Setting up aroma.affymetrix for analysis two directories : rawData and annotationData
from where aroma is launched

rawData/cancer_name/HuGene-1_0-st-v1/*.CEL
annotationData/chipTypes/HuGene-1_0-st-v1/*.cdf
annotationData/chipTypes/HuGene-1_0-st-v1/NetAffx/*.csv

RLE AND NUSE plots

Affy QC plots for exon arrays ...

One method of deciding whether or not an array is problematics from a quality standpoint is NUSE. The goal of NUSE is to identify any arrays which have elevated standard errors relative to other arrays in the dataset. This is done by standardizing the SE across arrays to have median 1 for each probeset. Our graphical tool consists of boxplots of these quantities for each array. A discordant boxplot indicates it is of poorer quality relative to the rest of the dataset. Instead of visually examining these quantities suitable numerical summaries such as the median and IQR NUSE could be used.

Another tool for making a decision about whether an array should be removed from subsequent analysis because of poor quality is RLE. These are the log-scale expression values relative to the median expression value computed on a probeset by probeset basis. a significantly different boxplot indicates problem.

aroma.affymetrix and other lots of packages allows to plot RLE and NUSE.

The NUSE is generally considered more sensitive than the RLE.

Bioconductor post on using arrayQualityMatrics with exon arrays

Looks like both simpleaffy and arrayQualityMetrics have problem with QCing Affy Exon 1.0 ST arrays.

Though following post does suggest a way to put custom CDF.


Hi Gard,

Sorry for the delay answering. I do not have much experience using
arrayQualityMetrics for Exon arrays, so I have talked with Crispin Miller
(simpleaffy package) about it and according to him "most of the Affymetrix
QC metrics for the 3' IVT arrays aren't directly applicable to the exon
arrays. They rely on MAS 5 and paired MM spots (neither of which are
applicable for exon arrays) and also make assumptions on 3'/5' ratios that
don't apply because the exon array chemistry is different."
I have now modified the package and version 2.4.3 of arrayQualityMetrics
should not perform the QC statistics from simpleaffy when "exon" is in the
cdfname.

Best wishes,
Audrey

> Hi.
>
> I am trying to get the arrayQualityMetrics package to run on a set of
CEL files from the Human Exon array from Affymetrix
>
> My problems begin when I want to run the arrayQualityMetrics function
and it gives the following error message :
>
> running R 2.9.2 and bioconductor version 2.4
>
> >library(affy)
> >library(simpleaffy)
> >ibrary(arrayQualityMetrics)
> >ecesbatch<-read.affybatch("H1.CEL", "H2.CEL", "H3.CEL", "H4.CEL",
> "H5.CEL", "H6.CEL", "H7.CEL", "H8.CEL", "H9.CEL", "H10.CEL",
> "H11.CEL", "H12.CEL", "H13.CEL", "H14.CEL", "H15.CEL", "H16.CEL")
>
> ## attach cdf to expr set
> ecesbatch cdfName <- "exon.pmcdf" ## this is a cdf file from the XMAP
website
>
> #Check the name is correct for the cdf file (unneccessary)
> > cdfname <- cleancdfname(cdfName(ecesbatch))
> > cdfname
> [1] "exon.pmcdf"
>
> >arrayQualityMetrics(expressionset = ecesbatch,outdir =
> "output",force = TRUE,do.logtransform = TRUE)
>
> This cmd runs for a very long time and generates a bunch of .pdfs and
.pngs and an empty QCReport.html file.
> And R says there is an error sonce the arrayQualityMetrics package does
not know the QCparameters of this chip.
>
>
> I have found an instruction from C. Miller (one of the persons behind
simpleaffy) about how use the three functions provided by the
> simpleaffy package, or to make the needed .qcdf file:
> I need alpha values (that is okay) and I need control and spike
> probeIDs.
>
> I am using the Human Exon array 1.0 from affymetrix, and I do not know
what to fill in in the .qcdf file,
> anyone who knows how to get by this problem?
> Trying to get the probenames to set the values I ran into another problem..
>
> > prbs <- ls(cdfname)
> Error in as.environment(pos) :
> no item called "exon.pmcdf" on the search list
> >
>
> crashes like this shown here.
>
>
> Please if anyone knows or has an idea, basicly what I need is
> the .qcdf file for the HUman Exon array from Affy.
> Best regards
> Gard
>
> #################################
> Gard Thomassen
> Ph.D student CMBN, Rikshospitalet, Oslo
> Bioinformatician, Radiumhospitalet, Oslo
> Norway
> Email : gardt@...
> Office: + 47 22781736
> Phone +47 93674926

Quality control of Affymetrix arrays

I am playing around with different Bioconductor packages for QC on affy exon arrays.

Here is a nice introduction to quality assessment and processing.

http://www.bioconductor.org/workshops/2009/GenentechNov2009/Module2/module2-affy-preprocess.pdf


Here are some files one should have:

• CEL : contain one observation per spot
• CDF : map from spot locations to probeset and ultimately to the identity of the
gene being probed
• Bioconductor annotation" packages map from probe sets to gene and other
annotations.
• Tab-delimited, database, or other les provide phenotypic information.



Some packages are

arrayQualityMetrics
SimpleAffy
yaqcaffy
estrogen package vignette also has some QC. [>openVignette("estrogen")]

Before we start normnalization process here are some quality matric that Affymetrix advises

1)Average background : should be similar for all chips
2) Scale Factor: should be within 3 fold
3) # of genes called present : For similar samples - number should be similar. May
be different for different tissue types.
4) 3' to 5' ratio of GAPDH and beta-actin : should be close to one up to 3 is fine.
1.25 is what "simpleaffy" recommands.
5)Value for spike in transcripts: present in atleast 70% of arrays
please see
http://bioconductor.org/packages/2.5/bioc/vignettes/simpleaffy/inst/doc/QCandSimpleaffy.pdf

for more info.

======

We are interested in looking at two different aspects : Per slide aspects and Between slide aspects. Per slide aspects are - intensity dependence of ratios and spatial effects on the array. This can be done by looking at MA plots and a false image of chip. Between slide aspects are Homogeneity, outlier samples and biological meanings. This can be done by Boxplots, density plots, Heatmap and PCA. Other plots include Variance-mean dependency, GC content and probe mapping studies. Other Affy only plots include NUSE, RLE, RNA degradation, QC stats, PM/MM. Finally, one should be able to identify outliers.

The image function allows us to look at the spatial distribution of the intensities on a chip.

Another way to visualize what is going on on a chip is to look at the histogram of
its intensity distribution. Because of the large dynamical range (O(104)), it is useful to look at the log-transformed values

To compare the intensity distribution across several chips, we can look at the boxplots, both of the raw intensities and the normalized probe set values

The scatterplot is a visualization that is useful for assessing the variation (or
reproducibility, depending on how you look at it) between chips. We can look at all probes, the perfect match probes only, the mismatch probes only, and of course also at the normalized, probe-set-summarized data

Differences between arrays in the shape or center of the distribution often highlight the need for normalization.

The MA plot is a rotated version of a scatter plot. The
rotation helps to detect patterns as deviations from horizontal,
rather than diagonal.

• Instead of ploting two vectors Y2;j versus Y1;j , we plot
Mj = Y2;j - Y1;j versus Aj = (Y2;j + Y1;j)=2.
• if Y1 and Y2 are logarithmic expression values, then
{ Mj represents fold change for gene j
{ Aj represents average log intensity for gene j.

Using Plotrix to create nice depletion/enrichment images

setwd ("/Users/shah/insulator_data/insulator_combination_numbers/clustering_visulaization")

# Read Table

insulatorFile <- read.table ("simulated_1fdr_no_250_pvalue_for_figure.txt", sep="\t")
presMatrix <- t(as.matrix (insulatorFile [,c (6,7,8,9,10,11)]))
#color2D.matplot (presMatrix, redrange=c(0.9,0), greenrange=c(0.9,0), bluerange=c(0.9,0), xlab = "presence", ylab ="insulator combination", axes=FALSE)

pv1 <- insulatorFile[,3]
changeVal <-function (x) { if (x <= 1e-16) { x <- 1e-16} else { x <- x}}

pv3 <- lapply (pv1, changeVal)
pv4 <- -1*log10(as.numeric(pv3))

direction <- insulatorFile[,4]

pv5 <- t(as.matrix(pv4*direction))

# Generate the colors

cellcol<-matrix(rep("#000000",63),nrow=63)
cellcol[pv5<0]<-color.scale(pv5[pv5<0], c(1,0),c(1,0),c(0,0))
cellcol[pv5>0]<-color.scale(pv5[pv5>0], c(0,0),c(0,0),c(0,1))

# Generate Legend (this one is yellow - black - blue)

legval<-seq(min(pv5),max(pv5),length.out=32)
legcol<-rep("#000000",32)

legcol[legval<0]<-color.scale(legval[legval<0], c(1,0),c(1,0), c(0,0))
legcol[legval>0]<-color.scale(legval[legval>0], c(0,0),c(0,0),c(0,1))


#color2D.matplot(pv5,cellcolors=cellcol,border=NA, axes=FALSE)
#color.legend(0,0,6,-4,round(c(min(pv5),0,max(pv5)),1),rect.col=legcol)
color2D.matplot(t(as.matrix(legcol)),cellcolors=legcol,border=NA, axes=FALSE)

xMAN commandline

Here is the xMAN how to from Tao Liu


Text formatted bpmap may not be processed by MAT directly. You need to convert your text bpmap to a so-called 'plainseq' format. Please type 'xMAN' without any argument to see the usage for detail.

And here is the commandline used to generate bpmap for fruitfly:

xMAN -f bpmap -i Dm_tiling2_MR_v01.bpmap -s Dm_tiling2_MR_v01.stderr -o Dm_tiling2_MR_v01_dm3 --nro=Dm_tiling2_MR_v01_dm3_all.NR --maxSeqCopy=10 --ProbeResolution=35 dm3/all/*.fa

The input file is a binary bpmap downloaded from affy website. And I put all the genome sequence files in dm3/all/ directory. Please remember to add '-nro' argument. After the run, the file 'Dm_tiling2_MR_v01_dm3_all.NR.bpmap' can be used by MAT.

Generating a Affymatrix "tpmap" file for TiMAT

R Code that uses the Bioconductor module Affxparser.

Note: The code has a bug and in normal circumstances returns an error

> writeTpmap(filename="Dm_tiling2_MR_v01.tpmap", bpmaplist=bpmapr)
> Error in as.vector(y) : argument "y" is missing, with no default

The "y" can be fixed as follows

> y <- readBpmap("Dm_tiling2_MR_v01.bpmap")
> writeTpmap(filename="Dm_tiling2_MR_v01.tpmap", bpmaplist='y', verbose = as.integer(0))

This seems to work but doesn't print anything into the tpmap file due to some "sequence skipping".

I will try the Affymatrix windows based program next.