I am still not able to get same as you for units (so if I continue
doing units2 <- fit(plm, verbose=-10) I am not getting the right thing
either).
I ran everything again and here is the output of the last part:
> plm <- AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=TRUE)
> units1 <- fitCnProbes(plm, verbose=verbose);
Estimating single-probe CN units...
Getting chip-effect set...
Getting chip-effect set from data set...
Getting chip-effect set from data set...done
Getting chip-effect set...done
Identifying CN units...
int [1:945826] 935590 935591 935592 935593 935594 935595 935596 935597
935598 935599 ...
Identifying CN units...done
Identifying subset to fit...
Identifying non-fitted units in chip-effect file...
Pathname:
plmData/test,ACC,ra,-XY,BPN,-XY,AVG,A+B/GenomeWideSNP_6/91228050T,chipEffects.CEL
Found indices cached on file
Reading data for these 1881415 cells...
Reading data for these 1881415 cells...done
Looking for stdvs <= 0 indicating non-estimated units:
int [1:1881415] 1 2 3 4 5 6 7 8 9 10 ...
Identifying non-fitted units in chip-effect file...done
int [1:945826] 935590 935591 935592 935593 935594 935595 935596 935597
935598 935599 ...
Identifying subset to fit...done
Identifying cell indices for CN units...
Identifying cell indices for CN units...done
Keeping only single-cell units...
Single-cell units:
int [1:945826] 935590 935591 935592 935593 935594 935595 935596 935597
935598 935599 ...
Cell indices:
int [1:945826] 3435872 3403913 422133 2456228 4388533 5792832 3788106
3660159 3524469 3264141 ...
Keeping only single-cell units...done
Identifying cell indices for estimates...
int [1:945826] 1876054 1876055 1876056 1876057 1876058 1876059 1876060
1876061 1876062 1876063 ...
Identifying cell indices for estimates...done
Fitting 1 arrays...
Array #1 ('91228050T') of 1...
Reading signals...
num [1:945826] 6448 7260 7878 6151 12230 ...
Reading signals...done
Transforming signals to estimates...
'data.frame': 945826 obs. of 4 variables:
$ cell : int 1876054 1876055 1876056 1876057 1876058 1876059
1876060 1876061 1876062 1876063 ...
$ theta : num 6448 7260 7878 6151 12230 ...
$ sdTheta : num 1.49e-08 1.49e-08 1.49e-08 1.49e-08 1.49e-08 ...
$ outliers: logi FALSE FALSE FALSE FALSE FALSE FALSE ...
Transforming signals to estimates...done
Writing estimates...
Storing flat data to file...
Updating file...
Updating file...done
Storing flat data to file...done
Writing estimates...done
Array #1 ('91228050T') of 1...done
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 483738 25.9 1733216 92.6 12913475 689.7
Vcells 28397870 216.7 83594476 637.8 130588565 996.4
Fitting 1 arrays...done
Estimating single-probe CN units...done
> str(units1);
int [1:945826] 935590 935591 935592 935593 935594 935595 935596 935597
935598 935599 ...
> units <- findUnitsTodo(plm);
> str(units);
int [1:935589] 1 2 3 4 5 6 7 8 9 10 ...
On Monday, 26 November 2012 23:08:47 UTC+2, Henrik Bengtsson wrote:
>
> On Mon, Nov 26, 2012 at 12:46 PM, Ioana Cutcutache
> <ioana.cu...@gmail.com <javascript:>> wrote:
> > Ok, so the first part ran without any errors, output below:
> >
> >> print(csN);
> > AffymetrixCelSet:
> > Name: test
> > Tags: ACC,ra,-XY,BPN,-XY
> > Path: probeData/test,ACC,ra,-XY,BPN,-XY/GenomeWideSNP_6
> > Platform: Affymetrix
> > Chip type: GenomeWideSNP_6,Full
> > Number of arrays: 1
> > Names: 91228050T [1]
> > Time period: 2011-07-19 05:08:51 -- 2011-07-19 05:08:51
> > Total file size: 65.83MB
> > RAM: 0.01MB
> >> print(getFullNames(csN));
> > [1] "91228050T"
> >> cdf2 <- getCdf(csN);
> >> stopifnot(identical(getChecksum(cdf2), getChecksum(cdf))); # IMPORTANT
> >>
> >
> > However I am not getting same as you for the units:
> >
> >> plm <- AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=TRUE);
> >> units <- findUnitsTodo(plm);
> >> str(units);
> > int [1:1881415] 1 2 3 4 5 6 7 8 9 10 ...
>
> Alright (but has nothing to do with the error you first reported). If so,
> do:
>
> units1 <- fitCnProbes(plm, verbose=verbose);
> str(units1);
> ## int [1:945826] 935590 935591 935592 935593 935594 935595 935596
> 935597 935598 935599 ...
>
> and you should get the same as:
>
> units <- findUnitsTodo(plm);
> str(units);
> ## int [1:934968] 622 623 624 625 626 627 628 629 630 631 ...
>
> If so, the try with
>
> units2 <- fit(plm, verbose=-10);
>
> and compare to the output I already pasted in my previous message (below).
>
> /Henrik
>
> >
> >
> > On Monday, 26 November 2012 21:15:32 UTC+2, Henrik Bengtsson wrote:
> >>
> >> Hmm... first run this in a fresh R session:
> >>
> >> library("aroma.affymetrix");
> >> cdf <- AffymetrixCdfFile$byChipType("GenomeWideSNP_6", tags="Full");
> >> csR <- AffymetrixCelSet$byName("test", cdf=cdf);
> >> acc <- AllelicCrosstalkCalibration(csR, model="CRMAv2");
> >> csC <- process(acc, verbose=verbose);
> >> bpn <- BasePositionNormalization(csC, target="zero");
> >> csN <- process(bpn, verbose=verbose);
> >>
> >> then what's the output of:
> >> print(csN);
> >> print(getFullNames(csN));
> >> cdf2 <- getCdf(csN);
> >> stopifnot(identical(getChecksum(cdf2), getChecksum(cdf))); # IMPORTANT
> >>
> >> If you get an error above, what is the output of print(cdf2) and
> >> getChecksum(cdf2)? If not, continue:
> >>
> >> plm <- AvgCnPlm(csN, mergeStrands=TRUE, combineAlleles=TRUE);
> >> units <- findUnitsTodo(plm);
> >> str(units);
> >> ## int [1:934968] 622 623 624 625 626 627 628 629 630 631 ...
> >>
> >> Do you get the same? If so, continue:
> >>
> >> units2 <- fit(plm, verbose=-10);
> >>
> >> You should get something like the following, except from the
> >> dataset/sample names (if not, let me know):
> >>
> >> ## Fitting model of class AvgCnPlm...
> >> ## AvgCnPlm:
> >> ## Data set: GSE13372
> >> ## Chip type: GenomeWideSNP_6,Full
> >> ## Input tags: testset,ACC,ra,-XY,BPN,-XY
> >> ## Output tags: testset,ACC,ra,-XY,BPN,-XY,AVG,A+B
> >> ## Parameters: {probeModel: chr "pm", shift: num 0, flavor: chr
> >> "median", mergeStrands: logi TRUE, combineAlleles: logi TRUE}
> >> ## Path:
> >> plmData/GSE13372,testset,ACC,ra,-XY,BPN,-XY,AVG,A+B/GenomeWideSNP_6
> >> ## RAM: 0.01MB
> >> ## Identifying non-estimated units...
> >> ## Identifying non-fitted units in chip-effect file...
> >> ## Pathname:
> >>
> >>
> plmData/GSE13372,testset,ACC,ra,-XY,BPN,-XY,AVG,A+B/GenomeWideSNP_6/GSM337708,chipEffects.CEL
>
>
> >> ## Found indices cached on file
> >> ## Reading data for these 1881415 cells...
> >> ## Reading data for these 1881415 cells...done
> >> ## Looking for stdvs <= 0 indicating non-estimated units:
> >> ## int [1:934968] 622 623 624 625 626 627 628 629 630 631 ...
> >> ## Identifying non-fitted units in chip-effect file...done
> >> ## Identifying non-estimated units...done
> >> ## Getting model fit for 934968 units.
> >> ## Identifying unit types:...
> >> ## Units:
> >> ## int [1:934968] 622 623 624 625 626 627 628 629 630 631 ...
> >> ## Unit types:
> >> ## atomic [1:934968] 2 2 2 2 2 2 2 2 2 2 ...
> >> ## - attr(*, "types")= Named int [1:8] 0 1 2 3 4 5 6 7
> >> ## ..- attr(*, "names")= chr [1:8] "unknown" "expression"
> >> "genotyping" "resequencing" ...
> >> ## unitTypes
> >> ## 2
> >> ## 934968
> >> ## Unique unit types:
> >> ## genotyping
> >> ## 2
> >> ## Identifying unit types:...done
> >> ## Setting up parameter sets...
> >> ## Creating CEL file...
> >> ## Chip type: GenomeWideSNP_6,Full
> >> ## Pathname:
> >>
> >>
> plmData/GSE13372,testset,ACC,ra,-XY,BPN,-XY,AVG,A+B/GenomeWideSNP_6/probeAffinities.CEL
>
>
> >> ## Returning already existing file.
> >> ## Creating CEL file...done
> >> ## used (Mb) gc trigger (Mb) max used (Mb)
> >> ## Ncells 549305 29.4 1606180 85.8 13015638 695.2
> >> ## Vcells 7692698 58.7 23004541 175.6 40154683 306.4
> >> ## Setting up parameter sets...done
> >> ## Fitting AvgCnPlm for each unit type separately...
> >> ## Unit types:
> >> ## genotyping
> >> ## 2
> >> ## Unit type #1 ('genotyping') of 1...
> >> ## Unit type: genotyping (code=2)
> >> ## Number of units of this type: 934968
> >> ## int [1:934968] 622 623 624 625 626 627 628 629 630 631 ...
> >> ## Fitting the model by unit dimensions (at least for the large
> >> classes)...
> >> ## Grouping units into equivalent (unit,group,cell) dimensions...
> >> [...]
> >>
> >> What is particular interesting are the outputs of your unit types and
> >> unit indices.
> >>
> >> /Henrik
> >>
> >> On Mon, Nov 26, 2012 at 10:45 AM, Ioana Cutcutache
> >> <ioana.cu...@gmail.com> wrote:
> >> >> library("aroma.affymetrix")
> >> >> log <- verbose <- Arguments$getVerbose(-8, timestamp=TRUE)
> >> >> options(digits=4)
> >> >> cdf <- AffymetrixCdfFile$byChipType("GenomeWideSNP_6", tags="Full")
> >> >> print(cdf)
> >> >> gi <- getGenomeInformation(cdf)
> >> >> print(gi)
> >> >> si <- getSnpInformation(cdf)
> >> >> print(si)
> >> >> acs <- AromaCellSequenceFile$byChipType(getChipType(cdf,
> >> >> fullname=FALSE))
> >> >> print(acs)
> >> >> cdf <- AffymetrixCdfFile$byChipType("GenomeWideSNP_6", tags="Full")
> >> >> csR <- AffymetrixCelSet$byName("test", cdf=cdf)
> >> >> print(csR)
> >> >> acc <- AllelicCrosstalkCalibration(csR, model="CRMAv2")
> >> >> print(acc)
> >> >> csC <- process(acc, verbose=verbose)
> >> >> print(csC)
> >> >> bpn <- BasePositionNormalization(csC, target="zero")
> >> >> print(bpn)
> >> >> csN <- process(bpn, verbose=verbose)
> >> >> print(csN)
> >
> > --
> > 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<javascript:>
> > To unsubscribe and other options, go to
> http://www.aroma-project.org/forum/
>
--
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-affymetrix@googlegroups.com
To unsubscribe and other options, go to http://www.aroma-project.org/forum/
