[BioC] seeking help on cghMCR pleeeease
Nathalie Conte
nac at sanger.ac.uk
Mon Dec 12 15:17:58 CET 2011
HI,
I am using cghMCR and I would need some advice on the software
behaviour. I don't seem to get any different results when I change the
gapAllowed parameter.
I have used DNAcopy to create a segmented data for all my samples. Then
I choose to use cghMCR to get common alterations between my samples.
These are all arguments from the package doc(latest version):
segments is a data frame extracted from the "output" element
of the object
returned by segment of the package DNAcopy or getSegments
gapAllowed is an integer specifying low threshold of base
pair number to
separate two adjacent segments, belower which the two
segments will be joined
as an altered span
alteredLow is a positive number between 0 and 1 specifying
the lower resh-
old percential value. Only segments with values falling
below this threshold are
considered as altered span
alteredHigh is a positive number between 0 and 1 specifying
the upper resh-
old percential value. Only segments with values falling over
this threshold are
considered as altered span
recurrence is an integer between 1 and 100 that specifies
the rate of occur-
rence for a gain or loss that are observed across sample.
Only gains/losses with
ocurrence rate grater than the threshold values are
declared as MCRs
spanLimit is an integer that defines the leangh of altered
spans that can be
considered as locus. It is not of any use at this time
thresholdType is a character string that can be either
"quantile" or "value"
indicating wether alteredLow or alteredHigh is quantial or
actual value
In my analysis, I have used cghMCR with threshold value log2R of low
-0.25 and high =0.25 , a spanLimit of 2.10^7 and I have tried several
gapAllowed values, 5, 500 , 5000 and 50000 .Each time I used the MCR
function to identify the minimum regions. I have put an example of each
results only on 25 lines of chromosome 4 to avoid massive sized files
(see attached) , but basically whatever the gapAllowed size i get the
same MCRs at the end for all postitions. I would expect this to vary as
segments should be fused together differently given this parameter..
Could you please help with this and advise on what I might be doing wrong?
Is spanLimit any use? in the doc, It is written "It is not of any use
at this time"??
Another point, the gapAllowed is specified " is an integer specifying
low threshold of base pair number " is the unit in base pair number, in
several examples it seems that the units are in kb???
thanks a lot ,
code for gapAllowed=5
##get the cghMCR function with these parameters
cghmcr0.25T_5k_2_20M_sdundo1.5 <- cghMCR(sdundo.segData_order1.5,
gapAllowed = 5, alteredLow = -0.25,alteredHigh = 0.25,
spanLimit=20000000,recurrence=2,thresholdType=c("value"))
##identify the MCRs
mcrs0.25T_5k_2_20M_sdundo1.5<-MCR(cghmcr0.25T_5k_2_20M_sdundo1.5)
mcrs0.25T_5k_2_20M_sdundo.bind1.5 <-cbind (
mcrs0.25T_5k_2_20M_sdundo1.5[, c ( "chromosome", "status", "loc.start",
"loc.end","mcr.start", "mcr.end","samples" ) ],z[1:9171] )
write.table(mcrs0.25T_5k_2_20M_sdundo.bind1.5,
file="PT_mcrs0.25T_5k_2_20M_sdundo.bind1.5.txt", sep="\t", row.names=F)
PT_mcrs0.25T_5k_2_20M_sdundo.bind1.5=read.table(file="PT_mcrs0.25T_5k_2_20M_sdundo.bind1.5.txt",
sep="\t", header=T)
##only 25 lines of chromsome 4
test.5k=head(PT_mcrs0.25T_5k_2_20M_sdundo.bind1.5[PT_mcrs0.25T_5k_2_20M_sdundo.bind1.5==4,],25)
##attached data
write.table(test.5k, file="test.5k.txt", sep="\t")
code for gapAllowed=500
cghmcr0.25T_500k_2_20M_sdundo1.5 <- cghMCR(sdundo.segData_order1.5,
gapAllowed = 500, alteredLow = -0.25,alteredHigh = 0.25,
spanLimit=20000000,recurrence=2,thresholdType=c("value"))
mcrs0.25T_500k_2_20M_sdundo1.5<-MCR(cghmcr0.25T_500k_2_20M_sdundo1.5)
mcrs0.25T_500k_2_20M_sdundo.bind1.5 <-cbind (
mcrs0.25T_500k_2_20M_sdundo1.5[, c ( "chromosome", "status",
"loc.start", "loc.end","mcr.start", "mcr.end","samples" ) ],z[1:9171] )
write.table(mcrs0.25T_500k_2_20M_sdundo.bind1.5,
file="PT_mcrs0.25T_500k_2_20M_sdundo.bind1.5.txt", sep="\t", row.names=F)
PT_mcrs0.25T_500k_2_20M_sdundo.bind1.5=read.table(file="PT_mcrs0.25T_500k_2_20M_sdundo.bind1.5.txt",
sep="\t", header=T)
test.500k=head(PT_mcrs0.25T_500k_2_20M_sdundo.bind1.5[PT_mcrs0.25T_500k_2_20M_sdundo.bind1.5==4,],25)
write.table(test.500k, file="test.500k.txt", sep="\t")
code for gapAllowed=5000
cghmcr0.25T_5000k_2_20M_sdundo1.5 <- cghMCR(sdundo.segData_order1.5,
gapAllowed = 5000, alteredLow = -0.25,alteredHigh = 0.25,
spanLimit=20000000,recurrence=2,thresholdType=c("value"))
mcrs0.25T_5000k_2_20M_sdundo1.5<-MCR(cghmcr0.25T_5000k_2_20M_sdundo1.5)
mcrs0.25T_5000k_2_20M_sdundo.bind1.5 <-cbind (
mcrs0.25T_5000k_2_20M_sdundo1.5[, c ( "chromosome", "status",
"loc.start", "loc.end","mcr.start", "mcr.end","samples" ) ],z[1:9171] )
write.table(mcrs0.25T_5000k_2_20M_sdundo.bind1.5,
file="PT_mcrs0.25T_5000k_2_20M_sdundo.bind1.5.txt", sep="\t", row.names=F)
PT_mcrs0.25T_5000k_2_20M_sdundo.bind1.5=read.table(file="PT_mcrs0.25T_5000k_2_20M_sdundo.bind1.5.txt",
sep="\t", header=T)
test.5000k=head(PT_mcrs0.25T_5000k_2_20M_sdundo.bind1.5[PT_mcrs0.25T_5000k_2_20M_sdundo.bind1.5==4,],25)
write.table(test.5000k, file="test.5000k.txt", sep="\t")
code for gapAllowed=500000
cghmcr0.25T_500000k_2_20M_sdundo1.5 <- cghMCR(sdundo.segData_order1.5,
gapAllowed = 500000, alteredLow = -0.25,alteredHigh = 0.25,
spanLimit=20000000,recurrence=2,thresholdType=c("value"))
mcrs0.25T_500000k_2_20M_sdundo1.5<-MCR(cghmcr0.25T_500000k_2_20M_sdundo1.5)
mcrs0.25T_500000k_2_20M_sdundo.bind1.5 <-cbind (
mcrs0.25T_500000k_2_20M_sdundo1.5[, c ( "chromosome", "status",
"loc.start", "loc.end","mcr.start", "mcr.end","samples" ) ],z[1:9171] )
write.table(mcrs0.25T_500000k_2_20M_sdundo.bind1.5,
file="PT_mcrs0.25T_500000k_2_20M_sdundo.bind1.5.txt", sep="\t", row.names=F)
PT_mcrs0.25T_500000k_2_20M_sdundo.bind1.5=read.table(file="PT_mcrs0.25T_500000k_2_20M_sdundo.bind1.5.txt",
sep="\t", header=T)
test.500000k=head(PT_mcrs0.25T_500000k_2_20M_sdundo.bind1.5[PT_mcrs0.25T_500000k_2_20M_sdundo.bind1.5==4,],25)
write.table(test.500000k, file="test.500000k.txt", sep="\t")
sessioninfo()
R version 2.13.0 (2011-04-13)
Platform: x86_64-unknown-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
[5] LC_MONETARY=C LC_MESSAGES=C
[7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] tools stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] cghMCR_1.10.0 limma_3.4.5 CNTools_1.6.0 genefilter_1.30.0
[5] DNAcopy_1.22.1
loaded via a namespace (and not attached):
[1] annotate_1.26.1 AnnotationDbi_1.10.1 Biobase_2.8.0
[4] DBI_0.2-5 RSQLite_0.9-1 splines_2.13.0
[7] survival_2.35-8 xtable_1.6-0
--
The Wellcome Trust Sanger Institute is operated by Genome Research
Limited, a charity registered in England with number 1021457 and a
company registered in England with number 2742969, whose registered
office is 215 Euston Road, London, NW1 2BE.
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: test.5k.txt
URL: <https://stat.ethz.ch/pipermail/bioconductor/attachments/20111212/2c70ceda/attachment.txt>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: test.500k.txt
URL: <https://stat.ethz.ch/pipermail/bioconductor/attachments/20111212/2c70ceda/attachment-0001.txt>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: test.5000k.txt
URL: <https://stat.ethz.ch/pipermail/bioconductor/attachments/20111212/2c70ceda/attachment-0002.txt>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: test.500000k.txt
URL: <https://stat.ethz.ch/pipermail/bioconductor/attachments/20111212/2c70ceda/attachment-0003.txt>
More information about the Bioconductor
mailing list