Hi. On Wed, Jul 22, 2009 at 5:56 PM, marco<[email protected]> wrote: > > Dear Henrik, > > thans for the useful answers! > > RAS = Relative Allele Scorel is normally used as a frequency estimate > for pooled DNA experiments. I saw around that SNP6.0 seems not to be > used for this scope before, so there are not really validated methods > to estimate frequencies from it. > > 1. The RAS score is defined as RAS = I_B/(I_A+I_B) , so it seems to me > to be equal to "fracB" as long as I am correct to assume that for the > snps, using combineAlleles=FALSE, thetaA = median(PMA) and > thetaB=median(PMB). This is my understanding from the CRMA2 paper. > Using the notation of the paper: > thataA_i_j = median_k(Y_i_j_k_A)
Almost: RAS is the relative allele *A* signal such that RAS ~= 0 for BB [1] whereas the Allele B Fraction (fracB) is the relative allele *B* signal such that fracB ~= 1 for BB. Thus, RAS = 1-beta. Also, RAS is defined for each probe/cell pair (PMA,PMB) [using MM:s too] whereas the Allele B Fraction is defined for the summarized signals [regardless what summarization method was used]. However, in practice they should give very similar outcomes. [1] Affymetrix, GeneChip DNA Analysis Software (GDAS) User's Guide (Version 3.0), 2004. Yes, plm <- AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=FALSE) will ignore strand information (mergeStrands=TRUE) and summarize each allele separately (combineAlleles=FALSE; name from dChip). Since Avg[Cn]Plm is used, this means that the probe signals are averaged (using median by default). You get equivalent results if you do: plm <- AvgSnpPlm(csN, mergeStrands=TRUE) The AvgSnpPlm class does not "understand" total CNs, which is why argument 'combineAlleles' does not make sense. You can still do: ces <- getChipEffectSet(plm); ce <- getFile(ces, 1); data <- extractTotalAndFreqB(ce, units=units); > > 2. By allelc crosstalk normalization I obtain quite a number of > negative probe intensities (around 3%). Do you have nay suggestion > about how to handle those? I guess this is a general problem since > from my understanding the matrix S is calculated without constraints. Yes, we do not constraint the crosstalk model - we simple believe it is alright to have slightly negative estimates assuming negative noise exist (conceptually). Note that it is rather unlikely that corrected PM = PM_A + PM_B is negative and even less so that theta = theta_A+theta_B is negative. However, the allele B fractions fracB=thetaB/theta will not be constraint to [0,1] but will be slightly outside. If this is a problem, you can shift the (theta_A,theta_B) or (theta,fracB) signals manually, or by: plm <- AvgSnpPlm(csN, mergeStrands=TRUE, shift=+300) I consider this an ad hoc solution and I prefer to accept the fact that you can have negative signals and avoid using the logarithmic scale (which is often the real reason why people don't like non-positive signals). > > 3. Is there a way to upload figures on this group? I have done some > diagnositc plots that might be interesting for other users. You can attach PNGs to email messages posted to the mailing list, but you cannot do this via the web interface. /Henrik > > > Best Regards > > Marco > > > On 26 Juni, 07:35, Henrik Bengtsson <[email protected]> wrote: >> Hi. >> >> 2009/6/24 marco <[email protected]>: >> >> >> >> > Dear Henrik, >> >> > I wonder if aroma.affymetrix is suitable for pooledGWAstudies and >> > if you have any experience with this subject. In principle this >> > consists in extracting the RAS scores, so should be doable. I have the >> > following specific questions and snippet: >> >> > 1. by using : >> >> > plm <- AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=FALSE) >> > data <- extractTotalAndFreqB(ce, units=units); >> >> > - the freqB column should be the equivalent of RAS? >> >> Yes, there should be a one-to-one relationship between fracB and RAS. >> I don't have the definition of RAS in front of me, but fracB is >> defined as fracB = thetaB / (thetaA+thetaB) [non logged signals!] so >> you should be able to find the transform. >> >> (Nowadays, we prefer to call it allele-B fraction instead of allele-B >> frequency; the package will later reflect this too). >> >> > - how is it possible to get the SNPs names in the data.frame from >> > extractTotalAndFreqB >> >> Note, that function does *not* return a 'data.frame' but a 'matrix' or >> an 'array' (see below). >> >> unitNames <- getUnitNames(cdf, units=units); >> >> > - can the freqB for all the arrays be extracted at once, or only array >> > by array? >> >> fracB <- extractTotalAndFreqB(ces, units=units); >> >> returns a Jx2xI array where J=length(units) and I=nbrOfArrays(ces), whereas, >> >> ce <- getFile(ces, 1); >> fracB <- extractTotalAndFreqB(ce, units=units); >> >> returns a Jx2 matrix. >> >> >> >> > 2. Is there a way to eliminate CN probes from the analysis from the >> > very beginning (just SNPs are need for pooling)? >> >> Don't know what you mean from the very beginning, but you can always >> find the types of a CDF unit using: >> >> unitTypes <- getUnitTypes(cdf); >> >> SNPs have value 2 (this is according to the Affymetrix CDF file format), e.g. >> >> snpUnits <- which(unitTypes == 2); >> >> will identify the SNP units for you. >> >> >> >> > 3. Will be the standard normalization steps of CRMA usefuls/ >> > deleterious/neutral? >> >> Can you restate this question; not sure I understand it. >> >> >> >> > I eliminated the crosstalk step so far, since I am not sure how this >> > behaves on a DNA pool. It seems to me that the crosstalk requires 3 >> > clusters for AA, AB BB to work properly. Is this correct? >> >> Please see CRMA and CRMA v2 papers [1,2] which explains this. The >> crosstalk estimate should be fairly robust against what is going on >> inside the genotype "cone". It tries to fit the crosstalk model by >> giving much more weight to data points there the (A,AA,AAA,...) and >> (B,BB,BBB,...) radials. >> >> >> >> > 4. Same question for the pcr and probe sequence normalization steps. >> > Seems to me that they should not disturb >> >> They should be applicable as is. >> >> >> >> >> >> > 5. Do you see any missing/unuseful/wrong step in the code below? >> >> > ############# >> > cdf <- AffymetrixCdfFile$byChipType("GenomeWideSNP_6", >> > tags="Full") >> > gi <- getGenomeInformation(cdf) >> > si <- getSnpInformation(cdf) >> > acs <- AromaCellSequenceFile$byChipType(getChipType(cdf, >> > fullname=FALSE)) >> > cdf <- AffymetrixCdfFile$byChipType("GenomeWideSNP_6", >> > tags="Full") >> > csR <- AffymetrixCelSet$byName("POOLS_SNP", cdf=cdf) >> >> > #acc <- AllelicCrosstalkCalibration(csR, model="CRMAv2") >> > #print(acc) >> > #csC <- process(acc, verbose=verbose) >> > #print(csC) >> >> > bpn <- BasePositionNormalization(csR, target="zero") >> > csN <- process(bpn, verbose=verbose) >> > plm <- AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=FALSE) >> >> > if (length(findUnitsTodo(plm)) > 0) { >> > units <- fitCnProbes(plm, verbose=verbose) >> > str(units) >> > units <- fit(plm, verbose=verbose) >> > str(units) >> > } >> > ces <- getChipEffectSet(plm) >> > fln <- FragmentLengthNormalization(ces, target="zero") >> > cesN <- process(fln, verbose=verbose) >> >> > ce <- getFile(cesN, 1) >> > data <- extractTotalAndFreqB(ce, units=units); >> >> Look alright to me. >> >> REFERENCES: >> [1] H. Bengtsson; R. Irizarry; B. Carvalho; T.P. Speed, Estimation and >> assessment of raw copy numbers at the single locus level, >> Bioinformatics, 2008. [pmid: 18204055] >> [2] H. Bengtsson; P. Wirapati & T.P, Speed, A single-array >> preprocessing method for estimating full-resolution raw copy numbers >> from all Affymetrix genotyping arrays including GenomeWideSNP 5 & 6, >> Biinformatics, 2009. >> >> /Henrik >> >> PS. Please try to split up unrelated questions in separate messages. >> >> > > > --~--~---------~--~----~------------~-------~--~----~ 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 [email protected] 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 -~----------~----~----~----~------~----~------~--~---
