[BioC] Is there an expected coefficient of variation for the hybridization probe intensities among technical replicates?
James W. MacDonald
jmacdon at uw.edu
Thu Feb 27 18:00:53 CET 2014
Hi Matt,
I don't use xps, so will have to go about looking at this from an oligo
perspective. It just so happens I am working on a study using these same
arrays.
library(oligo)
dat <- read.celfiles(filenames = list.files("../CEL", pattern =
"[Cc][Ee][Ll]$", full.names = TRUE))
eset <- rma(dat)
con <- db(dat)
> dbGetQuery(con, "select * from type_dict;")
type type_id
1 1 main
2 2 control->affx
3 3 control->chip
4 4 control->bgp->antigenomic
5 5 control->bgp->genomic
6 6 normgene->exon
7 7 normgene->intron
8 8 rescue->FLmRNA->unmapped
9 9 control->affx->bac_spike
10 10 oligo_spike_in
11 11 r1_bac_spike_at
> table(dbGetQuery(con, "select type from featureSet inner join
core_mps using(fsetid);")[,1])
1 2 4 7 9
215001 18 23 5083 18
So I see control->affx and control->affx->bac_spike controls.
affx <- exprs(eset)[as.character(dbGetQuery(con, "select
featureSet.transcript_cluster_id from featureSet inner join core_mps
using(fsetid) where type='2';")[,1]),]
affx2 <- exprs(eset)[as.character(dbGetQuery(con, "select
featureSet.transcript_cluster_id from featureSet inner join core_mps
using(fsetid) where type='9';")[,1]),]
> z <-round(cbind(apply(affx, 1, function(x)
sd(x)/mean(x)),apply(affx2, 1, function(x) sd(x)/mean(x)),
deparse.level=0), 2)
> row.names(z) <- NULL
> z
[,1] [,2]
[1,] 0.02 0.02
[2,] 0.03 0.03
[3,] 0.02 0.02
[4,] 0.03 0.04
[5,] 0.02 0.02
[6,] 0.03 0.04
[7,] 0.01 0.01
[8,] 0.03 0.03
[9,] 0.01 0.01
[10,] 0.03 0.03
[11,] 0.01 0.01
[12,] 0.01 0.04
[13,] 0.01 0.01
[14,] 0.02 0.03
[15,] 0.01 0.00
[16,] 0.02 0.04
[17,] 0.01 0.00
[18,] 0.03 0.03
So much lower CVs than what you see. However, you should note that you
are not using logged data, which is fairly nonsensical in this situation
(e.g., computing a CV on highly skewed data is a recipe for biased results).
Best,
Jim
On 2/26/2014 6:02 PM, Thornton, Matthew wrote:
> Hello!
>
> I am digging through the control intensities for technical replicates obtained with an Affymetrix GeneChip Rat Gene 2.0 ST array.
>
> I have 3 technical replicates of the same sample. The kernel density plot overlay of the raw data shows differences in the intensity distributions.
>
> When I pull out the raw intensities with grepl("^AFFX", GeneName), I get the recorded intensities for each pixel mapped to a TranscriptClusterID corresponding to the particular AFFX control, some are 20 pixels others are 11. The groups of AFFX controls I get are the hybridization controls: AFFX-BioB, AFFX-BioC, AFFX-BioDn,
> AFFX-CreX; then the PolyA controls:AFFX-DapX, AFFX-LysX, AFFX-PheX, AFFX-ThrX, AFFX-TrpnX (?! not listed anywhere, is it a PolyA or hyb?), AFFX-r2-Bs-dap, AFFX-r2-Bs-lys, AFFX-r2-Bs-phe, and AFFX-r2-Bs-thr. Another set of PolyA controls: AFFX-r2-Ec-bioB, AFFX-r2-Ec-bioC, AFFX-r2-Ec-bioD, AFFX-r2-P1-cre. I also get the ERCC spike-in controls (which we did use). The maximum intensity from all controls is 25,069 grey levels the average is 1500 grey levels. What is a "good" average intensity?
>
> If I calculate the coefficient of variation across each row in my data set, which corresponds to each X,Y intensity measurement across my replicates for each probe. I see an average coefficient of variation at 22% just for the hybridization controls. Is this consistent with the experience of the list? What should be the expected coefficient of variation for the hybridization controls across technical replicates? the coefficient of variation is very similar for samples that show very low intensity values which are probably not expressed and for samples close to the average intensity. Even though I haven't done it, I don't think excluding probes based on a threshold with improve the average, because the numbers are uniformly similar.
>
> Here is my code:
> #!/usr/bin/Rscript --vanilla --slave
>
> # Load Libraries
> library(Biobase)
> library(xps)
>
> # Directory to store ROOT scheme files
> scmdir <- "/data/met/scmdir/"
> scheme_RaGene <- root.scheme(file.path(scmdir, "scheme_RaGene20stv1.root"))
>
> # Collect scheme information - Pulls out information from the Root scheme
> ann <- export(scheme_RaGene, treetype="ann", outfile="scheme_RaGene_ann.txt", as.dataframe=TRUE, verbose=FALSE)
> scm <- export(scheme_RaGene, treetype="scm", outfile="scheme_RaGene_scm.txt", as.dataframe=TRUE, verbose=FALSE)
>
> # more file locations
> datdir <- "/data/met/26Feb14_TR_CTRs/"
> celdir <- "/data/met/26Feb14_TR_CTRs/CEL/"
> celfiles <- c("CTR_TR1.CEL", "CTR_TR2.CEL", "CTR_TR3.CEL")
> celnames <- c("CTR_TR1", "CTR_TR2", "CTR_TR3")
> outdir <- "/data/met/26Feb14_TR_CTRs/"
>
> # Import CEL files and begin
> data_raw <- import.data(scheme_RaGene, "raw_data_CTR", filedir=outdir, celdir=celdir, celfiles=celfiles, celnames=celnames, verbose=TRUE)
>
> # Attach info
> data_raw <- attachMask(data_raw)
> data_raw <- attachInten(data_raw)
>
> tmp <- intensity(data_raw)
>
> write.csv(tmp, file="raw_data.csv", row.names=FALSE)
>
> # Save an image
> save.image(file="24Feb14_1.RData")
>
> ##
> ## Manipulation of "Raw Data"
> ##
>
> # plots
> png(file="Raw_Data_Density_Plot.png", width=600, height=600)
> par(mar=c(6,3,1,1));
> hist(data_raw, add.legend=TRUE)
> dev.off()
>
> # Get AFFX controls for assessment of PolyA, ERCC, and Hyb
> annot_raw <- merge(tmp, scm, by=c("X", "Y"), all.x=TRUE)
> annot_raw2 <- merge(annot_raw, ann, all.x=TRUE)
> raw_data <- annot_raw2[ order(annot_raw2[,3], annot_raw2[,2]), ]
> write.csv(raw_data, file="raw_data_annotated.csv", row.names=FALSE)
>
> # get the only handle for ERCC
> set.affx <- subset(raw_data, grepl("^AFFX", GeneName))
> write.csv(set.affx, file="affx_controls_raw.csv", row.names=FALSE)
>
> # I usually open the csv file with a spreadsheet program delete the columns that I don't need and calculate the coefficient of variation for each set.
>
> Any comments, advice, or information are greatly appreciated!!
>
> Matt
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
--
James W. MacDonald, M.S.
Biostatistician
University of Washington
Environmental and Occupational Health Sciences
4225 Roosevelt Way NE, # 100
Seattle WA 98105-6099
More information about the Bioconductor
mailing list