Re: [aroma.affymetrix] ACNE for matched pairs
Hi Henrik,
Thanks for the advice. Now it works OK but I have problems in getting ACNE
results. I¹m using the X Titration samples from Affymetrix.
This is the code:
## I could use doASCRMAv2() instead for this part
library(aroma.affymetrix)
cdf - AffymetrixCdfFile$byChipType(GenomeWideSNP_6, tags=Full)
cs.comb- AffymetrixCelSet$byName(XTitration, cdf = cdf)
samplenames - getNames(cs.comb)
ref.normals - grep(2X, samplenames, fixed=TRUE)
acc - AllelicCrosstalkCalibration(cs.comb,model=CRMAv2)
csC - process(acc, verbose=verbose)
bpn - BasePositionNormalization(csC, target=zero)
csN - process(bpn, verbose=verbose)
## Now ACNE. I use as reference the arrays with 2 copies of X chromosome.
plm -NmfSnpPlm(csN, mergeStrands=TRUE, refs=ref.normals);
if (length(findUnitsTodo(plm)) 0) {
+ units - fitCnProbes(plm, verbose=verbose)
+ str(units)
+ units - fit(plm, verbose=verbose)
}
ces - getChipEffectSet(plm)
## And save the data
for (s in 1:length(samplenames)) {
+ cf - getFile(ces, s)
+ alldata - NULL
+ for (chr in 1:24) {
+chromosome - chr
+ units - getUnitsOnChromosome(gi, chromosome=chromosome)
+pos - getPositions(gi, units=units)
+chrom - rep(chr, length(pos))
+data - extractTotalAndFreqB(cf, units=units)
+data - cbind(chrom, pos, data)
+alldata - rbind(alldata, data)
+}
+write.table(alldata, file=paste(/ResultsACNEXTitration/,
samplenames[s], .txt, sep=\t))
}
However, I get results that look OK for BAF (although a little bit noisy)
but weird for the total copy number:
summary(alldata)
chrompostotal freqB
Min. : 1.00 Min. : 527 Min. : -136.139 Min. :0.000e+00
1st Qu.: 4.00 1st Qu.: 32618193 1st Qu.:2.107 1st Qu.:4.868e-16
Median : 8.00 Median : 69754406 Median : 563.033 Median :4.674e-01
Mean : 9.36 Mean : 78882883 Mean : 3354.815 Mean :4.748e-01
3rd Qu.:14.00 3rd Qu.:115310488 3rd Qu.: 5680.981 3rd Qu.:9.933e-01
Max. :24.00 Max. :247191012 Max. :85284.445 Max. :1.000e+00
NA's :9.458e+05
That is, total has ranges from -136 to 85,284.
Do you know what I'm doing wrong here?
Thanks a lot,
Oscar
On 5/11/10 00:12, Henrik Bengtsson henrik.bengts...@aroma-project.org
wrote:
Hi.
On Thu, Nov 4, 2010 at 3:11 AM, Oscar Rueda oscar.ru...@cancer.org.uk wrote:
Hi Henrik,
Many thanks for your detailed response even when I didn¹t include some
important information.
See comments below.
On 3/11/10 23:00, Henrik Bengtsson henrik.bengts...@gmail.com wrote:
Hi.
On Wed, Nov 3, 2010 at 6:03 AM, Oscar Rueda oscar.ru...@cancer.org.uk
wrote:
Hi all,
I'm trying to normalize several thousands of tumor samples using aroma.
The
BAF profiles I get using
AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=FALSE)
are quite noisy, so I have started using ACNE for getting 'cleaner' ones.
Yes, BAF profiles are quite noisy if not taking into account
SNP-specific crosstalk effects, which ACNE and some other methods do.
It is not clear from your message what chip type you are working on,
but here I will assume GenomeWideSNP_6.
Yes, This is GenomeWideSNP_6.
The problem is that I don't want to use the tumors as a reference, because
some of them have a lot of alterations.
I understand this concern. Though, note that with a robust estimator
together with the assumption that at any given SNP the majority of
samples/arrays are normal/diploid, - (CA,CB) = (2,0), (1,1) or (0,2) -
then the estimates of the crosstalk parameters should be ok. However,
I agree that if you know a sample is not-diploid it is best to exclude
it. The safe version of the latter is to only use normal samples.
Before answering your specific ACNE-related questions/problems, are
you aware of the TumorBoost method which is designed for the case
where you have a matched normal for each tumor?
Bengtsson, H.; Neuvial, P. Speed, T. P. TumorBoost: Normalization of
allele-specic tumor copy numbers from a single pair of tumor-normal
genotyping microarrays BMC Bioinformatics, 2010. (Downloadable via
http://aroma-project.org/publications/)
Using that, you can normalize the tumor BAFs utilizing the match
normal so that the normalized BAFs have much(!) greater
signal-to-noise ratios. It won't/cannot normalize your normal samples
- only the tumors - but if that is all you need it may be easier to
use. It also needs a matched normal for *each* tumor; you cannot
normalize tumors without a matched normal. Another advantage is that
if you preprocess with AS-CRMAv2 and BAF normalize with TumorBoost you
have a truly single-pair processing method, i.e. you can process your
tumor-normal pairs one by one independently of all other pairs you
have.
You will find vignettes from TumorBoost at
http://aroma-project.org/vignettes/. FYI, we are
Re: [aroma.affymetrix] ACNE for matched pairs
On Wed, Nov 10, 2010 at 7:20 AM, Oscar Rueda oscar.ru...@cancer.org.uk wrote:
Hi Henrik,
Thanks for the advice. Now it works OK but I have problems in getting ACNE
results. I¹m using the X Titration samples from Affymetrix.
This is the code:
## I could use doASCRMAv2() instead for this part
library(aroma.affymetrix)
cdf - AffymetrixCdfFile$byChipType(GenomeWideSNP_6, tags=Full)
cs.comb- AffymetrixCelSet$byName(XTitration, cdf = cdf)
samplenames - getNames(cs.comb)
ref.normals - grep(2X, samplenames, fixed=TRUE)
acc - AllelicCrosstalkCalibration(cs.comb,model=CRMAv2)
csC - process(acc, verbose=verbose)
bpn - BasePositionNormalization(csC, target=zero)
csN - process(bpn, verbose=verbose)
## Now ACNE. I use as reference the arrays with 2 copies of X chromosome.
plm -NmfSnpPlm(csN, mergeStrands=TRUE, refs=ref.normals);
if (length(findUnitsTodo(plm)) 0) {
+ units - fitCnProbes(plm, verbose=verbose)
+ str(units)
+ units - fit(plm, verbose=verbose)
}
ces - getChipEffectSet(plm)
## And save the data
for (s in 1:length(samplenames)) {
+ cf - getFile(ces, s)
+ alldata - NULL
+ for (chr in 1:24) {
+ chromosome - chr
+ units - getUnitsOnChromosome(gi, chromosome=chromosome)
+ pos - getPositions(gi, units=units)
+ chrom - rep(chr, length(pos))
+ data - extractTotalAndFreqB(cf, units=units)
+ data - cbind(chrom, pos, data)
+ alldata - rbind(alldata, data)
+ }
+ write.table(alldata, file=paste(/ResultsACNEXTitration/,
samplenames[s], .txt, sep=\t))
}
However, I get results that look OK for BAF (although a little bit noisy)
but weird for the total copy number:
summary(alldata)
chrom pos total freqB
Min. : 1.00 Min. : 527 Min. : -136.139 Min. :0.000e+00
1st Qu.: 4.00 1st Qu.: 32618193 1st Qu.: 2.107 1st Qu.:4.868e-16
Median : 8.00 Median : 69754406 Median : 563.033 Median :4.674e-01
Mean : 9.36 Mean : 78882883 Mean : 3354.815 Mean :4.748e-01
3rd Qu.:14.00 3rd Qu.:115310488 3rd Qu.: 5680.981 3rd Qu.:9.933e-01
Max. :24.00 Max. :247191012 Max. :85284.445 Max. :1.000e+00
NA's :9.458e+05
That is, total has ranges from -136 to 85,284.
That is because they are total signals (theta), not total copy-number
ratios (C). In order to get CNs you need to take the ratio toward a
reference (thetaR), e.g. the robust average of all arrays, e.g. C = 2
* theta / thetaR. This is, by definition, not an issue for BAF.
That you get some negative values is not strange either, because they
may occur after subtracting offset (background) in the signals.
Think about them as noise around zero.
/Henrik
Do you know what I'm doing wrong here?
Thanks a lot,
Oscar
On 5/11/10 00:12, Henrik Bengtsson henrik.bengts...@aroma-project.org
wrote:
Hi.
On Thu, Nov 4, 2010 at 3:11 AM, Oscar Rueda oscar.ru...@cancer.org.uk
wrote:
Hi Henrik,
Many thanks for your detailed response even when I didn¹t include some
important information.
See comments below.
On 3/11/10 23:00, Henrik Bengtsson henrik.bengts...@gmail.com wrote:
Hi.
On Wed, Nov 3, 2010 at 6:03 AM, Oscar Rueda oscar.ru...@cancer.org.uk
wrote:
Hi all,
I'm trying to normalize several thousands of tumor samples using aroma.
The
BAF profiles I get using
AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=FALSE)
are quite noisy, so I have started using ACNE for getting 'cleaner' ones.
Yes, BAF profiles are quite noisy if not taking into account
SNP-specific crosstalk effects, which ACNE and some other methods do.
It is not clear from your message what chip type you are working on,
but here I will assume GenomeWideSNP_6.
Yes, This is GenomeWideSNP_6.
The problem is that I don't want to use the tumors as a reference, because
some of them have a lot of alterations.
I understand this concern. Though, note that with a robust estimator
together with the assumption that at any given SNP the majority of
samples/arrays are normal/diploid, - (CA,CB) = (2,0), (1,1) or (0,2) -
then the estimates of the crosstalk parameters should be ok. However,
I agree that if you know a sample is not-diploid it is best to exclude
it. The safe version of the latter is to only use normal samples.
Before answering your specific ACNE-related questions/problems, are
you aware of the TumorBoost method which is designed for the case
where you have a matched normal for each tumor?
Bengtsson, H.; Neuvial, P. Speed, T. P. TumorBoost: Normalization of
allele-specic tumor copy numbers from a single pair of tumor-normal
genotyping microarrays BMC Bioinformatics, 2010. (Downloadable via
http://aroma-project.org/publications/)
Using that, you can normalize the tumor BAFs utilizing the match
normal so that the normalized BAFs have much(!) greater
signal-to-noise ratios. It
Re: [aroma.affymetrix] ACNE for matched pairs
Hi Henrik,
Many thanks for your detailed response even when I didn't include some
important information.
See comments below.
On 3/11/10 23:00, Henrik Bengtsson henrik.bengts...@gmail.com wrote:
Hi.
On Wed, Nov 3, 2010 at 6:03 AM, Oscar Rueda oscar.ru...@cancer.org.uk wrote:
Hi all,
I'm trying to normalize several thousands of tumor samples using aroma. The
BAF profiles I get using
AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=FALSE)
are quite noisy, so I have started using ACNE for getting 'cleaner' ones.
Yes, BAF profiles are quite noisy if not taking into account
SNP-specific crosstalk effects, which ACNE and some other methods do.
It is not clear from your message what chip type you are working on,
but here I will assume GenomeWideSNP_6.
Yes, This is GenomeWideSNP_6.
The problem is that I don't want to use the tumors as a reference, because
some of them have a lot of alterations.
I understand this concern. Though, note that with a robust estimator
together with the assumption that at any given SNP the majority of
samples/arrays are normal/diploid, - (CA,CB) = (2,0), (1,1) or (0,2) -
then the estimates of the crosstalk parameters should be ok. However,
I agree that if you know a sample is not-diploid it is best to exclude
it. The safe version of the latter is to only use normal samples.
Before answering your specific ACNE-related questions/problems, are
you aware of the TumorBoost method which is designed for the case
where you have a matched normal for each tumor?
Bengtsson, H.; Neuvial, P. Speed, T. P. TumorBoost: Normalization of
allele-specic tumor copy numbers from a single pair of tumor-normal
genotyping microarrays BMC Bioinformatics, 2010. (Downloadable via
http://aroma-project.org/publications/)
Using that, you can normalize the tumor BAFs utilizing the match
normal so that the normalized BAFs have much(!) greater
signal-to-noise ratios. It won't/cannot normalize your normal samples
- only the tumors - but if that is all you need it may be easier to
use. It also needs a matched normal for *each* tumor; you cannot
normalize tumors without a matched normal. Another advantage is that
if you preprocess with AS-CRMAv2 and BAF normalize with TumorBoost you
have a truly single-pair processing method, i.e. you can process your
tumor-normal pairs one by one independently of all other pairs you
have.
You will find vignettes from TumorBoost at
http://aroma-project.org/vignettes/. FYI, we are aware that the
high-level version of TumorBoost for the aroma framework is still a
bit tedious to apply. We simply haven't decided on the final design
where to put the result files etc. Regardless, the results are still
correct. Depending on what you are doing downstream, you might find
the low-level/all-in-memory version much easier to use.
Thanks for this! I remember reading that paper some time ago, so I'll
definitely give it
a try for the matched pairs.
With aroma I used the 'target=0'
option to get log intensities
That option I do not follow?!? Where do you use 'target=0'? Even if
you get signals that are on average near zero, I suspect that you are
misinterpreting the results/output, which most likely are *not*
log-intensities. Please show some code - actually show all you code
to avoid ambiguities.
This is the code I use:
library(aroma.affymetrix)
cdf - AffymetrixCdfFile$byChipType(GenomeWideSNP_6, tags=Full)
gi - getGenomeInformation(cdf)
si - getSnpInformation(cdf)
geneinform - getGenomeInformation(cdf)
snpinform - getSnpInformation(cdf)
cs.comb- AffymetrixCelSet$byName(Project1, chipType =
GenomeWideSNP_6,cdf = cdf)
getFullName(cs.comb)
getPath(cs.comb)
getNames(cs.comb)
samplenames = getNames(cs.comb)
setCdf(cs.comb, cdf)
acc - AllelicCrosstalkCalibration(cs.comb,model=CRMAv2)
csC - process(acc, ram=28,verbose=verbose)
bpn - BasePositionNormalization(csC, target=zero)
csN - process(bpn, ram=28,verbose=verbose)
plm - AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=FALSE)
if (length(findUnitsTodo(plm)) 0) {
units - fitCnProbes(plm, verbose=verbose)
str(units)
units - fit(plm, ram=28, verbose=verbose)
str(units)
}
ces - getChipEffectSet(plm)
fln - FragmentLengthNormalization(ces, target=zero)
cesN - process(fln, verbose=verbose)
cdf.mono - getCdf(cesN);
for (kk in 1:24) {
units - getUnitsOnChromosome(geneinform, chromosome=kk);
pos - getPositions(geneinform, units=units);
unitNames - getUnitNames(cdf.mono, units=units);
theta - extractMatrix(cesN, units=1, field = theta);
theta - colnames(theta)
Vals - extractTheta(cesN, units=units)
Vals - apply(Vals, c(1,2), cbind)
rownames(Vals) - unitNames;
colnames(Vals) - rep(theta, rep(2, length(theta)))
pos.log2I.val-cbind(pos,Vals)
write.table(pos.log2I.val,
paste(Project1_log2I.pos_Chrm,kk,.txt,sep=))
}
So my understanding is that these are intensities, not ratios. I can obtain
Re: [aroma.affymetrix] ACNE for matched pairs
Hi. On Thu, Nov 4, 2010 at 3:11 AM, Oscar Rueda oscar.ru...@cancer.org.uk wrote: Hi Henrik, Many thanks for your detailed response even when I didn’t include some important information. See comments below. On 3/11/10 23:00, Henrik Bengtsson henrik.bengts...@gmail.com wrote: Hi. On Wed, Nov 3, 2010 at 6:03 AM, Oscar Rueda oscar.ru...@cancer.org.uk wrote: Hi all, I'm trying to normalize several thousands of tumor samples using aroma. The BAF profiles I get using AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=FALSE) are quite noisy, so I have started using ACNE for getting 'cleaner' ones. Yes, BAF profiles are quite noisy if not taking into account SNP-specific crosstalk effects, which ACNE and some other methods do. It is not clear from your message what chip type you are working on, but here I will assume GenomeWideSNP_6. Yes, This is GenomeWideSNP_6. The problem is that I don't want to use the tumors as a reference, because some of them have a lot of alterations. I understand this concern. Though, note that with a robust estimator together with the assumption that at any given SNP the majority of samples/arrays are normal/diploid, - (CA,CB) = (2,0), (1,1) or (0,2) - then the estimates of the crosstalk parameters should be ok. However, I agree that if you know a sample is not-diploid it is best to exclude it. The safe version of the latter is to only use normal samples. Before answering your specific ACNE-related questions/problems, are you aware of the TumorBoost method which is designed for the case where you have a matched normal for each tumor? Bengtsson, H.; Neuvial, P. Speed, T. P. TumorBoost: Normalization of allele-specic tumor copy numbers from a single pair of tumor-normal genotyping microarrays BMC Bioinformatics, 2010. (Downloadable via http://aroma-project.org/publications/) Using that, you can normalize the tumor BAFs utilizing the match normal so that the normalized BAFs have much(!) greater signal-to-noise ratios. It won't/cannot normalize your normal samples - only the tumors - but if that is all you need it may be easier to use. It also needs a matched normal for *each* tumor; you cannot normalize tumors without a matched normal. Another advantage is that if you preprocess with AS-CRMAv2 and BAF normalize with TumorBoost you have a truly single-pair processing method, i.e. you can process your tumor-normal pairs one by one independently of all other pairs you have. You will find vignettes from TumorBoost at http://aroma-project.org/vignettes/. FYI, we are aware that the high-level version of TumorBoost for the aroma framework is still a bit tedious to apply. We simply haven't decided on the final design where to put the result files etc. Regardless, the results are still correct. Depending on what you are doing downstream, you might find the low-level/all-in-memory version much easier to use. Thanks for this! I remember reading that paper some time ago, so I’ll definitely give it a try for the matched pairs. With aroma I used the 'target=0' option to get log intensities That option I do not follow?!? Where do you use 'target=0'? Even if you get signals that are on average near zero, I suspect that you are misinterpreting the results/output, which most likely are *not* log-intensities. Please show some code - actually show all you code to avoid ambiguities. This is the code I use: library(aroma.affymetrix) cdf - AffymetrixCdfFile$byChipType(GenomeWideSNP_6, tags=Full) gi - getGenomeInformation(cdf) si - getSnpInformation(cdf) geneinform - getGenomeInformation(cdf) snpinform - getSnpInformation(cdf) FYI, the above 4 lines are just there for assertions to make sure the annotation files are there before starting. cs.comb- AffymetrixCelSet$byName(Project1, chipType = GenomeWideSNP_6,cdf = cdf) You only need to specify one of 'chipType' and 'cdf' - don't specify both. If you do: cs.comb - AffymetrixCelSet$byName(Project1, cdf=cdf) you don't have to do the following/below setCdf(cs.comb, cdf) statement. getFullName(cs.comb) getPath(cs.comb) getNames(cs.comb) samplenames = getNames(cs.comb) setCdf(cs.comb, cdf) acc - AllelicCrosstalkCalibration(cs.comb,model=CRMAv2) csC - process(acc, ram=28,verbose=verbose) Instead of explicitly specifying the 'ram' argument here, which will cause issues if you bring the script somewhere else, you can set an aroma setting as: setOption(aromaSettings, memory/ram, 28); at the beginning of the script. Actually, you only have to do it once on the computer if you also save it afterward with saveAnywhere(aromaSettings). That way you don't have to put in the script at all, making you script totally transparent to system properties. See http://aroma-project.org/settings/ for more details. bpn - BasePositionNormalization(csC, target=zero) csN - process(bpn, ram=28,verbose=verbose) Same here. plm - AvgCnPlm(csN, mergeStrands=TRUE,
Re: [aroma.affymetrix] ACNE for matched pairs
Hi. On Wed, Nov 3, 2010 at 6:03 AM, Oscar Rueda oscar.ru...@cancer.org.uk wrote: Hi all, I'm trying to normalize several thousands of tumor samples using aroma. The BAF profiles I get using AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=FALSE) are quite noisy, so I have started using ACNE for getting 'cleaner' ones. Yes, BAF profiles are quite noisy if not taking into account SNP-specific crosstalk effects, which ACNE and some other methods do. It is not clear from your message what chip type you are working on, but here I will assume GenomeWideSNP_6. The problem is that I don't want to use the tumors as a reference, because some of them have a lot of alterations. I understand this concern. Though, note that with a robust estimator together with the assumption that at any given SNP the majority of samples/arrays are normal/diploid, - (CA,CB) = (2,0), (1,1) or (0,2) - then the estimates of the crosstalk parameters should be ok. However, I agree that if you know a sample is not-diploid it is best to exclude it. The safe version of the latter is to only use normal samples. Before answering your specific ACNE-related questions/problems, are you aware of the TumorBoost method which is designed for the case where you have a matched normal for each tumor? Bengtsson, H.; Neuvial, P. Speed, T. P. TumorBoost: Normalization of allele-specic tumor copy numbers from a single pair of tumor-normal genotyping microarrays BMC Bioinformatics, 2010. (Downloadable via http://aroma-project.org/publications/) Using that, you can normalize the tumor BAFs utilizing the match normal so that the normalized BAFs have much(!) greater signal-to-noise ratios. It won't/cannot normalize your normal samples - only the tumors - but if that is all you need it may be easier to use. It also needs a matched normal for *each* tumor; you cannot normalize tumors without a matched normal. Another advantage is that if you preprocess with AS-CRMAv2 and BAF normalize with TumorBoost you have a truly single-pair processing method, i.e. you can process your tumor-normal pairs one by one independently of all other pairs you have. You will find vignettes from TumorBoost at http://aroma-project.org/vignettes/. FYI, we are aware that the high-level version of TumorBoost for the aroma framework is still a bit tedious to apply. We simply haven't decided on the final design where to put the result files etc. Regardless, the results are still correct. Depending on what you are doing downstream, you might find the low-level/all-in-memory version much easier to use. With aroma I used the 'target=0' option to get log intensities That option I do not follow?!? Where do you use 'target=0'? Even if you get signals that are on average near zero, I suspect that you are misinterpreting the results/output, which most likely are *not* log-intensities. Please show some code - actually show all you code to avoid ambiguities. and then I normalized them against my pool of normals or against them matched normal if available just subtracting the values. I need to see your code to know what you are doing. For ACNE it seems I have to specify the reference, so I created folders for each of my matched pairs and used NmfSnpPlm(csN, mergeStrands=TRUE, refs=c(FALSE, TRUE)) To indicate that the second array is the normal and the first the tumor. Is this the best way to proceed in my situation? Nope. When you find yourself having to do tedious tricks such as splitting up your data set to data sets containing a single tumor-normal pair take it as sign for there must be a better way to do this (which I guess is one reason why you posted this message). The former command produces the following values for chromosome 2: summary(data) total freqB Min. : 0.0 Min. : NA 1st Qu.: 0.0 1st Qu.: NA Median : 905.6 Median : NA Mean : 2345.2 Mean : NaN 3rd Qu.: 4192.8 3rd Qu.: NA Max. :53200.6 Max. : NA NA's :153732 So it's clear that I'm doing something wrong. Does anyone know what is it? The problem is that you are asking the ACNE method to estimate the (SNP-specific) crosstalk parameters using only two arrays, actually only one since you specify that it is only the normal sample that should be used for the estimation. The model parameters are non-identifiable with only one sample, i.e. it is not possible to estimate the ACNE parameter. It should probably be added to the ACNE method/estimator (somewhere in the NmfSnpPlm class) to assert that not too few reference samples are used, instead of returning missing values. Instead, you should keep all your samples in the same directory and process the whole data set at once. You can specify which samples should be used as the reference set by using argument 'refs' to NmfSnpPlm, just as above. For example: refs - seq(from=2, to=length(csN), by=2); nmf - NmfSnpPlm(csN, mergeStrands=TRUE,
