Question: HTqPCR
0
7.2 years ago by
Heidi Dvinge40
Heidi Dvinge40 wrote:
qpcr go glad htqpcr • 1.2k views
modified 7.2 years ago by Heidi Dvinge2.0k • written 7.2 years ago by Heidi Dvinge40
0
7.2 years ago by
Heidi Dvinge40
Heidi Dvinge40 wrote:
0
7.2 years ago by
Heidi Dvinge40
Heidi Dvinge40 wrote:
0
7.2 years ago by
Heidi Dvinge2.0k
Heidi Dvinge2.0k wrote:
Hi Simon, thanks for your email, and sorry if the example files are shoddy. I'm travelling the next couple of days, but will have a look as soon as I'm back. Off the top of my head, I seem to remember that for some of the BioMark files I originally had access to, the sample names were in slightly different formats. Therefore, as a temporary measure I just ignored sample names in the file, in which case the sample names can be added with the 'samples' parameter to readCtData(). Or alternative, added with sampleNames(object) <- c("",...) later. Obviously, that's not optimal though, and I should fix that for future versions. In the case of the example file, a (somewhat messy) workaround may be to say: > exPath <- system.file("exData", package = "HTqPCR") > temp <- read.csv(file.path(exPath, "BioMark_sample.csv"), as.is=TRUE, skip=11) > raw1 <- readCtData(files = "BioMark_sample.csv", path = exPath, format = "BioMark", n.features = 48, n.data = 48, samples=temp$Name[seq(1,nrow(temp), 48)]) > head(sampleNames(raw1)) [1] "no preamp" "preamp neat" "no preamp.1" "preamp neat.1" "no preamp.2" [6] "preamp neat.2" Apart from the missing sample names, that what exactly are the problems you're seeing in importing your own data into a qPCRset object? HTH \Heidi > Hi Heidi, > I'm working on getting your package routinely going for for our biomark > data, and it looks really good. However, I've been trying to use your > example Biomark file in the package as an example, and I'm running into a > few problems. Firstly, its pretty normal to have uneven numbers of samples > per plate. We never do a single sample per plate! Remember, on the small > plates, there are 2304 PCR reactions, and on the large plates, 9216. Due > to the large volume of samples which can be run, we typically do multiple > different samples with the same suite of genes per plate, with varying > replicates. However, it can be tough to make things symetrical due to > running of controls (no template controls, etc). We also strive to have > the samples fields present in the CSV files as it makes things easier for > our existing analysis. > > In the file you supply, the samples do not appear to be recognized by your > scripts. If you do > pData(raw1) > > on the newly imported Biomark file you supply, you just get a generic list > of sample ID's for each of the 48 chambers, and the > > raw1 <- readCtData(files = "BioMark_sample2.csv",path = exPath, format = > "BioMark", n.features = 48,n.data = 48) > > command does not appear to import the sample names associated with each of > the genes which are as follows: > > > no preamp > preamp neat > preamp 1:10 > preamp 1:100 > preamp 1:1000 > > In the biomark file, there are also a number of samples which do not have > any identifiers (blank), so I filled these in as "Blank". Further, the > number of assays being run for each of these samples is quite variable, > ranging form a low of 48 assays, up to 1152 assays. > > I'm having trouble using the examples you provided to reformat the data so > that I can use your excellent tools. I'm also puzzled why the sample names > within the file are not being imported along with the Ct and gene data. > > I'd really like to be able to try various normalizing methods, but I cant > get there till we can import the data so it makes sense. > > best > > Simon. > > Simon Melov Ph.D. > Associate Professor & > Director of Genomics > Buck Institute for Research on Aging > 8001 Redwood Blvd > Novato, CA 94945 > > Office: 415 209 2068 > Cell: 415 827 4979 > Fax: 415 209 9920 > > > > > > > ADD COMMENTlink written 7.2 years ago by Heidi Dvinge2.0k Answer: HTqPCR 0 7.2 years ago by Heidi Dvinge40 Heidi Dvinge40 wrote: Hi Silvia, On 20 Jun 2012, at 19:25, Silvia Halim wrote: > Hi Heidi, > > I've tried sample.reps and feature.reps parameters on plotCtVariation, however, I couldn't make it work. It keeps saying 'figure margins too large'. This isn't related to HTqPCR as such, but rather a general R statement saying that you're trying to plot something that's too large for your current plotting device. Most devices (e.g. quartz(), X11()) will let you set a 'width' and 'height' parameter, or you can perhaps plot directly to a file, where you can set the dimensions as required. Just for the record, plotting a 96x96 object works for me. What it the output of the following on your computer: # Load 48x48 example data exPath <- system.file("exData", package = "HTqPCR") raw <- readCtData(files = "BioMark_sample.csv", path = exPath, format = "BioMark", n.features = 48, n.data = 48) # Create 'pseudo' 96x96 data set tmp <- cbind(raw, raw) raw96 <- rbind(tmp, tmp) # Plot plotCtVariation(tmp) > Actually how should I set either of these 2 parameters? > In your example in the help page, your command is plotCtVariation(qPCRraw[1:40,], sample.reps=rep(1:2,3)). Do you specify 'rep(1:2,3)' for sample.reps because each feature is replicated exactly twice only and you want to plot only 3 samples? > Typing rep(1:2, 3) in the console will give you a vector c(1,2,1,2,1,2), which in this case indicates that the 6 samples in my object falls into 2 different groups. Here they're just named '1' and '2' for the sake of ease, but I might as well have said e.g. sample.reps=c("control", "treatment", "control", "treatment", "control", "treatment") or whatever applies to a given experiment. The same goes for feature.reps. In the example, feature.reps=paste("test", rep(1:96, each=4)) simply refers to that on a 384-well assay there are 96 individual features, and each feature is present 4 times after each other. C.f. > paste("test", rep(1:96, each=4))[1:10] [1] "test 1" "test 1" "test 1" "test 1" "test 2" "test 2" "test 2" "test 2" "test 3" [10] "test 3" Typically, feature.reps would simply be the featureNames of your object given that each feature name appears multiple times. This is also the default behaviour of plotCtVariation. HTH \Heidi > And how should I interpret ' feature.reps=paste("test", rep(1:96, each=4))' in the other example of plotCtVariation? > > Thanks, > Silvia > > -----Original Message----- > From: Heidi Dvinge > Sent: 20 June 2012 9:46 AM > To: Silvia Halim > Cc: bioconductor at r-project.org > Subject: Re: HTqPCR > > Hi SIlvia, > > On 19 Jun 2012, at 18:48, Silvia Halim wrote: > >> Hi Heidi, >> >> Thanks for your tips. I figured I could probably use plotCtVariation(). I am able to use the function to plot variation across samples but how can I use it to plot across genes (or features)? >> I tried following commands: >> plotCtVariation(temp[,1:10], variation = "sd", log = TRUE, main = "SD >> of replicated features", col = "lightgrey") >> plotCtVariation(temp[1:10,], variation = "sd", log = TRUE, main = "SD >> of replicated features", col = "lightgrey") > > Here you're just subsetting your qPCRset object before plotting, but you're not changing the actual plots. > >> There's a difference in the plots but both plots give me same labels on x-axis, i.e. sample names, though I was expecting the second command would give me gene names on x-axis label. >> > If you look at the plotCtVariation help files (especially the 'Examples' and 'Details' section), the parameters sample.reps and feature.reps controls whether you plot the variation for each gene across or within samples. In order to get gene names, you have to set sample.reps, to indicate which samples are replicates of each other. > > Per default, the function calculates the variation between replicated features within each of your samples, and plots the distribution (boxplot) of this variation for each sample. If you want to check individual features or samples more specifically, you ahve to use type="detail" and possibly add.featurenames=TRUE. There are some examples included in the plotCtVariation help file. > >> Also, the manual says we can exclude unreliable or undetermined data by setting the Ct values to NA using filterCategory. I am wondering how I can get rid of NA data from the plate. I also cannot exclude this kind of data or those having 'Failed' flags the very first time before reading in the input as a qPCRset object because the input has to be something like 48 x 48 or 96 x 96. >> > The question is, why do you want to remove the NA values? If you just leave them as NAs, then they're ignored during e.g. the calculation of differential expression and for most plotting purposes. > > You can't remove them as such, since (as you note), the object has to be in a certain features x samples format. If you want, you can replace them though, if you e.g. want to set all NA values to Ct=40: exprs(temp)[is.na(exprs(temp))] <- 40. But beware, because in that case the value '40' will be include into all numerical calculations, which may not be what you want. > > HTH > \Heidi > >> Many thanks, >> Silvia >> >> -----Original Message----- >> From: Heidi Dvinge >> Sent: 19 June 2012 2:05 PM >> To: Silvia Halim >> Cc: bioconductor at r-project.org >> Subject: Re: HTqPCR >> >> Hi Silvia, >> >> On 18 Jun 2012, at 17:51, Silvia Halim wrote: >> >>> Hi Heidi, >>> >>> The function breaks at plotCtReps. >>>> traceback() >>> 1: plotCtReps(temp, card = 2, percent = 20, xlim = c(0, 100), ylim = c(0, >>> 100)) >>> >>> You've pointed out the problem about the duplicates as I have 3 replicates on my assay. I got confused reading the manual as it says plotCtReps can be used for a sample containing duplicate measurements (which I thought to be 2 or more measurements). >>>> table(table(featureNames(temp))) >>> >>> 3 6 >>> 30 1 >>> >> If you try running the examples for plotCtReps, you'll see that the function directly plots two replicates of a feature against each other on the (x,y) axis. 3D (x,y,z) plots aren't implemented, so features that are replicated 3 times can't be plotted. I'll try to clarify the text for the function. >> >> Perhaps something like plotCtVariation() will give you what you're after? If you only want to visually inspect your data, then grep("plot", ls("package:HTqPCR"), value=TRUE) will list all the plotting functions available in HTqPCR. >> >> HTH >> \Heidi >> >>> Btw there's no NA in my data. >>>> sumis.na(temp)) >>> [1] 0 >>> Warning message: >>> In is.na(temp) : is.na() applied to non-(list or vector) of type 'S4' >>>> >>> >>> Thanks, >>> Silvia >>> >>> -----Original Message----- >>> From: Heidi Dvinge >>> Sent: 15 June 2012 9:06 PM >>> To: Silvia Halim >>> Cc: bioconductor at r-project.org >>> Subject: Re: HTqPCR >>> >>> Hi Silvia, >>> >>> On 15 Jun 2012, at 18:45, Silvia Halim wrote: >>> >>>> Hi Heidi, >>>> >>>> I ran into below problem when using plotCtReps. >>>> >>>>> plotCtReps(temp, card = 1, percent = 20, xlim = c(0,50), ylim = >>>>> c(0,50)) >>>> Error in split.data[[s]] : subscript out of bounds In addition: >>>> Warning messages: >>>> 1: In min(x, na.rm = na.rm) : >>>> no non-missing arguments to min; returning Inf >>>> 2: In max(x, na.rm = na.rm) : >>>> no non-missing arguments to max; returning -Inf >>>>> plotCtReps(temp, card = 1, percent = 20, xlim = c(0,50), ylim = >>>>> c(0,50)) >>>> Error in split.data[[s]] : subscript out of bounds In addition: >>>> Warning messages: >>>> 1: In min(x, na.rm = na.rm) : >>>> no non-missing arguments to min; returning Inf >>>> 2: In max(x, na.rm = na.rm) : >>>> no non-missing arguments to max; returning -Inf >>>>> plotCtReps(temp, card = 2, percent = 20, xlim = c(0,100), ylim = >>>>> c(0,100)) >>>> Error in split.data[[s]] : subscript out of bounds In addition: >>>> Warning messages: >>>> 1: In min(x, na.rm = na.rm) : >>>> no non-missing arguments to min; returning Inf >>>> 2: In max(x, na.rm = na.rm) : >>>> no non-missing arguments to max; returning -Inf >>> >>> What's the output from traceback(), i.e. exactly where does the function break? >>>> >>> A couple of things you can try: >>> >>> - plotCtReps is meant to be used in cases where there are exactly 2 replicates of the features on your assay. Is this the case? For example, with the data below there are 190 features that will be plotted, and 1 that will be skipped: >>>> data(qPCRraw) >>>> table(table(featureNames(qPCRraw))) >>> 2 4 >>> 190 1 >>> >>> - are there any NAs in your data? E.g. sumis.na(qPCRraw))>0. >>> >>> HTH >>> \Heidi >>> >>>> Here is how 'temp' looks like >>>>> temp >>>> An object of class "qPCRset" >>>> Size: 96 features, 96 samples >>>> Feature types: Reference, Test >>>> Feature names: b-Actin b-Actin b-Actin ... >>>> Feature classes: >>>> Feature categories: OK >>>> Sample names: NTC_4 PMPT352 NTC_3 ... >>>> >>>> Do you know why it is complaining about split.data? >>>> >>>> Thanks, >>>> Silvia >>>> >>>> -----Original Message----- >>>> From: Heidi Dvinge >>>> Sent: 11 June 2012 6:11 PM >>>> To: Silvia Halim >>>> Subject: Re: HTqPCR >>>> >>>> Ok, so you already have a 96 by 96 matrix, so you don't need changeCtLayout. >>>> Good luck with the rest, and let me know if you encounter any problems. >>>> >>>> On 11 Jun 2012, at 19:05, Silvia Halim wrote: >>>> >>>>> Hi Heidi, >>>>> >>>>> Thank you for your clarification. >>>>> >>>>> Btw this is how it looks like when I type 'temp' >>>>>> temp >>>>> An object of class "qPCRset" >>>>> Size: 96 features, 96 samples >>>>> Feature types: Reference, Test >>>>> Feature names: b-Actin b-Actin b-Actin ... >>>>> Feature classes: >>>>> Feature categories: OK >>>>> Sample names: NTC_4 PMPT352 NTC_3 ... >>>>> >>>>> Cheers, >>>>> Silvia >>>>> >>>>> -----Original Message----- >>>>> From: Heidi Dvinge >>>>> Sent: 08 June 2012 7:12 PM >>>>> To: Silvia Halim >>>>> Subject: Re: HTqPCR >>>>> >>>>> Hi Silvia, >>>>> >>>>> what are the dimensions of the "temp" object that you read in? I.e. >>>>> what does it look like if you just type >>>>>> temp >>>>> >>>>> If you read in the data with n.features=96 and n.data=96, then you should already have an object with 96 rows and 96 columns, in which case you don't need to change the layout. >>>>> >>>>> Best, >>>>> \Heidi >>>>> >>>>> On 8 Jun 2012, at 19:13, Silvia Halim wrote: >>>>> >>>>>> Hi Heidi, >>>>>> >>>>>> I finally have time to try out your HTqPCR bioconductor package again and I was trying to use 'changeCtLayout' function. However, I got following error message: >>>>>> >>>>>>> qPCRnew <- changeCtLayout(temp, sample.order = sample_order) >>>>>> Error in data.frame(..., check.names = FALSE) : >>>>>> arguments imply differing number of rows: 0, 96 In addition: >>>>>> Warning >>>>>> message: >>>>>> In split.default(x = seq_len(nrow(x)), f = f, drop = drop, ...) : >>>>>> data length is not a multiple of split variable >>>>>> >>>>>> The commands that I run are as follows: >>>>>>> temp <- readCtData("110614 BENIGN_1 DATA 96X96.csv", path = >>>>>>> getwd(), n.features = 96, n.data=96, flag = 9, feature = 5, type= >>>>>>> 6, Ct = 7, position = 1, skip = 12, sep = ",") sample_order <- >>>>>>> rep(sampleNames(temp), each = 96) qPCRnew <- changeCtLayout(temp, >>>>>>> sample.order = sample_order) >>>>>> >>>>>> I've tried to follow what's written in changeCtLayout function description. Can you please advise what went wrong? >>>>>> >>>>>> Thanks, >>>>>> Silvia >>>>>> >>>>>> -----Original Message----- >>>>>> From: Heidi Dvinge >>>>>> Sent: 29 April 2012 8:18 PM >>>>>> To: Silvia Halim >>>>>> Subject: Re: HTqPCR >>>>>> >>>>>> HI Silvia, >>>>>> >>>>>> I'm glad you got it working. Depending on what you're supposed to do with the data, you may need to tweak some functions slightly, as you mention. Let me know if you run into any more trouble. >>>>>> >>>>>> Cheers >>>>>> \Heidi >>>>>> >>>>>> On 26 Apr 2012, at 18:37, Silvia Halim wrote: >>>>>> >>>>>>> Hi Heidi, >>>>>>> >>>>>>> Thanks for the help! It's working for me now. Right now I'm figuring it out how I can use the functions that you described in the vignette. I might have to tweak the parameters for using the functions on Fluidigm data. >>>>>>> >>>>>>> Cheers, >>>>>>> Silvia >>>>>>> >>>>>>> -----Original Message----- >>>>>>> From: Heidi Dvinge >>>>>>> Sent: 25 April 2012 8:56 AM >>>>>>> To: Silvia Halim >>>>>>> Subject: Re: HTqPCR >>>>>>> >>>>>>> Hiya, >>>>>>> >>>>>>> sorry, I only just now realised that you'd attached a file. When I saved as csv, the following command worked: >>>>>>> >>>>>>>> raw <- readCtData("110614 BENIGN_1 DATA 96x96.csv", >>>>>>>> format="BioMark", >>>>>>>> n.features=96*96) raw >>>>>>> An object of class "qPCRset" >>>>>>> Size: 9216 features, 1 samples >>>>>>> Feature types: >>>>>>> Feature names: b-Actin b-Actin b-Actin ... >>>>>>> Feature classes: >>>>>>> Feature categories: OK >>>>>>> Sample names: 110614 BENIGN_1 DATA 96x96 ... >>>>>>> >>>>>>> The data isn't transformed into a 96x96 format immediately though (in case you read in multiple arrays, and want to normalise them independently). If you want to change this, you can use changeCtLayout(). Alternatively you can say: >>>>>>> >>>>>>>> raw <- readCtData("110614 BENIGN_1 DATA 96x96.csv", >>>>>>>> format="BioMark", n.features=96, n.data=96) raw >>>>>>> An object of class "qPCRset" >>>>>>> Size: 96 features, 96 samples >>>>>>> Feature types: >>>>>>> Feature names: b-Actin b-Actin b-Actin ... >>>>>>> Feature classes: >>>>>>> Feature categories: OK >>>>>>> Sample names: Sample1 Sample2 Sample3 ... >>>>>>>> plotCtArray(raw) >>>>>>> >>>>>>> HTH >>>>>>> \Heidi >>>>>>> >>>>>>> On 24 Apr 2012, at 17:55, Silvia Halim wrote: >>>>>>> >>>>>>>> Hi Heidi, >>>>>>>> >>>>>>>> I have some problems updating R on lustre. Therefore, I chose to run HTqPCR on my desktop for the moment. >>>>>>>> >>>>>>>> Reading in your sample file looks fine, however, reading in the >>>>>>>> file that I showed you just now gave me below error message. >>>>>>>> (The file is as attached) >>>>>>>> >>>>>>>>> temp <- readCtData("110614 BENIGN_1 DATA 96x96.xlsx", path = >>>>>>>>> getwd() , n.features = 96*96, flag = 9, feature = 5, type= 6, >>>>>>>>> Ct = 7,position = 1, skip = 12, sep = ",") >>>>>>>> Error in read.table(file = file, header = header, sep = sep, quote = quote, : >>>>>>>> no lines available in input >>>>>>>> In addition: Warning message: >>>>>>>> In readLines(file, skip) : >>>>>>>> incomplete final line found on 'C:/Users/halim01/Documents/20110627_RossAdamsH_DN_Fluid/110614 BENIGN_1 DATA 96x96.xlsx' >>>>>>>>> sessionInfo() >>>>>>>> R version 2.14.0 (2011-10-31) >>>>>>>> Platform: x86_64-pc-mingw32/x64 (64-bit) >>>>>>>> >>>>>>>> locale: >>>>>>>> [1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United Kingdom.1252 LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C LC_TIME=English_United Kingdom.1252 >>>>>>>> >>>>>>>> attached base packages: >>>>>>>> [1] stats graphics grDevices utils datasets methods base >>>>>>>> >>>>>>>> other attached packages: >>>>>>>> [1] Biostrings_2.22.0 IRanges_1.12.6 BiocInstaller_1.2.1 marray_1.32.0 HTqPCR_1.8.0 limma_3.10.3 RColorBrewer_1.0-5 Biobase_2.14.0 gdata_2.8.2 >>>>>>>> >>>>>>>> loaded via a namespace (and not attached): >>>>>>>> [1] affy_1.32.1 affyio_1.22.0 gplots_2.10.1 gtools_2.6.2 preprocessCore_1.16.0 tools_2.14.0 zlibbioc_1.0.1 >>>>>>>>> >>>>>>>> >>>>>>>> I did a quick check on the file and it only has 9228 lines including 12 header lines which I had skipped when reading in the file. Do you know what could possibly go wrong? >>>>>>>> >>>>>>>> Cheers, >>>>>>>> Silvia >>>>>>>> >>>>>>>> -----Original Message----- >>>>>>>> From: Heidi Dvinge >>>>>>>> Sent: 24 April 2012 5:09 PM >>>>>>>> To: Silvia Halim >>>>>>>> Subject: Re: HTqPCR >>>>>>>> >>>>>>>> Hm, that looks like it may be x11 acting up. I often have similar issues when I work on a remote server. >>>>>>>> >>>>>>>> Actually, the processing of Fluidigm files is very computationally light. So you can easily do it on your desktop, if you can't update on lustre. >>>>>>>> >>>>>>>> I can also email you and older version of the vignette if you want to have a look. However, in HTqPCR 1.2.0 I don't even think I had a dedicated function for plotting the Fluidigm assays yet (the plotCtArray shown in the vignette). >>>>>>>> >>>>>>>> Cheers >>>>>>>> \Heidi >>>>>>>> >>>>>>>> On 24 Apr 2012, at 16:39, Silvia Halim wrote: >>>>>>>> >>>>>>>>> Hi Heidi, >>>>>>>>> >>>>>>>>> This is what I got when accessing the vignette. >>>>>>>>> >>>>>>>>>> openVignette(package="HTqPCR") >>>>>>>>> Please select a vignette: >>>>>>>>> >>>>>>>>> 1: HTqPCR - qPCR analysis in R >>>>>>>>> >>>>>>>>> Selection: 1 >>>>>>>>> Opening >>>>>>>>> /home/mib-cri/local/lib64/R/library/HTqPCR/doc/HTqPCR.pdf >>>>>>>>>> xprop: unable to open display '' >>>>>>>>> /usr/local/bin/xdg-open: line 370: firefox: command not found >>>>>>>>> /usr/local/bin/xdg-open: line 370: mozilla: command not found >>>>>>>>> /usr/local/bin/xdg-open: line 370: netscape: command not found >>>>>>>>> xdg-open: no method available for opening '/home/mib- cri/local/lib64/R/library/HTqPCR/doc/HTqPCR.pdf' >>>>>>>>> >>>>>>>>> Sorry for the confusion, you are right that I was looking at a newer version of HTqPCR than the one installed on lustre. I think that's because I have different installations of HTqPCR on lustre and on my desktop. If I can update the one on lustre, I'll go ahead with the update. >>>>>>>>> >>>>>>>>> Thank you, >>>>>>>>> Silvia >>>>>>>>> >>>>>>>>> -----Original Message----- >>>>>>>>> From: Heidi Dvinge >>>>>>>>> Sent: 24 April 2012 4:28 PM >>>>>>>>> To: Silvia Halim >>>>>>>>> Subject: Re: HTqPCR >>>>>>>>> >>>>>>>>> Ah, right, it looks like you have an older version of R, and therefore also HTqPCR. >>>>>>>>> >>>>>>>>> The most current release version is 1.10.0. In that version, readCtData() was modified to accept different types of input data, including from Fluidigm. Before that, this sort of data had to be read in 'manually'. >>>>>>>>> >>>>>>>>> I guess the vignette that you were looking at comes from a >>>>>>>>> version of HTqPCR that's newer than the one you have installed? >>>>>>>>> If you access the vignette corresponding to your HTqPCR version >>>>>>>>> via >>>>>>>>>> openVignette(package="HTqPCR") >>>>>>>>> what do you get then? >>>>>>>>> >>>>>>>>> If you get an older version, then depending on how old it is, there may be a section towards the end giving an example of how to process Fluidigm data more 'manually'. If not, an update may be your best bet. >>>>>>>>> >>>>>>>>> Cheers >>>>>>>>> \Heidi >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> On 24 Apr 2012, at 16:15, Silvia Halim wrote: >>>>>>>>> >>>>>>>>>> Hi Heidi, >>>>>>>>>> >>>>>>>>>> Thanks for looking into the matter. Below is the output of my >>>>>>>>>> sessionInfo() >>>>>>>>>> >>>>>>>>>>> sessionInfo() >>>>>>>>>> R version 2.13.0 (2011-04-13) >>>>>>>>>> Platform: x86_64-unknown-linux-gnu (64-bit) >>>>>>>>>> >>>>>>>>>> locale: >>>>>>>>>> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C >>>>>>>>>> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 >>>>>>>>>> [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 >>>>>>>>>> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C >>>>>>>>>> [9] LC_ADDRESS=C LC_TELEPHONE=C >>>>>>>>>> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C >>>>>>>>>> >>>>>>>>>> attached base packages: >>>>>>>>>> [1] stats graphics grDevices utils datasets methods base >>>>>>>>>> >>>>>>>>>> other attached packages: >>>>>>>>>> [1] marray_1.26.0 Biostrings_2.20.1 IRanges_1.10.3 HTqPCR_1.2.0 >>>>>>>>>> [5] limma_3.6.9 RColorBrewer_1.0-2 Biobase_2.12.1 gdata_2.8.0 >>>>>>>>>> >>>>>>>>>> loaded via a namespace (and not attached): >>>>>>>>>> [1] affy_1.26.1 affyio_1.20.0 gplots_2.8.0 >>>>>>>>>> [4] gtools_2.6.2 preprocessCore_1.14.0 >>>>>>>>>>> >>>>>>>>>> >>>>>>>>>> Cheers, >>>>>>>>>> Silvia >>>>>>>>>> >>>>>>>>>> -----Original Message----- >>>>>>>>>> From: Heidi Dvinge >>>>>>>>>> Sent: 24 April 2012 4:07 PM >>>>>>>>>> To: Silvia Halim >>>>>>>>>> Subject: HTqPCR >>>>>>>>>> >>>>>>>>>> Hi Silvia, >>>>>>>>>> >>>>>>>>>> I just tested the read fluidigm from the vignette, and it works on both my mac and a single unix system that I've tested. Although from the errors you were getting, it seemed like the headers weren't been read correctly/at all. >>>>>>>>>> >>>>>>>>>> Would you mind sending me the output of your sessionInfo(), so I can compare which package versions we have? >>>>>>>>>> >>>>>>>>>> Best, >>>>>>>>>> \Heidi >>>>>>>>>> >>>>>>>>>>> sessionInfo() >>>>>>>>>> R version 2.15.0 (2012-03-30) >>>>>>>>>> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) >>>>>>>>>> >>>>>>>>>> locale: >>>>>>>>>> [1] >>>>>>>>>> en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8 >>>>>>>>>> >>>>>>>>>> attached base packages: >>>>>>>>>> [1] tools stats graphics grDevices utils datasets methods base >>>>>>>>>> >>>>>>>>>> other attached packages: >>>>>>>>>> [1] HTqPCR_1.10.0 limma_3.12.0 RColorBrewer_1.0-5 Biobase_2.16.0 >>>>>>>>>> [5] BiocGenerics_0.2.0 >>>>>>>>>> >>>>>>>>>> loaded via a namespace (and not attached): >>>>>>>>>> [1] affy_1.34.0 affyio_1.24.0 BiocInstaller_1.4.3 >>>>>>>>>> [4] gdata_2.8.2 gplots_2.10.1 gtools_2.6.2 >>>>>>>>>> [7] preprocessCore_1.18.0 stats4_2.15.0 zlibbioc_1.2.0 >>>>>>>>>> >>>>>>>>>> >>>>>>>>> >>>>>>>> >>>>>>>> <110614 BENIGN_1 DATA 96x96.xlsx> >>>>>>> >>>>>> >>>>> >>>> >>> >> > NOTICE AND DISCLAIMER This e-mail (including any attachments) is intended for ...{{dropped:16}} ADD COMMENTlink written 7.2 years ago by Heidi Dvinge40 Answer: HTqPCR 0 7.2 years ago by Heidi Dvinge2.0k Heidi Dvinge2.0k wrote: Hello Deborah, > Good morning Heidi, > > thank you for your response. > > I have 6 samples with the "LightCycler" format (96 genes one each sample) > : one control and its duplicate, one "thirty minutes" and its duplicate > and one "three hours" and its duplicate. > I did the ttestCtdata, the mannwhitneyCtData and the limmaCtdata on my > samples for compare all the results. > I did the commands : >>PlaqueControle<-readCtData("essai_CT.txt", >> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96, >> format="LightCycler") >>Plaque30min<-readCtData("essai_30min.txt", >> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96, >> format="LightCycler") >>Plaque3h<-readCtData("essai_3h.txt", >> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96, >> format="LightCycler") >>PlaqueControleB<-readCtData("essai_CT_BIS.txt", >> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96, >> format="LightCycler") >>Plaque30minB<-readCtData("essai_30min_BIS.txt", >> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96, >> format="LightCycler") >>Plaque3hB<-readCtData("essai_3h_BIS.txt", >> column.info=list(feature="Name",position="Pos",Ct="Cp"), n.features=96, >> format="LightCycler") > >>essai<-cbind(PlaqueControle,Plaque30min,Plaque3h,PlaqueControleB,Pla que30minB,Plaque3hB) >>files_essai<-read.table("files_essai.txt",header=T) This all looks okay. Incidentally, you can also read in all the files together, using readCtData(files=c("essai_CT.txt", "essai_30min.txt", ...), n.data=6), in case you don't what to do the cbind() afterwards. >>essai.cat<- setCategory(essai, groups = files_essai$Treatment,quantile = >> 0.8) >>deltaCtnorm <- normalizeCtDataessai.cat, norm = "deltaCt",deltaCt.genes >> =? c("gene85", "gene86", "gene87", "gene88", "gene89", "gene90", >> "gene91")) Just out of curiosity, did you try other normalisation methods as well? I've seen a few cases where one of hte supposed housekeeper genes has been really off the chart. > > Then I applied the ttestCtData, I compared the "control" samples with the > "three hours" samples : >>essai_ttest <- ttestCtData(deltaCtnorm[,c(1,3,6,4)], groups = >> files_essai$Treatment[1:4],calibrator = "Controle") > ? > mannwhitneyCtData on the same samples : >>essaimwtest<-mannwhitneyCtData(deltaCtnorm[,c(1,3,6,4)], groups = >> files_essai$Treatment[1:4],calibrator = "Controle") > Looks okay. I assume the order of your samples is different in your qPCRset compared to your files_essai$Treatment. > > And finally, I did the limmaCtData : >>essai_design<-model.matrix(~0 + files_essai$Treatment) >>colnames(essai_design) <- c("Controle", "Trente_m", "Trois_h") >>print(essai_design) > Again, your sample order in files_essai is definitely different from your object, but you do seem to use c(1,3,6,4,2,5) further down, so I guess it should work. >>essai_contrasts <- makeContrasts(Trois_h-Controle,Trois_h - Trente_m, >> Trente_m - Controle, (Trente_m +Trois_h)/2 - Controle, levels = >> essai_design) >>colnames(essai_contrasts) <- c("3h-CT", "3h-30min", "30min-CT", >> "30min&3h-CT") >>print(essai_contrasts) > > >>deltaCtnorm2 <- deltaCtnorm[order(featureNames(deltaCtnorm)),] > You don't actually have to reorder when ndups=1. Of coruse, it doesn't hurt though. >>essai_limma <- limmaCtData(deltaCtnorm2[,c(1,3,6,4,2,5)], design = >> essai_design,contrasts = essai_contrasts, ndups = 1) > > > I didn't obtain the same conclusions for the same genes in the three > tests. So I think that I did a lot of error on the commands... > >From just glancing over it, it looks okay, although I of course can't actually see what your initial and resulting objects look like. > I give you the head screenshots of my results for the three tests in > attachments. > > For example, if I take the "gene10", with limmaCtData, I obtained that > "gene10" in the "control" samples is significatively different to "gene10" > in "three hours" samples (p-value = 0.0057). With ttestCtData, "gene10" is > not significatively different to "gene 10" in the "three hours" samples > (p-value = 0.063) ; ditto with mannwhitneyCtData (p-value =0.25). > > In the part 10.3 of "HTqPCR.pdf", you say about limma that "The result is > a list with one component per comparison. Each component is similar to the > result from using ttestCtData." ; so I suppose that my results are not > consistent. > The results from limma are *in principle* similar to a t-test, i.e. each component is the result from one of the individual tests specified by your contrasts. However since you're using 3 different types of statistical tests, it's not surprising that the results you get vary. Mann-Whitney is mainly used when the data aren't normally distributed. It's a non-parametric test, so it has less statistical power than the other two. It's therefore not surprising that the p-values are less significant. limma uses a modified, more advanced version of the standard t-test. It always considers data from all samples, not just the ones being compared for any given contrast, and it "borrows" data across all features on your assay to achieve a more robust estimate. It is thus expected to give more significant values that a standard t-test, which only considers the samples and genes being compared in each individual test in isolation. There's no hard'n'fast rule for which of the 3 tests to use, since it also depends on the (distribution of) your starting data etc. It also depends on what you're using your data for. Do you want something where you're 100% sure that the genes do indeed differ between your conditions, e.g. as a validation. Or is this preliminary data, that will be used for further studies in the lab, so you're more interested in e.g. the top 1-10 hits. Depending on this, you can always have a look at the actual data for some of your genes (either the values themselves or some of the plots from HTqPCR). The tests mainly differ in how conservative they are. With this in mind, you can check whether your results are relatively consistent by for example comparing the order of the resulting p-value. For example, if you plot the p-values from one test versus another, are they then completely randomly scattered, or do they follow a general trend. At the moment it sounds like your results mainly differ in the level of the p-value, not whether a gene is e.g. up- or down-regulated in a given comparison. If the latter is the case, then it sounds like there's definitely something wrong, and I'll have to look into that. I'm sorry that I can't give you any simple reply, or tell you which of the 3 tests to use. But it really depends on both the purpose of your study, and your data. \Heidi > Furthermore, all my results with the mannwhitneyCtData are not > significant... I don't know if it was the good way to use these tests. > > Sincerely, > > Deborah. > > > ________________________________ > De?: Heidi Dvinge <heidi at="" ebi.ac.uk=""> > ??: Deborah Ung <ung.deborah at="" yahoo.fr=""> > Envoy? le : Vendredi 22 juin 2012 11h33 > Objet?: Re: HTqPCR > > Hello Deborah, >> >> Good morning Heidi Dvinge, >> >> I am a french student and I am currently a trainee in a biotechnology >> company. I would like to know if you have other documentations about >> Limma >> applied to HTqPCR because I have some problems to analyze my results. >> > I'm afraid the only documentation is the examples in the help files for > limmaCtData, and the examples in the vignette section 10, i.e. > > ?limmaCtData > openVignette("HTqPCR") > > If you can email the commands you used and the resulting error messages, > (or why you think the results don't make sense), we can try to dissect the > problem. > > Best > \Heidi > >> Sincerely, >> >> Deborah. >>
0
7.2 years ago by
Heidi Dvinge2.0k
Heidi Dvinge2.0k wrote:
Hi Deborah, > Good morning Heidi, > > Yes, the order of my samples is different in my qPCRset compared to my > files_essai$Treatment. Do I have to order them in the same way ? > The order doesn't have to be the same, as long as you make sure you always remember re-order one of them when you use them both in the same function, such as the limmaCtData examples you provided last time. I don't trust myself to always remember such things (or to do it correctly!), so I always order the samples in the same way ;). But that's not a requirement. > <snip> > > And with normalizeCtData I obtained : >> deltaCtnorm <- normalizeCtDataessai.cat, norm = >> "deltaCt",deltaCt.genes=? c("gene85", "gene86", "gene87", "gene88", >> "gene89", "gene90","gene91")) > Calculating deltaCt values > ??????? Using control gene(s): gene85 gene86 gene87 gene88 gene89 gene90 > gene91 > ??????? Card 1: Mean=20.93????? Stdev=2.74 > ??????? Card 2: Mean=21.58????? Stdev=2.73 > ??????? Card 3: Mean=21.73????? Stdev=2.81 > ??????? Card 4: Mean=20.96????? Stdev=2.73 > ??????? Card 5: Mean=21.69????? Stdev=2.69 > ??????? Card 6: Mean=21.73????? Stdev=2.83 > I chose this method because my director gave me a file where he has chosen > only seven housekeeper genes on the eight (one of them has different Cp > results in almost each sample) so I did the same thing. > In this case the values of all these housekeepers do look robust across your samples, although the Ct values are higher than for typical housekeepers such as b-actin. (Which BTW isn't necessarily a bad thing). So as long as you/your boss is happy with it, I guess that's fine. > And about the other normalization methods, I tried them only for seeing > the difference between them. > I did : >>par(mfrow=c(3,2)) >>plot(exprs(essai),exprs(essai_g.mean),pch=20,main="Normalisation avec >> geometric.mean",col=rep(brewer.pal(6,"Spectral"),each=96)) >>plot(exprs(essai),exprs(essai_scale.rank),pch=20,main="Normalisation avec >> scale.rankinvariant", col=rep(brewer.pal(6,"Spectral"), each=96)) >>plot(exprs(essai),exprs(essai_deltaCt),pch=20,main="Normalisation avec >> deltaCt",col=rep(brewer.pal(6,"Spectral"),each=96)) >>plot(exprs(essai),exprs(essai_q.norm),pch=20,main="Normalisation avec >> quantile",col=rep(brewer.pal(6,"Spectral"),each=96)) >>plot(exprs(essai),exprs(essai_norm.rank),pch=20,main="Normalisation avec >> norm.rankinvariant", col=rep(brewer.pal(6,"Spectral"), each=96)) > In your case you have relatively few genes (96), which may not be quite enough for some of the methods. If the deltaCt-normalised data doesn't look too discrepant from all the other methods you're probably fine. > And then I wanted to compare only one of the sample as your example and I > used abline() but it didn't work. >>plot(exprs(DU145)[,3],exprs(essai_g.mean)[,3],pch=20,col="magenta") >>abline(exprs(DU145)[,3],exprs(essai_scale.rank)[,3],pch=20,col="blue ") >>abline(exprs(DU145)[,3],exprs(essai_deltaCt)[,3],pch=20,col="purple" ) > I did'nt get error message but there was only one plot. I also changed the > xlim value... > Well, you only use one plot() command, hence only one plot gets produced. Do you perhaps want to add data from multiple objects using points()? plot(exprs(DU145)[,3],exprs(essai_g.mean)[,3],pch=20,col="magenta") points(exprs(DU145)[,3],exprs(essai_norm.rank)[,3],pch=20,col="blue") ...etc.. > About the p.value, I tried to plot them as you had suggest me to do. > Here are the p.values of each test. >> MWTEST<-read.csv("essai_3h_CT_mwest1.csv", sep = ";", >> dec=",",header=TRUE) >> MWTEST$p.value > ?[1] 0.2452781 0.2452781 0.2452781 1.0000000 0.2452781 1.0000000 0.2452781 > 0.2452781 0.2452781 0.2452781 0.2452781 0.2452781 0.2452781 0.2452781 > 0.2452781 > [16] 0.2452781 1.0000000 0.2452781 0.2452781 0.6985354 0.2452781 0.2452781 > 0.2452781 0.6985354 0.2452781 0.6985354 0.2452781 0.6985354 0.2452781 > 0.2452781 > [31] 0.2452781 0.2452781 0.2452781 0.2452781 0.2452781 0.2452781 0.2452781 > 0.2452781 0.2452781 1.0000000 0.2452781 1.0000000 0.2452781 0.2452781 > 0.2452781 > [46] 0.2452781 0.2452781 0.2452781 1.0000000 0.2452781 0.6985354 0.2452781 > 0.6985354 0.2452781 0.2452781 0.2452781 0.2452781 0.2452781 1.0000000 > 0.2452781 > [61] 0.2452781 0.2452781 0.2452781 0.2452781 0.2452781 1.0000000 0.2452781 > 0.2452781 0.2452781 0.6985354 0.2452781 0.2452781 0.2452781 0.2452781 > 0.2452781 > [76] 0.2452781 0.2452781 0.2452781 0.2452781 1.0000000 0.2452781 0.2452781 > 0.2452781 0.2452781 > >> TTEST<-read.csv("essai_3h_CT_ttestFINAL.csv",sep=";",dec=",") >> TTEST$p.value > ?[1] 0.0000445916 0.0007203397 0.0011614062 0.0013499661 0.0021628924 > 0.0026335148 0.0028568178 0.0031900227 0.0047685877 0.0048604855 > 0.0058520334 > [12] 0.0061487450 0.0101863324 0.0101863324 0.0101863324 0.0101863324 > 0.0101863324 0.0101863324 0.0101863324 0.0101863324 0.0118721483 > 0.0130066108 > [23] 0.0166399589 0.0169334479 0.0170523949 0.0209596603 0.0364914159 > 0.0411074450 0.0427607922 0.0448349199 0.0476675408 0.0494563227 > 0.0512518277 > [34] 0.0514256601 0.0572253657 0.0625769920 0.0656449529 0.0797613537 > 0.0804359345 0.0819791697 0.0821601769 0.0909080090 0.0918788852 > 0.0986962901 > [45] 0.0993682993 0.1208709609 0.1261907225 0.1331849915 0.1338834931 > 0.1590798074 0.1611312123 0.1657960803 0.1718536210 0.1844114022 > 0.2035267116 > [56] 0.2092967748 0.2111576859 0.2192894000 0.2223339619 0.2393817321 > 0.2416781885 0.2479843103 0.2570206534 0.2570800840 0.2610404909 > 0.2755461365 > [67] 0.2886380998 0.3133822666 0.4574691996 0.4790123864 0.4963391483 > 0.5780428714 0.5827604076 0.6029711831 0.6738622120 0.6905548966 > 0.7800699292 > [78] 0.8045384637 0.8399336418 0.9347460142 0.9531859762 0.9719743053 > 0.9774759886 0.9934469629 > >> LIMMATEST<-read.csv("essai_limmaFINAL.csv",sep=";",dec=",") >> LIMMATEST$X3h.CT.p.value > ?[1] 4.538818e-01 9.722424e-01 9.681357e-01 1.478327e-02 9.765722e-01 > 1.224899e-01 4.579647e-01 1.035570e-03 6.190137e-02 6.862192e-03 > 1.828032e-02 > [12] 2.354413e-02 2.027634e-01 3.924912e-03 1.245438e-03 9.538478e-07 > 5.915158e-04 2.714424e-01 2.449646e-04 9.943747e-02 8.115928e-02 > 1.014429e-01 > [23] 1.730959e-04 2.283943e-01 4.106429e-02 8.292733e-01 7.384857e-01 > 9.053543e-04 3.031922e-05 4.381594e-02 8.697809e-05 1.730959e-04 > 9.949012e-01 > [34] 3.584419e-03 3.713434e-04 7.691588e-01 1.336464e-01 3.141131e-01 > 3.500428e-02 5.853026e-06 6.234777e-02 1.096195e-01 5.065608e-01 > 1.425943e-02 > [45] 7.720779e-01 2.074906e-05 2.596116e-04 6.080595e-02 6.472036e-01 > 1.730959e-04 6.924510e-05 8.243564e-03 2.010885e-01 9.367344e-01 > 2.535135e-01 > [56] 9.788777e-01 1.730959e-04 6.717992e-02 1.041109e-01 3.951307e-04 > 1.152792e-01 2.552804e-04 8.276034e-01 6.578508e-03 3.226937e-02 > 1.730959e-04 > [67] 1.730959e-04 1.052211e-01 6.826300e-05 1.730959e-04 2.939883e-01 > 1.116254e-02 2.997326e-01 5.701757e-02 2.319393e-03 3.023084e-02 > 8.304573e-01 > [78] 4.892519e-01 6.178556e-01 4.863336e-01 8.506124e-02 1.730959e-04 > 1.380221e-01 7.850957e-03 > > So I did : > par(mfrow=c(1,3)) > plot(LIMMATEST$X3h.CT.p.value,col="green",main="p-value LIMMA") > plot(TTEST$p.value,col="red",main="p-value Student") > plot(MWTEST$p.value,col="blue",main="p-value Mann-Whitney") > Considering the graphs that I obtained, I can say that the p.values don't > follow a general trend... So there is a real problem somewhere... Is that > alright ? > I'm sorry for not being clear, I meant plot them against each other. For example plot(TTEST$p.value, MWTEST\$p.value). > Actually, I have a question about the "Summary" component of the > limmaCtData : how do it do to calculate if it is up- or down- regulation ? > Because when I calculated the expression (2^ddCt) of the gene18, I > obtained expression = 11.41. > So the gene18 is up-regulated whereas in the "Summary" there is no change. > Is there a link between the "Summary" and the expression ? > In Summary, -1/0/1 should correspond to down-regulation/no difference/up-regulation respectively. The summary is linked to the expression, but it requires that the change in expression is statistically significant at p<0.05. Otherwise it's just "0" in the output. Best, \Heidi > Thank you again for your help and your advice, > > Deborah. >