VarScan-Somatic Copy Number Alteration (CNA) Calling

Question

[2nelly]

Hi all, reading carefully the documentation of VarScan CNA pipeline I noticed in step 4 the following suggestion:

If all of the data and segments are consistently above or below the neutral value (0.0), you can re-center the data points with VarScan copyCaller.

My data seem to belong in this category after plotting in R using DNAcopy package. All are consistently below 0.0. So, my question is how should I know how much I do need to re-center my data? I mean I can calculate this by eye, but how accurate can be that? Is there any proper way to calculate that?

Thank you in advance.

Answer 1

[Chris Miller]

Here's an old script I have laying around that does the recentering. Run it on the varscan copyCaller output, then repeat segmentation on the new recentered file.

[2nelly]

Many thanks Chris!!! Just a question to be sure. The output of this script is the recenter up or down value, right? Nothing more.

[mittu1602]

If you read the script it has recentre up and down both.

[mittu1602]

I am using your script for recentre and I get the following error,

no LOH regions to consider. Falling back to use the mean genome-wide CN value 
(may be dangerous in polyploid or heavily CN-altered tumors)
;recenter_baseline: 1.350223039787
submitting job: java -cp ~dkoboldt/Software/VarScan net.sf.varscan.VarScan copyCaller Output2 --output-file Output2.recentered --recenter-up 1.350223039787
sh: 1: bsub: not found

What possible reason could it be, I do not understand. Thank you

[Chris Miller]

Those are two separate errors. The first tells you that you haven't looked for regions of loss of heterozygosity, which is an optional step. The second tells you that it's trying to submit a job to a cluster via LSF. If you don't have one of those, then you'll have to alter the script accordingly (instead of submitting the command, just run it)

[mittu1602]

I tried to run the command, but failed to understand on how do you decide the "recenter_baseline value" it will be a great help.

[Chris Miller]

If you trace through the code you can see that it comes from the get_genome_mean_cn() function. So, it's operating under the assumption that the mean CN of the genome is 2, or at least, not very far off from 2. (median would probably be better, but that's besides the point).

[mittu1602]

Thank you for the explanation. :)