Cutoff-based calling

[ilarischeinin]

callBins() uses CGHcall, which seems to expect the input to be a set of cancer samples, i.e. that there are a reasonable amount of aberrations in the set. When running with a single sample and/or non-cancer samples, it frequently fails.

I'm working on a simple cutoff-based calling that could be used instead in those cases.

It's in branch cutoff-calling of my fork.

[jiading]

Hi,

do you know why I still get such error msg, even if I try to use the cutoff method?
Thanks!

callBins(copyNumbersSegmented,method=c("cutoff"), cutoffs=c(-0.1, -0.1, 0.1, 0.1))
Error: Command CGHcall() returned the following error message:
Error in CGHcall(seg, organism = organism, ...): unused arguments (method = "cutoff", cutoffs = c(-0.1, -0.1, 0.1, 0.1))
Please contact maintainer of package CGHcall: Mark van de Wiel mark.vdwiel@vumc.nl

[ilarischeinin]

Your installed QDNAseq package version must be too old, as it does not recognize the method and cutoffs arguments. Cutoff-based calling was added in QDNAseq 1.8, which was part of Bioconductor 3.3. The current Bioconductor release is 3.4 (and includes QDNAseq 1.10).

By the way, if you do not want to distinguish between gains and amplifications, you should specify only one positive value for cutoffs, not the same value twice. And similarly for losses and homozygous deletions, and negative values. So, if you want to make calling according to these cutoffs: loss < -0.1 ≤ normal ≤ 0.1 < gain, you should use cutoffs=c(-0.1, 0.1).

[jiading]

I would also like to store the segments into output file. (e.g. the orange line in the final plots) Could you pls briefly type the command here? Thanks!

…

On Tue, Mar 7, 2017 at 3:14 PM, Ilari Scheinin ***@***.***> wrote: Your installed QDNAseq package version must be too old, as it does not recognize the method and cutoffs arguments. Cutoff-based calling was added in QDNAseq 1.8, which was part of Bioconductor 3.3. The current Bioconductor release is 3.4 (and includes QDNAseq 1.10). By the way, if you do not want to distinguish between gains and amplifications, you should specify only one positive value for cutoffs, not the same value twice. And similarly for losses and homozygous deletions, and negative values. So, if you want to make calling according to these cutoffs: loss < -0.1 < normal < 0.1 < gain, you should use cutoffs=c(-0.1, 0.1). — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#20 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAc-rQRo03lyePEwxuHWWAG0AsrpzDjVks5rjWZbgaJpZM4HHihe> .

[ilarischeinin]

exportBins(object, file="output.tsv", type="segments")

[jiading]

Since my binsize = 100kb exportBins(copyNumbersSegmented,file = OutFile, format="bed", filter=T,type=c("segments")) gives me: track name=".freeze1.4" description="segments" 1 900000 1000000 1:900001-1000000 -0.021 + 1 1000000 1100000 1:1000001-1100000 -0.021 + 1 1100000 1200000 1:1100001-1200000 -0.021 + 1 1200000 1300000 1:1200001-1300000 -0.021 + 1 1300000 1400000 1:1300001-1400000 -0.021 + 1 1400000 1500000 1:1400001-1500000 -0.021 + 1 1700000 1800000 1:1700001-1800000 -0.021 + 1 1900000 2000000 1:1900001-2000000 -0.021 + 1 2000000 2100000 1:2000001-2100000 -0.021 + What I want is the final segment line (the orange line). Is it possible to recover the data from that line? Thanks!

…

On Tue, Mar 7, 2017 at 3:59 PM, Ilari Scheinin ***@***.***> wrote: exportBins(object, file="output.tsv", type="segments") — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#20 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAc-rcjBr96KseyWqSBSz0k6Rl2oOmf9ks5rjXDGgaJpZM4HHihe> .

[ilarischeinin]

I know what you mean, but that's not implemented, at least not currently. It's something for the special case of single samples only, and not for data sets of two or more samples. Starting from the table you pasted above, you could use the rle() function and some basic R coding to achieve what you're after.

(An approach for larger data sets that includes this kind of dimensionality reduction could be for example something along the lines of:

library(CGHregions)

calls <- callBins(object)
cgh <- makeCgh(calls)
regions <- CGHregions(cgh)

Now, object regions contains a data set of reduced dimensions.)

[jiading]

I notice in the plot() from plot-methods.R, there is a function related to segmentation "*doSegments*". My guess is that QDNAseq did that analysis. Is it possible to play some tricks there to get the segmented file?

…

On Tue, Mar 7, 2017 at 4:26 PM, Ilari Scheinin ***@***.***> wrote: I know what you mean, but that's not implemented, at least not currently. It's something for the special case of single samples only, and not for data sets of two or more samples. Starting from the table you pasted above, you could use the rle() function and some basic R coding to achieve what you're after. (An approach for larger data sets that includes this kind of dimensionality reduction could be for example something along the lines of: library(CGHregions) calls <- callBins(object)cgh <- makeCgh(calls)regions <- CGHregions(cgh) Now, object regions contains a data set of reduced dimensions.) — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#20 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAc-rVkaLeohWeiz8-kG2JK9CqPA8Icuks5rjXcMgaJpZM4HHihe> .

[jiading]

I've also tried:

regions <- CGHregions(cgh)

Error in CGHregions(cgh) : object 'ncolm' not found Is there way to avoid this error? Thanks.

…

On Tue, Mar 7, 2017 at 4:26 PM, Ilari Scheinin ***@***.***> wrote: I know what you mean, but that's not implemented, at least not currently. It's something for the special case of single samples only, and not for data sets of two or more samples. Starting from the table you pasted above, you could use the rle() function and some basic R coding to achieve what you're after. (An approach for larger data sets that includes this kind of dimensionality reduction could be for example something along the lines of: library(CGHregions) calls <- callBins(object)cgh <- makeCgh(calls)regions <- CGHregions(cgh) Now, object regions contains a data set of reduced dimensions.) — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#20 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAc-rVkaLeohWeiz8-kG2JK9CqPA8Icuks5rjXcMgaJpZM4HHihe> .

[ilarischeinin]

Your object cgh is not of a type that CGHregions() is expecting. You must have missed one of the other lines of code I gave. This command class(cgh) should return "cghCall", but probably doesn't right now. If you cannot figure it out, include the outputs of print(cgh) and sessionInfo() here.

By the way, this discussion really is in the wrong place. It has very little to do with the QDNAseq package. There are e.g. mailing lists for general support on how to use R or Bioconductor.

[jiading]

My data is not the standard. All I tried to do is to get the segmentation from the *callBins*. I am afraid my next question is related to "CGHregions" pacakge :)

copyNumbersCalled <-

try(callBins(copyNumbersSegmented,method="cutoff",cutoffs=c(-0.1, 0.1)),silent=T) Calling aberrations with the following cutoffs: loss < -0.1 < normal < 0.1 gain

copyNumbersCalled

QDNAseqCopyNumbers (storageMode: lockedEnvironment) assayData: 30970 features, 1 samples element names: calls, copynumber, probgain, probloss, probnorm, segmented protocolData: none phenoData sampleNames: GC033602.freeze1.4 varLabels: name total.reads ... loess.family (6 total) varMetadata: labelDescription featureData featureNames: 1:1-100000 1:100001-200000 ... Y:59300001-59373566 (30970 total) fvarLabels: chromosome start ... use (9 total) fvarMetadata: labelDescription experimentData: use 'experimentData(object)' Annotation:

cgh <- makeCgh(copyNumbersCalled) cgh

cghCall (storageMode: lockedEnvironment) assayData: 26150 features, 1 samples element names: calls, copynumber, probgain, probloss, probnorm, segmented protocolData: none phenoData sampleNames: GC033602.freeze1.4 varLabels: name total.reads ... loess.family (6 total) varMetadata: labelDescription featureData featureNames: 1:900001-1000000 1:1000001-1100000 ... Y:59000001-59100000 (26150 total) fvarLabels: Chromosome Start ... use (9 total) fvarMetadata: labelDescription experimentData: use 'experimentData(object)' Annotation:

regions <- CGHregions(cgh)

*Error in datareg[nrow(datareg), ] : incorrect number of dimensions*

…

On Fri, Mar 10, 2017 at 1:26 PM, Ilari Scheinin ***@***.***> wrote: Your object cgh is not of a type that CGHregions() is expecting. You must have missed one of the other lines of code I gave. This command class(cgh) should return "cghCall", but probably doesn't right now. If you cannot figure it out, include the outputs of print(cgh) and sessionInfo() here. By the way, this discussion really is in the wrong place. It has very little to do with the QDNAseq package. There are e.g. mailing lists for general support on how to use R or Bioconductor. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#20 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAc-rQba5Ia8QCrFE2XHw9-xaJeOYkXYks5rkUF0gaJpZM4HHihe> .

[ilarischeinin]

I see that fixed the previous error message you saw. But when I said that using CGHregions was "an approach for larger data sets", I meant that; it does not work for individual samples (for which I suggested using rle() plus some manual R coding to achieve what you're after). If you want to read more about CGHregions, here is the paper.