Hi Sabrina.
On 08/12/2008, at 3:25 PM, sabrina wrote: > > Hi, Mark: > Thanks for the help! I think now I have better understanding of LIMMA > now after I read the related chapter of LIMMA user's guide :) > > > On Dec 4, 4:10 pm, "Mark Robinson" <[EMAIL PROTECTED]> wrote: >> Sabrina. >> >> When asking questions like this, please post your full set of >> commands. >> Its hard to decode what variables have been set without this. :) >> >> I'm guessing that affyAnnotation is what came from the 'read.csv' >> call and >> x2 is the output of topTable? Have you set the rownames of exFirma? > > x2 is the output of topTable > I do not know how to set rownames of exFirma, when I check exFirma > [1:3,1:6], it gave me the following: > > unitName groupName unit group cell Array1 > 1 6838637 4304927 1 1 1 1.0372756 > 2 6838637 4330595 1 2 2 -1.1679380 > 3 6838637 4356771 1 3 3 1.2217411 > > > Can you show me how to set the rownames? where can I find the row > names? Is it correct that groupName is the exon id ? A good place to start: ?rownames > I used the following code to import affymetrix annotation file: > > affyAnnotation<- read.csv(file="MoEx-1_0-st- > v1.na27.mm9.transcript.csv",head=TRUE,sep=",",skip=20); > > and >> >> Is this just a typo then? It should be >> affyAnnotation$transcript_cluster_id instead of >> affyAnnotation$transcript_cluster_Id. >> >> What does exFirma[x2$ID,1] actually give you? If the >> rownames(exFirma) >> are not set, this could mislead you. > > I used the suggestion you gave to me to redo the lmFit, and here is > what I used : > > mm <- model.matrix(~cls) > fit <- lmFit(exFirma[,6:ncol(exFirma)], mm) > fit<-eBayes(fit) > x<-topTable(fit,coef=2,adj="fdr") > > > Then I used exFirma[x$ID,1], it gave me: > [1] "6996667" "6965348" "6982181" "6917790" "6899766" "7005797" > "6880577" "6962054" "6897349" > [10] "7018774" > > when I printed x, it gave me: > > ID logFC t P.Value adj.P.Val B > 11253 11253 -4.856783 -8.524568 2.223806e-06 0.4982258 -4.323317 > 183160 183160 4.058556 7.676406 6.413325e-06 0.7184271 -4.332955 > 32235 32235 -3.692429 -6.555475 2.963704e-05 1.0000000 -4.350176 > 132610 132610 2.315536 6.508698 3.170265e-05 1.0000000 -4.351036 > 165279 165279 3.269404 6.311540 4.225007e-05 1.0000000 -4.354809 > 171777 171777 3.845620 6.174738 5.172731e-05 1.0000000 -4.357574 > 45380 45380 1.971388 6.088251 5.886544e-05 1.0000000 -4.359388 > 132873 132873 -4.569042 -6.021868 6.505114e-05 1.0000000 -4.360815 > 86037 86037 3.206732 5.890372 7.943255e-05 1.0000000 -4.363738 > 27761 27761 -2.651559 -5.825434 8.774484e-05 1.0000000 -4.365229 > > > > when I typed exFirma[11253,], it gave me: > > exFirma[11253,] > unitName groupName unit group cell array1, ... > 11253 6996667 5092691 535 6 11253 2.718783 ... > > so it appeared to me that x$ID matches the index of exFirma entry. > > Then I used; > index<-match(exFirma[x$ID,1], affyAnnotation$transcript_cluster_id) > to find the index of these entries (significant different between two > groups) in Affy annotation, then I checked the entries in Affy > Annotation file: > affyAnnotation[index, 1:5] > > transcript_cluster_id probeset_id seqname strand start > 71605 6996667 6996667 chr9 - 70631116 > 62004 6965348 6965348 chr7 + 151009070 > 67351 6982181 6982181 chr8 - 47660950 > ... This all looks good to me. > But then I don't know how to check whether it is correct or not , or > how to find and plot corresponding probeset (exon ) as what was shown > in the FIRMA paper supplementary . Can you give me some hints? Thanks! Elizabeth made the plots in the paper and has some code for it. For starters, familiarize yourself with the 'GenomeGraphs' package -- note that in the GenomeGraphs vignette, it has an exon array example. I have a vignette for this, but haven't finished it yet. HTH, Mark > Sabrina > >> >> Cheers, >> Mark >> >> >> >>> Ok, I just used the following line to read Affy Annotation >> >>> read.csv(file="MoEx-1_0-st- >>> v1.na27.mm9.transcript.csv",head=TRUE,sep=",",skip=20) >> >>> then I used >>> match(exFirma[x2$ID,1],affyAnnotation$transcript_cluster_Id) >> >>> but I only got NAs >> >>> I can not figure out what went wrong because I checked the first ID, >>> it is in Affy annotation file. Can you point out what I did >>> incorrectly? Thanks ! >> >>> Sabrina >> >>> On Dec 3, 4:53 pm, Mark Robinson <[EMAIL PROTECTED]> wrote: >>>> Sabrina. >> >>>> See below. >> >>>> On 04/12/2008, at 8:41 AM, sabrina wrote: >> >>>>> Hi, Mark: >> >>>>> Thank you so much for the detailed explanation! >>>>> Actually I am working on Mouse Exon arrays. I did check the UCSC >>>>> genome browser, but they do not have the Exon Array available as >>>>> Human Exon. >> >>>> Right! Good point. >> >>>>> My study is trying to find alternative splicing events between two >>>>> groups with identical genotypes but different treatments. each >>>>> group >>>>> has 4 or 5 biological replicates. So I did all as described in the >>>>> examples from Human Exon array, then used the following codes to >>>>> find >>>>> the significant genes with potential splicing events: >> >>>>> ... >>>>> exFirma<-extractDataFrame(firmaScore,addNames=TRUE); >>>>> exFirma[,6:ncol(exFirma)]<-log2(exFirma[,6:ncol(exFirma)]) >> >>>>> library(limma) >> >>>>> design <- >>>>> cbind(Grp1=c(1,1,1,1,1,0,0,0,0),Grp2=c(0,0,0,0,0,1,1,1,1)) >> >>>>> fit<-lmFit(exFirma[,6:ncol(exFirma)],design) >>>>> fit<-eBayes(fit) >>>>> x<-topTable(fit) >> >>>> Careful here as to what you are fitting and testing with limma. >>>> You >>>> are fitting a parameterization that would require a contrast to >>>> be fit >>>> ( ... such as 'contrasts.fit' ... ) in order to assess the >>>> differences >>>> between your 2 groups. Did you get this example from a web page >>>> or a >>>> thread somewhere? If so, I should correct it. >> >>>> Probably what you want is: >> >>>> design <- cbind(Grp1=1,Grp2=c(0,0,0,0,0,1,1,1,1)) >> >>>> ... >> >>>> topTable(fit,coef=2) >> >>>> Under this parameterization, setting coef=2 will test the >>>> difference >>>> between the 2 groups. Note that topTable by default gives the >>>> top 10 >>>> values, regardless of statistics. >> >>>> Does that make sense? >> >>>>> However I only got about 10 IDs with all of which adjusted P >>>>> values as >>>>> 1. I wonder if there is anything I did wrong here. Also when I >>>>> checked >>>>> the IDs back to Affymetrix annotation file, I could only find one >>>>> match, I wonder what went wrong. >> >>>> How exactly did you match them back? >> >>>> Mark >> >>>>> Thanks! >> >>>>> BTW, I used the coreR1 , A20080718, MR, cdf version >> >>>>> Sabrina >> >>>>> On Dec 2, 6:27 pm, Mark Robinson <[EMAIL PROTECTED]> wrote: >>>>>> Hi Sabrina. >> >>>>>> Indeed, you'll want to set 'mergeGroups=TRUE' for the probe level >>>>>> model (ExonRmaPlm object) that you send to 'FirmaModel' ... >> >>>>>> ---------http://groups.google.com/group/aroma-affymetrix/web/human-exon-array- >>>>>> ... >>>>>> says: >>>>>> ... >>>>>> plmTr <- ExonRmaPlm(csN, mergeGroups=TRUE) >>>>>> print(plmTr) >>>>>> ... >>>>>> firma <- FirmaModel(plmTr) >>>>>> fit(firma, verbose=verbose) >>>>>> --------- >> >>>>>> NO, the unitName is not the same as gene id in NCBI. If you are >>>>>> using >>>>>> the CDF from: >> >>>>>> http://groups.google.com/group/aroma-affymetrix/web/affymetrix- >>>>>> define... >> >>>>>> (e.g. HuEx-1_0-st-v2,coreR3,A20071112,EP.cdf) the unitName is an >>>>>> Affymetrix-defined identifier. >> >>>>>> If you go to: >> >>>>>> http://www.affymetrix.com/products_services/arrays/specific/ >>>>>> exon.affx... >> >>>>>> you will see a 'HuEx-1_0-st-v2 Probeset Annotations' file that >>>>>> you >>>>>> can >>>>>> download. In there, you will find the annotations you need >>>>>> (unitName >>>>>> corresponds to the 'transcript_cluster_id' column and groupName >>>>>> corresponds to the 'probeset_id' column) >> >>>>>> Also, you may be interested to know that the UCSC genome browser >>>>>> has a >>>>>> track called 'Affy HuEx 1.0' under the set of 'Expression' >>>>>> tracks. >>>>>> Unfortunately, its not searchable on the browser, but quite >>>>>> useful if >>>>>> you can locate it in the browser and find out which exon is >>>>>> differentially spliced. >> >>>>>> Hope that helps. >>>>>> Mark >> >>>>>> On 03/12/2008, at 5:56 AM, sabrina wrote: >> >>>>>>> Hi, all: >>>>>>> I just started learning to use firma analysis. My impression >>>>>>> from >>>>>>> previous discussions of other users is that we should only use >>>>>>> plm >>>>>>> with mergeGroup=True, is that correct? Another question is : Is >>>>>>> unitName same as gene uid in NCBI database? how can I map >>>>>>> groupName to >>>>>>> EXON name from NCBI? Thanks !!! >> >>>>>>> Sabrina >> >>>>>> ------------------------------ >>>>>> Mark Robinson >>>>>> Epigenetics Laboratory, Garvan >>>>>> Bioinformatics Division, WEHI >>>>>> e: [EMAIL PROTECTED] >>>>>> e: [EMAIL PROTECTED] >>>>>> p: +61 (0)3 9345 2628 >>>>>> f: +61 (0)3 9347 0852 >>>>>> ------------------------------ >> >>>> ------------------------------ >>>> Mark Robinson >>>> Epigenetics Laboratory, Garvan >>>> Bioinformatics Division, WEHI >>>> e: [EMAIL PROTECTED] >>>> e: [EMAIL PROTECTED] >>>> p: +61 (0)3 9345 2628 >>>> f: +61 (0)3 9347 0852 >>>> ------------------------------ > > ------------------------------ Mark Robinson Epigenetics Laboratory, Garvan Bioinformatics Division, WEHI e: [EMAIL PROTECTED] e: [EMAIL PROTECTED] p: +61 (0)3 9345 2628 f: +61 (0)3 9347 0852 ------------------------------ --~--~---------~--~----~------------~-------~--~----~ When reporting problems on aroma.affymetrix, make sure 1) to run the latest version of the package, 2) to report the output of sessionInfo() and traceback(), and 3) to post a complete code example. You received this message because you are subscribed to the Google Groups "aroma.affymetrix" group. To post to this group, send email to aroma-affymetrix@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/aroma-affymetrix?hl=en -~----------~----~----~----~------~----~------~--~---
