[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:58:48 CET 2014
Hi Matt,
Please don't take conversations off list.
On 2/27/2014 12:37 PM, Thornton, Matthew wrote:
> Hi Jim!
>
> Thank you for your reply, you have helped me so many times..
>
> So this data that I am pulling out is raw, no background correction or normalization. If we were to process the data with RMA or VSN, then that wouldn't give us the same information. I am trying to see how similar technical replicates as raw data should be before processing. This probably more of an Affymetrix technical question. I figured that people who use Bioconductor to process these chips in large numbers might have a better real world consensus. I am thinking that because the hybridization cocktail is a master-mix added exactly the same to each sample before hybridization that the coeffcient of variation should be very similar and that the intensity distribution should be practically identical. How realistic is it to have practically identical technical replicates?
Without background correction and normalization? Not very.
You have to think about what we are doing here. We have this silicone
wafer with a bunch of little short oligos attached, and we are splashing
some labeled cDNA (of various lengths) on, and letting it swish around
for a while, hoping that the cDNA will bind to complementary sequences
on the array. We then wash that off, and splash on some SAPE, let that
bind for a while and then wash off the excess. We then put it in a
scanner that takes a picture and digitizes the intensities.
You then want to take these digitized intensities and see how consistent
they are between arrays. But those intensities arise from a combination
of actual binding of cDNA to the oligo as well as biotin-streptavidin
binding, non-specific binding, quality of the various washes, pipetting
acumen of the tech, manufacturing inconsistencies, batch effects,
differences in how 'warmed up' the scanner is, hobgoblins of the mind
(ok, well maybe not that one), etc, etc.
And of all those things, only one has to do with what you are trying to
measure. We can't really quantify the contribution(s) of the other
things, but instead just try to make some assumptions and strip out the
technical variation as best we can. But if you don't even try to strip
out the technical variation, then you are looking at signal+noise and
expecting consistency.
Best,
Jim
>
> Thank you!!
>
> Matt
>
> ________________________________________
> From: James W. MacDonald [jmacdon at uw.edu]
> Sent: Thursday, February 27, 2014 9:00 AM
> To: Thornton, Matthew
> Cc: bioconductor at r-project.org
> Subject: Re: [BioC] Is there an expected coefficient of variation for the hybridization probe intensities among technical replicates?
>
> 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
>
--
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