CNANorm: A normalization method for Copy Number Aberration in cancer samples.

[stefanoberri]

CNAnorm performs ratio, GC content correction and normalization of data obtained using low coverage (one read every 100-10,000 bp) high throughput sequencing. It performs a "discrete" normalization looking for the ploidy of the genome. It will also provide tumour content if at least two ploidy states can be found.

Now published on Bioinformatics

ABSTRACT:
MOTIVATION: Comparison of read depths from next generation sequencing between cancer and normal cells makes the estimation of copy number alteration (CNA) possible, even at very low coverage. However, estimating CNA from patients' tumour samples poses considerable challenges due to infiltration with normal cells and aneuploid cancer genomes. Here we provide a method that corrects contamination with normal cells and adjusts for genomes of different sizes so that the actual copy number of each region can be estimated.
RESULTS: The procedure consists of several steps. First, we identify the multi-modality of the distribution of smoothed ratios. Then we use the estimates of the mean (modes) to identify underlying ploidy and the contamination level, and finally we perform the correction. The results indicate that the method works properly to estimate genomic regions with gains and losses in a range of simulated data as well as in two datasets from lung cancer patients. It also proves a powerful tool when analysing publicly available data from two cell lines (HCC1143 and COLO829).
AVAILABILITY: An R package, called CNAnorm, is available at http://www.precancer.leeds.ac.uk/cnanorm or from Bioconductor.

Further information, original data and a Perl script to produce input files from sam/bam are available on the author's website

I have also created the Wiki entry

[dmacmillan]

Hi stefanoberri, I wonder if you would be able to assist me with a bug I am encountering. It says there was a concatenation error at line 302 in bam2windows.pl. Has this bug been reported yet? Can you suggest anything?

[stefanoberri]

Hi.
This is the first time this error is reported.
you'd need to help me a bit more
- what platform are you using (operative system, version of Perl...)
- Report the exact error as provided by Perl (it usually point to a certain line of code or so on)
- What input files are you using? What parameters?
- Can you obtain a correct output using the bam files and parameter of the web page? (http://www.precancer.leeds.ac.uk/sof...asets/cnanorm/)

you can also use private message or email ([email protected]) to get in touch with me.

Thanks

[dmacmillan]

Unfortunately the error message was buried in the command line after "vi"ing into the bam2windows.pl to see the line that was causing the issue. What I said earlier was fairly close to the output, but since I am using very large .bam files it will take at least another day to run the code on them again.

As you probably figured from above I am running linux, specifically centOS 2.16.0
Perl version v5.8.8 built for x86_64-linux-thread-multi

The exact command I used was as follows:
perl bam2windows.pl /projects/MYUSERNAME/DLBCL/CNV/RG065/tumour/A01419_9_lanes_dupsFlagged.bam /projects/MYUSERNAME/DLBCL/CNV/RG065/normal/A01440_8_lanes_dupsFlagged.bam > A01419.tab

Note that these are actually LINKS to the bam files and not the bam files themselves. I discovered this afterwards and I am currently testing them on the included .bam files that were supplied.

I'll let you know how it goes, thanks for the very prompt response!

-cheers

[stefanoberri]

Hi.

How big are they? CNAnorm was designed for low coverage, however I am now using it for high coverage (30X) without problems. I would suggest you do not use the default 30 reads per window with high coverage, not to start at least. try

--window 100000

to test

As a suggestion, add

--saveTest testFile.tmp --saveControl controlFile.tmp

arguments. It will save intermediate files, so it will be easier to figure out when the error occurs. Also, as I was suggesting earlier, try with the little bam files on the web site first.

[dmacmillan]

Good News, I successfully got a .tab file produced using the sample .bam files provided. I am now running bam2windows on two smaller .bam files and also an error log is being written to a file so if there are any errors I will be able to share them. I'll let you know how it goes, thank you!

[dmacmillan]

Ok so I ran it on two smaller .bam files and successfully created a .tab file. However, the error log is full of two errors, each of which repeat multiple times over and over, they are (respectively):

Use of uninitialized value in numeric lt (<) at bam2windows.pl line 498.
Use of uninitialized value in concatenation (.) or string at bam2windows.pl line 302.

Do you have any advice?

Was I supposed to get a .RData file as well?

[stefanoberri]

Ok so I ran it on two smaller .bam files and successfully created a .tab file. However, the error log is full of two errors

Good, we are making progress. If you provide 2 very small bam files that still produce the errors (I guess it is a warning) I'll try and fix it.

I suspect, however, this is unrelated to the first error you reported.

Was I supposed to get a .RData file as well?

No. the bam2window.pl script only produces a tab file.

Thanks for helping with the debugging. I'll contact people in IT here to see if they can suggest a way to send large files across which is "approved" by the university.

[dmacmillan]

I will see if I can release them, but it is doubtful.
This may sound like a dumb question but I am not sure where I go from here, how do I use the .tab file? I can't seem to find any documentation or usage on this anywhere.

[stefanoberri]

You now need to use Bioconductor, see the webpage and vignette for instruction. The difference is that, instead of loaing built-in data with

data(LS041)

you load the data frame with

myData <- read.delim("path/to/table.tab")
CN <- dataFrame2object(myData)

[dmacmillan]

Oh I see, slowly but surely its starting to make sense. When I do the dataFrame2object(mydata) I get this

Error in ratio[normNoZ] <- obj@Test[normNoZ]/obj@Norm[normNoZ] :
NAs are not allowed in subscripted assignments

Perhaps it is due to the errors in the .tab file. I'll send you a download link to the .bam files I am using if I can.

[dmacmillan]

I can't share the full .bam file, but I can give you the header and the first few reads. I'll email them to you directly.

[stefanoberri]

Hi dmacmillan,

I looked at the files you sent me to figure out why they didn't work. The problem had to do with the reference genome you use that is different to what I was using. I had implemented a way to provide a different genome annotation, but it wasn't very straight forward. I have now changed the code so that it reads from the bam file header directly.

It should now work with your bam files.
I tried with the small ones you sent and it gave no errors.

Still, keep in mind that, if you perform GC correction, you will need to download the GC content file for your reference genome and the chromosomes will have to have the same names.

I strongly recommend to perform GC correction.

I have now uploaded bam2windows.pl to googlecode so that it is easier to get the latest version of the script and people could more easily contribute if they wish.

I also update the web page about CNAnorm and uploaded the actual bam files used to produce the tab file then used in CNAnorm (LS041)

Thank you for reporting the bug.

Best of luck with your analysis.

[dmacmillan]

Hi stefanoberri,

You are by far the most helpful developer I have ever dealt with, I truly appreciate your effort in helping me! Thank you very much.

[dmacmillan]

Hi stefanoberri, I am having a different issue now.

Firstly, this comes up in stderr:

Counting reads over window 10161 bp wide...
  Calculating GC content...
Argument "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0..." isn't numeric in division (/) at bam2windows.pl line 256, <IN> line 576287280.
Use of uninitialized value in addition (+) at bam2windows.pl line 262, <IN> line 576287280.

I still get a .tab file in the end, though it is only 10mb. Seems like it should be larger especially considering that I am running it on a tumour and normal bam file each of which is ~11gigs. So I tried running R on it anyways...

However, when running R, I put my tab file into a CN object and follow the steps in the manual. But when I get to the plotPeaks this happens

plotPeaks(CN, special1 = 'chrX', special2 = 'chrY')
     chrX does not have enough values to estimate density. Skipping.
     chrY does not have enough values to estimate density. Skipping.

Do you have any suggestions? I made sure the gc_reference I was using had the same names for the chromosomes as my data this time as well.