Hi.
2009/6/24 marco <[email protected]>:
>
> Dear Henrik,
>
> I wonder if aroma.affymetrix is suitable for pooled GWA studies 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
-~----------~----~----~----~------~----~------~--~---
