Mouse diversity array --building the required files for aroma.affymetrix UGP, UFL
Ivanek, Robert
I would like to analyse the copy number variation data from Affymetrix
Mouse Diversity Array. I have not found any information on your website
about this particular array. I have tried to build the annotation files
which are required by aroma but without success. I have few questions
regarding that:
1: Is aroma.affymetrix able to analyse the "Mouse Diversity Array" ?
2: I tried to build the "UGP" file directly from NetAffx annotation
files using the code on your website, however I am getting the following
error.
##
library("aroma.affymetrix")
##
## create UGP from NetAffx files
cdf <- AffymetrixCdfFile$byChipType("MOUSEDIVm520650")
##
## Creates an empty UGP file for the CDF, if missing.
ugp <- AromaUgpFile$allocateFromCdf(cdf, tags=c("na30", "RI20100526"))
##
## Import NetAffx unit position data
csv <- AffymetrixNetAffxCsvFile$byChipType("MOUSEDIVm520650",
otags=".na30")
Error in list(`AffymetrixNetAffxCsvFile$byChipType("MOUSEDIVm520650",
tags = ".na30")` = <environment>, :
[2010-05-26 15:11:00] Exception: File format error of the tabular file
('annotationData/chipTypes/MOUSEDIVm520650/NetAffx/MOUSEDIVm520650.na30.annot.csv'): \
line 1 did not have 12 elements
at throw(Exception(...))
at throw.default("File format error of the tabular file ('",
getPathname(this), "'): ", ex$message)
at throw("File format error of the tabular file ('",
getPathname(this), "'): ", ex$message)
at value[[3]](cond)
at tryCatchOne(expr, names, parentenv, handlers[[1]])
at tryCatchList(expr, classes, parentenv, handlers)
at tryCatch({
at verify.TabularTextFile(this, ...)
at verify(this, ...)
at this(...)
at newInstance.Class(clazz, ...)
at newInstance(clazz, ...)
at newInstance.Object(static, pathname)
at newInstance(static, pathname)
at method(static, ...)
at AffymetrixNetAffxCsvFile$byChipType("MOUSEDIVm520650", tags =
".na30")
In addition: Warning message:
In read.table(3L, header = TRUE, colClasses = c(NA_character_,
NA_character_, :
not all columns named in 'colClasses' exist
3. I tried it also by using the "manual" approach using the
tab=delimited file, however it seems to me that the mitochondria probes
are skipped (NA values in ugp[,1] but valid values in ugp[,2]). Another
problem is that some positions for other chromosomes are not loaded in
properly (valid values in ugp[,1] but NA values in ugp[,2]).
Here is the sessionInfo:
R version 2.11.0 (2010-04-22)
x86_64-unknown-linux-gnu
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=C
[6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C
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] aroma.affymetrix_1.6.0 aroma.apd_0.1.7 affxparser_1.20.0
R.huge_0.2.0 aroma.core_1.6.0 aroma.light_1.16.0
[7] matrixStats_0.2.1 R.rsp_0.3.6 R.cache_0.3.0
R.filesets_0.8.1 digest_0.4.2 R.utils_1.4.0
[13] R.oo_1.7.2 R.methodsS3_1.2.0 RColorBrewer_1.0-2
loaded via a namespace (and not attached):
[1] tools_2.11.0
Best Regards
Robert
--
Robert Ivanek
Postdoctoral Fellow Schuebeler Group
Friedrich Miescher Institute
Maulbeerstrasse 66
4058 Basel / Switzerland
Office phone: +41 61 697 6100
Henrik Bengtsson
On Wed, May 26, 2010 at 3:24 PM, Ivanek, Robert <robert...@fmi.ch> wrote:
> Dear Sir or Madam,
>
> I would like to analyse the copy number variation data from Affymetrix
> Mouse Diversity Array. I have not found any information on your website
> about this particular array.
I have created page for this:
http://aroma-project.org/chipTypes/MOUSEDIVm520650
> I have tried to build the annotation files
> which are required by aroma but without success. I have few questions
> regarding that:
>
> 1: Is aroma.affymetrix able to analyse the "Mouse Diversity Array" ?
Yes, because there should be no reason why it shouldn't - it uses a
standard CDF etc. As you've noted, UGP (and UFL) files have not been
created by anyone yet.
For CN analysis, at least the UGP (genome positions) annotation data
file needs to be there.
I had a look at the MOUSEDIVm520650.na30.annot.csv file. The line
containing column names, that is:
"Probe Set ID","dbSNP RS ID","Chromosome","Physical
Position","Strand","Cytoband","Allele A","Allele B","Associated
Gene","Genetic Map","Fragment Enzyme Type Length Start Stop",
contains a trailing comma (,) that shouldn't be there ("file format
error"). This cause R to think there should be 12 and not 11 columns
in the data set. Open the file in an editor and remove that trailing
comma and any whitespace after "Fragment Enzyme Type Length Start
Stop". Then save the file. That should solve the problem.
The other CSV file - MOUSEDIVm520650.cn.na30.annot.csv - does not have
this problem.
>
>
> 3. I tried it also by using the "manual" approach using the
> tab=delimited file, however it seems to me that the mitochondria probes
> are skipped (NA values in ugp[,1] but valid values in ugp[,2]).
The Affymetrix NetAffx CSV files use s "M" for the mitochondria
chromosome. In aroma we encode this by integer 25.
> Another
> problem is that some positions for other chromosomes are not loaded in
> properly (valid values in ugp[,1] but NA values in ugp[,2]).
You don't show how you read the data "manually", so it is hard to say
what you are doing wrong here. But note that there are quite a few
arguments in read.table() that you need to set correctly in order to
read Affymetrix NetAffx CSV files (it doesn't make easier that
Affymetrix changes the file format once in a while and have stray
erroneous symbols such as the above comma).
Also, search our forum for 'MOUSEDIVm520650', because about a year ago
David Rosenberg disscussed this chip type and I think he did create
various annotation data files for the chip type. This was before the
chip type was publicly announced by Affymetrix.
/Henrik
> --
> 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 with website http://www.aroma-project.org/.
> To post to this group, send email to aroma-af...@googlegroups.com
> To unsubscribe and other options, go to http://www.aroma-project.org/forum/
>
Ivanek, Robert
Thanks for the answer and also the ACS file.
I have one more question regarding the UFL file generation.
I tried it by using the NettAffx and I got the following error:
R> ufl <- AromaUflFile$allocateFromCdf(cdf, nbrOfEnzymes=2,
tags=c("na30", "RI20100601"))
R> csv <- AffymetrixNetAffxCsvFile$byChipType(chipType, tags=".na30");
R> units <- importFrom(ufl, csv);
Warning messages:
1: In updateDataColumn.AromaTabularBinaryFile(this, .con = con, .hdr =
hdr, :
Values to be assigned were out of range [-32767,32768] and therefore
censored to fit the range.
2: In updateDataColumn.AromaTabularBinaryFile(this, .con = con, .hdr =
hdr, :
Values to be assigned were out of range [-32767,32768] and therefore
censored to fit the range.
R> csv <- AffymetrixNetAffxCsvFile$byChipType(chipType,
tags=".cn.na30");
R> units <- importFrom(ufl, csv);
Warning messages:
1: In updateDataColumn.AromaTabularBinaryFile(this, .con = con, .hdr =
hdr, :
Values to be assigned were out of range [-32767,32768] and therefore
censored to fit the range.
2: In updateDataColumn.AromaTabularBinaryFile(this, .con = con, .hdr =
hdr, :
Values to be assigned were out of range [-32767,32768] and therefore
censored to fit the range.
And the summary produce the following
R> summary(ufl)
length length.02
Min. :-32768 Min. :-32768
1st Qu.: 614 1st Qu.: 541
Median : 1146 Median : 997
Mean : 1601 Mean : 1466
3rd Qu.: 2195 3rd Qu.: 2000
Max. : 22095 Max. : 30002
NA's :230775 NA's :230775
Would you be so kind and build also the UFL and UGP files?
Best Regards
Robert
Ivanek, Robert
I was a little bit investigating the error and I found out that some
of the fragments reported in NetAffx files are really long.
Why they got a negative value of -32768 and not a positive one?
Robert
Henrik Bengtsson
On Wed, Jun 2, 2010 at 11:16 AM, Ivanek, Robert <robert...@fmi.ch> wrote:
> HI Henrik,
>
> I was a little bit investigating the error and I found out that some
> of the fragments reported in NetAffx files are really long.
> Why they got a negative value of -32768 and not a positive one?
Thanks for reporting. It turns out to be a bug in aroma.core causing
it to censor values into [-32767,32768], whereas it should have been
[-32768,32767]. Thus, the fragment lengths that are too large where
written as 32768, which when read back became -32768 (that's how
signed integers loops around when output of range). That should have
been written as 32767.
I have fixed this in the next release of aroma.core. Until that is
released, you can install a patch as explained in:
http://aroma-project.org/howtos/updateOrPatch
With the patch, you will get correct censoring and more informative
warnings, e.g.
Warning messages:
1: In updateDataColumn.AromaTabularBinaryFile(this, .con = con, .hdr = hdr, :
33 values to be assigned were out of range [-32768,32767] and
therefore censored to fit the range. Of these, 33 values in
[35102,655381] were too large.
2: In updateDataColumn.AromaTabularBinaryFile(this, .con = con, .hdr = hdr, :
21 values to be assigned were out of range [-32768,32767] and
therefore censored to fit the range. Of these, 21 values in
[50496,56758] were too large.
About the very large fragment lengths: My guess is that they are
"theoretical" fragments lengths. After running the PCR in the assay,
very long fragments are not amplified and hence filtered out. For the
specific enzyme, you should not get any hybrization signal for very
long fragments. It is possible that you have signal from the cuts of
the other enzyme. Maybe someone else has a better explanation of why
they are so long and still on the array? You could also drop a
message on the Affymetrix forums and ask.
/Henrik
> You received this message because you are subscribed to the Google Groups "aroma.affymetrix" group with website http://www.aroma-project.org/.
Ivanek, Robert
I think you are right, the fragment sizes are theoretical ones. I
would guess that the reason why also the long fragments are reported
is because the same SNP is present in short fragment produced by the
other enzyme.
Thank you very much for the patch.
Would you mind to update the MOUSEDIVm520650 chipType page and add
there the UGP and UFL files?
Best Regards
Robert
>
> read more »
Henrik Bengtsson
> Hi Henrik,
>
> I think you are right, the fragment sizes are theoretical ones. I
> would guess that the reason why also the long fragments are reported
> is because the same SNP is present in short fragment produced by the
> other enzyme.
>
> Thank you very much for the patch.
>
> Would you mind to update the MOUSEDIVm520650 chipType page and add
> there the UGP and UFL files?
Ideally users contribute with UGP and UFL too, though this time I've
done it since I've already done most of the work. Please compare to
what you got when you did.
/Henrik
Ivanek, Robert
one more question to this array.
During the analysis (Step 4 - Normalization for PCR fragment-length
effect) I am getting the following error:
[2010-06-09 11:51:49] Exception: Cannot fit normalization function to
enzyme, because there are no (finite) data points that are unique to
this enzyme: 1
at throw(Exception(...))
at throw.default("Cannot fit normalization function to enzyme,
because there are no (finite) data points that are unique to this
enzyme: ", ee)
at throw("Cannot fit normalization function to enzyme, because there
are no (finite) data points that are unique to this enzyme: ", ee)
at normalizeFragmentLength.default(y, fragmentLengths = fl,
targetFcns = targetFcns, subsetToFit = subset, onMissing =
onMissing, ...)
at normalizeFragmentLength(y, fragmentLengths = fl, targetFcns =
targetFcns, subsetToFit = subset, onMissing = onMissing, ...)
at process.FragmentLengthNormalization(fln, verbose = verbose)
I think that this is caused by the fact that the UFL file contains
also "theoretical" values for all fragments and therefore there is no
single SNP/CN probe annoatated as produced by only one enzyme.
What do you think about filtering out fragments bigger than 5kb (or
even smaller)? They should not be amplified anyway during the PCR.
Best Regards
Robert
On Jun 3, 3:22 pm, Henrik Bengtsson <h...@stat.berkeley.edu> wrote:
>
> read more »
Henrik Bengtsson
On Wed, Jun 9, 2010 at 11:57 AM, Ivanek, Robert <robert...@fmi.ch> wrote:
> Hi Henrik,
>
> one more question to this array.
> During the analysis (Step 4 - Normalization for PCR fragment-length
> effect) I am getting the following error:
>
> [2010-06-09 11:51:49] Exception: Cannot fit normalization function to
> enzyme, because there are no (finite) data points that are unique to
> this enzyme: 1
> at throw(Exception(...))
> at throw.default("Cannot fit normalization function to enzyme,
> because there are no (finite) data points that are unique to this
> enzyme: ", ee)
> at throw("Cannot fit normalization function to enzyme, because there
> are no (finite) data points that are unique to this enzyme: ", ee)
> at normalizeFragmentLength.default(y, fragmentLengths = fl,
> targetFcns = targetFcns, subsetToFit = subset, onMissing =
> onMissing, ...)
> at normalizeFragmentLength(y, fragmentLengths = fl, targetFcns =
> targetFcns, subsetToFit = subset, onMissing = onMissing, ...)
> at process.FragmentLengthNormalization(fln, verbose = verbose)
>
> I think that this is caused by the fact that the UFL file contains
> also "theoretical" values for all fragments and therefore there is no
> single SNP/CN probe annoatated as produced by only one enzyme.
You are completely right here; the algorithm for estimating PCR
fragment-length normalization effects assumes that for each
restriction enzyme (if more than one is used as here), there exists at
least some loci (SNPs or CN units) that are only for that enzyme. It
is these data points that are used to estimate the fragment-length
effect per enzyme, which is then used for correcting for the effect.
More precisely, this is Step (1) in the algorithm described in Section
2.2.4 'Normalization for fragment-length effects' in the CRMAv2 paper
(Bengtsson, Wirapati et al. 2009;
http://www.aroma-project.org/publications/).
Now looking at the UFL information imported from the NetAffx CSV files:
library("aroma.affymetrix");
cdf <- AffymetrixCdfFile$byChipType("MOUSEDIVm520650");
ufl <- getAromaUflFile(cdf);
print(ufl);
AromaUflFile:
Name: MOUSEDIVm520650
Tags: na30,mm9,HB20190603
Full name: MOUSEDIVm520650,na30,mm9,HB20190603
Pathname: annotationData/chipTypes/MOUSEDIVm520650/MOUSEDIVm520650,na30,mm9,HB20190603.ufl
File size: 9.38 MB (9835569 bytes)
RAM: 0.00 MB
Number of data rows: 2458697
File format: v1
Dimensions: 2458697x2
Column classes: integer, integer
Number of bytes per column: 2, 2
Footer: <createdOn>20100603 05:11:33 PDT</createdOn><platform>Affymetrix</platfo
rm><chipType>MOUSEDIVm520650</chipType><createdBy><fullname>Henrik Bengtsson</fu
llname><email>h...@aroma-project.org</email></createdBy><srcFiles><srcFile1><filen
ame>MOUSEDIVm520650.CDF</filename><filesize>668362285</filesize><checksum>1d83d7
dde5e4816f4be315f6863f6a7c</checksum></srcFile1><srcFile2><filename>MOUSEDIVm520
650.na30.annot.csv</filename><filesize>358437027</filesize><checksum>e3ef3524595
4d5a08ae299b95666fc3a</checksum></srcFile2><srcFile3><filename>MOUSEDIVm520650.c
n.na30.annot.csv</filename><filesize>795251151</filesize><checksum>5b5bbdb60ed44
3b2616358beea06ddaf</checksum></srcFile3></srcFiles>
Chip type: MOUSEDIVm520650
Platform: Affymetrix
print(getChecksum(ufl));
[1] "31cc2be0708c34a977f456dfb4b4d3b9"
hasEnzyme1 <- is.finite(ufl[,1,drop=TRUE]);
hasEnzyme2 <- is.finite(ufl[,2,drop=TRUE]);
hasEnzyme1Only <- (hasEnzyme1 & !hasEnzyme2);
hasEnzyme2Only <- (!hasEnzyme1 & hasEnzyme2);
hasEither <- (hasEnzyme1 | hasEnzyme2);
print(summary(hasEither));
Mode FALSE TRUE NA's
logical 230775 2227922 0
print(summary(hasEnzyme1Only));
Mode FALSE NA's
logical 2458697 0
print(summary(hasEnzyme2Only));
Mode FALSE NA's
logical 2458697 0
So, yes, neither enzyme has data points that are unique to that
enzyme. This is a problem with the annotation data - not the CEL
file data.
Querying say probesets "JAX00289207" the online Affymetrix NetAffx
Analysis Center:
and the click on the small 'Details' label (not the check box) you get to:
https://www.affymetrix.com/analysis/netaffx/mappingfullrecord.affx?pk=MOUSEDIVm520650:JAX00289207
and confirms that the fragment-length information is:
Enzyme Name, Length
NspI, 11635
StyI, 406
Compare to:
unit <- indexOf(cdf, "JAX00289207");
print(ufl[unit,]);
length length.02
1 11635 406
> What do you think about filtering out fragments bigger than 5kb (or
> even smaller)? They should not be amplified anyway during the PCR
That might be a solution - though a bit ad hoc. In order to get this
correct, I would spend the time to carefully do the following:
1) Check what the Affymetrix's Assay documentation says about what
range of fragment lengths the PCR amplification is supposed to keep
and what are supposed to be filtered out. You should be able to fine
that documentation via the Affymetrix Product Page for
MOUSEDIVm520650. You find links to the latter at
http://www.aroma-project.org/chipTypes/MOUSEDIVm520650
2) Ask Affymetrix about the fragment-length information available in
the NetAffx CSV files (and the above). The NetAffx group are quite
responsive and are happy to get feedback so they can correct any
mistakes. If you are quick enough, they might be able to correct it
before the NetAffx Release 31 scheduled for "June/July". You can post
a message to one of the forums at Affymetrix Scientific Community
Forums [https://www.affymetrix.com/community/forums/index.jspa]. You
may also want send a message to Affymetrix Support
[www.affymetrix.com/support/] - in my experience they are also very
responsive. I would point them also to the above NetAffx query
results.
It's worth getting this right once and for all. Feel free to update
us here on the list.
Cheers,
Henrik
Henrik Bengtsson
On Thu, Jun 16, 2011 at 10:58 AM, Anguraj Sadanandam
<anguraj.s...@epfl.ch> wrote:
> Hello,
>
> I am trying to process several (~800) .CEL files of Hug133a affymetrix and
> having trouble processing due to memory error. I tried to look through aroma
> affymetrix code to do this but not fruitful.
Exactly what did you try and what did not work?
The Aroma framework can process any number of arrays. You'll find vignettes at:
http://aroma-project.org/vignettes/
and you may also be interested in the one-line doRMA() command:
http://aroma-project.org/blocks/doRMA
/Henrik
>
> Could you please point me to the link to memory efficient RMA analysis?
>
> Thanks,
>
> Anguraj
>
>
Pierre Neuvial
---------- Forwarded message ----------
From: Pierre Neuvial <pie...@stat.berkeley.edu>
Date: Wed, Nov 16, 2011 at 10:16 PM
Subject: Re: [aroma.affymetrix] Processing multiple Affymetrix arrays
To: Anguraj Sadanandam <anguraj.s...@epfl.ch>
Hi Anguraj,
On Wed, Nov 16, 2011 at 4:13 PM, Anguraj Sadanandam
<anguraj.s...@epfl.ch> wrote:
>
> Hi Henrik,
>
> I am interested in processing about 228 Affy SNP 6.0 array samples. I used the following code (attached below). Though it was working well, I was trying to get the raw data through a loop (as you can see in the code below). I was executing the loop for 228 times for that many samples (may be not a good idea). The code stopped somewhere around 162. This happened 2 times but not exactly at 162. May be the loop is not the best way. Could you please suggest the better way? I like to increase the number of samples up to 1070 and is it good idea to do it by batch or do it all together? I have the normal samples for normalization.
>
> My sessioninfo() is below.
>
> Thanks,
>
> Anguraj
>
>
>
> dataSet <- "Panc";
> chipType <- "GenomeWideSNP_6,Full";
> expNbrOfArrays <- 228;
>
>
> library("aroma.affymetrix");
> verbose <- Arguments$getVerbose(-8, timestamp=TRUE);
>
> cdf <- AffymetrixCdfFile$byChipType(chipType);
> print(cdf);
>
> # Assert that all annotation data files exist
> gi <- getGenomeInformation(cdf);
> print(gi);
>
> si <- getSnpInformation(cdf);
> print(si);
>
> acs <- AromaCellSequenceFile$byChipType(getChipType(cdf, fullname=FALSE));
> print(acs);
>
> # Setup data set
> cs <- AffymetrixCelSet$byName(dataSet, cdf=cdf);
> print(cs);
> # Sanity check
> stopifnot(nbrOfArrays(cs) == expNbrOfArrays);
>
>
> # Step 1: Crosstalk calibration
> acc <- AllelicCrosstalkCalibration(cs);
> print(acc);
> csC <- process(acc, verbose=verbose);
> print(csC);
>
>
> # Step 2 - Normalization for nucleotide-position probe sequence effects
> bpn <- BasePositionNormalization(csC, target="zero");
> print(bpn);
> csN <- process(bpn, verbose=verbose);
> print(csN);
>
>
> # Step 3 - Probe summarization
> plm <- AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=TRUE);
> print(plm);
>
> if (length(findUnitsTodo(plm)) > 0) {
> # Fit CN probes quickly (~5-10s/array + some overhead)
> units <- fitCnProbes(plm, verbose=verbose);
> str(units);
> # int [1:945826] 935590 935591 935592 935593 935594 935595 ...
> # Fit remaining units, i.e. SNPs (~5-10min/array)
> units <- fit(plm, verbose=verbose);
> str(units);
> }
>
> # Get the estimated chip effects
> ces <- getChipEffectSet(plm);
> print(ces);
>
>
> # Step 4 - Normalization for PCR fragment-length effects
> fln <- FragmentLengthNormalization(ces);
> print(fln);
> cesN <- process(fln, verbose=verbose);
> print(cesN);
>
FYI, you could replace Step 1 through 4 by
cesN <- doCRMAv2(cs)
see http://www.aroma-project.org/blocks/doCRMAv2
>
> # Sanity check
> stopifnot(nbrOfArrays(cesN) == expNbrOfArrays);
>
> # Identify the index of all normal samples by mapping to names in a text file (
> or whatever you wish)
> normalNames <- readLines("normal.txt");
> #normals <- indexOf(cesN, normalNames);
> normals <- match(normalNames, getNames(cesN))
> str(normals);
> stopifnot(all(is.finite(normals))); # Sanity check
>
> # Extract pool of reference samples (normals only)
> cesR <- extract(cesN, normals);
> print(cesR);
>
> # Calculate the robust average across this pool
> ceR <- getAverageFile(cesR, verbose=verbose);
> print(ceR);
>
> #cesR <- cesNHapMap;
> #print(cesR);
> ceR <- getAverageFile(cesR, verbose=verbose);
> print(ceR);
>
> # Set up the CBS model with ceR as the reference for all samples
> cbs <- CbsModel(cesN, ceR);
> print(cbs);
>
> # Fit the segmentation
> #fit(cbs, verbose=verbose);
>
> gi <- getGenomeInformation(cdf)
> print(gi)
>
> l<-NULL
> for(i in 1:23) {
> units <- getUnitsOnChromosome(gi, chromosome=i);
> pos <- getPositions(gi, units=units);
> l[i]<-length(pos)
> }
>
> length(l)
>
> l_sum<-sum(l)
>
> raw<-NULL
>
> for(i in 1:228) {
> cat("\n",i,"\n")
> for(j in 1:23) {
> rawCNs <- extractRawCopyNumbers(cbs, array=i, chromosome=j)
> raw<-rbind(raw,as.data.frame(rawCNs))
> cat(j,"\n")
> }
> }
>
Here you are trying to concatenate copy number estimates for all
samples in a huge data.frame.
Remember that the SNP 6 array has roughly 2 million loci ! After 162
iterations, this data.frame has about 324 million lines !!!
It probably breaks down when your computer runs out of memory.
I see two problems here:
1) I don't think it's a good idea to try to create a text file with
228*2e6=456 million lines. Why don't you create only one file per
sample ? That would also be much easier to manipulate afterwards.
2) In case you'd really want to do that (again, I really don't see why
if you have that many samples), you should write the results to file
for each chromosome directly, using write.table(..., append=TRUE)
within the loop on index "j". This way you would avoid using so much
memory, and also loosing a lot of time with this "rbind" which
re-creates a bigger and bigger data.frame at each iteration.
And my suggestion is: look at writeDataFrame: it should give you what
you need, see:
http://www.aroma-project.org/howtos/writeDataFrame
Cheers,
Pierre
> write.table(raw,"all-data-PANC-40-normal.txt",sep="\t")
>
> raw_1<-matrix(raw[,2],nrow=l_sum,ncol=228)
> raw_2<-cbind(raw[1:l_sum,1],raw_1)
>
> write.table(raw_1,"all-data-matrix-Panc-40-normal.txt", sep="\t")
> write.table(raw_2,"all-data-matrix-probes-PANC-40-normal.txt", sep="\t")
>
>
> R version 2.13.0 (2011-04-13)
> Platform: x86_64-redhat-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
>
>
>
>