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

stefanoberri · 11-18-2011, 02:10 AM

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 · 01-06-2012, 08:49 AM

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 · 01-06-2012, 08:58 AM

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 (s.berri@leeds.ac.uk) to get in touch with me.

Thanks

dmacmillan · 01-06-2012, 09:08 AM

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 · 01-06-2012, 09:20 AM

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

Code:

--window 100000

to test

As a suggestion, add

Code:

--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 · 01-06-2012, 10:44 AM

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 · 01-06-2012, 12:03 PM

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 · 01-09-2012, 02:20 AM

Quote:

Originally Posted by 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

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.

Quote:

Originally Posted by dmacmillan

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 · 01-09-2012, 08:16 AM

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 · 01-09-2012, 08:22 AM

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

Code:

data(LS041)

you load the data frame with

Code:

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

dmacmillan · 01-09-2012, 08:41 AM

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 · 01-09-2012, 10:31 AM

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 · 01-13-2012, 07:19 AM

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 · 01-13-2012, 09:34 AM

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 · 01-17-2012, 12:19 PM

Hi stefanoberri, I am having a different issue now.

Firstly, this comes up in stderr:

Code:

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

Code:

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.

stefanoberri · 01-18-2012, 07:35 AM

Quote:

Originally Posted by dmacmillan

Code:

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.

It looks like you have something strange on line 576287280 of your GC file (maybe last line?) probably "\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..."

Quote:

Originally Posted by dmacmillan

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.

The size of the final tab file is largely independent on the size of the input sam/bam, but on the size of the window. If you set '--window 100000' it will produce approx 30000 lines, no matter how large the original bam files. It will provide, however, higher reads per window.

If you set '--readNum 50' instead, it will produce the number of windows required to have, on average, 50 reads per window in the sample with least number of reads. In this case, larger sam/bam files will produce larger tab files

WARNING! CNAnorm was designed for low coverage data. In particular, if you have large bam files and set parameters to have MANY windows, function peakPloidy and addDNACopy will take very long.
There are ways to solve this, I might write a tutorial in future...

Quote:

Originally Posted by dmacmillan

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

Code:

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.

This is just a warning, and it is beacause your reference chromosomes are not called 'chrX' and 'chrY' but, I guess, 'X' and 'Y', similarly, you should change

Code:

toSkip <- c("chrY", "chrM")

to

Code:

toSkip <- c("Y", "M")

I hope this helps

dmacmillan · 01-25-2012, 08:38 AM

Hey Stefano,

Do you have time to explain your ways to work with large .BAM files? I am managing to run your program on my large files without any errors. However, when following the manual I end up with ~20% tumour content estimation whereas the oncologist estimated our tumour content to be 96%. I feel I may be doing something wrong in the post processing of the tab file in R. Any suggestions would be appreciated.

stefanoberri · 01-25-2012, 09:40 AM

Hi.
To start, I would suggest to skip the smoothing step.

Code:

CN <- addSmooth(CN, lambda = 7)

If you have very high coverage, the signal should be already well smoothed, and further decresing noise might over-fit mixture models.

Lower estimation of tumour content can be due to several reasons such as high tumur heterogeneity (see the example of COLO829 as described in the sup material of the publication) or full genome duplication without homozygous deletions. In both these cases, the peak detection might be inaccurate or too "conservative" (most of the genome still diploid)

If you trust the oncologist estimation, you could use that piece of information to obtain a more reliable normalisation. If it is really almost pure, then the main peaks should get a ratio of 0, 0.5, 1, 1.5...

dmacmillan · 02-02-2012, 08:51 AM

Hi Stefano, I had another question. If I wanted to view the ploidy at each region (like in the tab file) how would I go about doing this?

stefanoberri · 02-02-2012, 09:00 AM

I am not sure I understand what you mean. How about

Code:

exportTable(CN, file = "CNAnorm_table.tab", show = 'ploidy')

It will produce a tab file with chromosome, position, ratio... and segmented normalised copy number (last column).

alternatively, you can do

Code:

ploidy <- 2 * segMean.n(CN)

However, it is the second time this question get asked. I'll put in the todo list an easier and documented way to get the ploidy from CN.

Also, I found a bug in exportTable. Please download the latest version (1.1.4) from the developers page.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
how to check copy number from certain gene	cliff	Bioinformatics	2	07-14-2014 03:52 PM
exome copy number and LOH	rcorbett	Bioinformatics	20	11-26-2013 01:22 AM
Copy number normalization in capture experiments?	anjulka	General	2	01-23-2012 01:07 PM
tools to detect copy number variants	nans_bn	Bioinformatics	2	01-12-2012 05:25 AM
Copy number analysis using 454 data	ps376	Bioinformatics	3	10-08-2011 05:56 AM

11-18-2011, 02:10 AM	#1
stefanoberri Member Location: Cambridge area, UK Join Date: Jan 2010 Posts: 35	CNANorm: A normalization method for Copy Number Aberration in cancer samples. 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 Last edited by stefanoberri; 11-18-2011 at 03:11 AM.

01-06-2012, 09:20 AM	#5
stefanoberri Member Location: Cambridge area, UK Join Date: Jan 2010 Posts: 35	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 Code: --window 100000 to test As a suggestion, add Code: --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.

01-06-2012, 12:03 PM	#7
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	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? Last edited by dmacmillan; 01-06-2012 at 12:26 PM.

01-09-2012, 08:16 AM	#9
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	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. Last edited by dmacmillan; 01-09-2012 at 08:20 AM.

01-06-2012, 08:49 AM	#2
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	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?

01-06-2012, 08:58 AM	#3
stefanoberri Member Location: Cambridge area, UK Join Date: Jan 2010 Posts: 35	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 (s.berri@leeds.ac.uk) to get in touch with me. Thanks

01-06-2012, 09:08 AM	#4
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	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

01-06-2012, 10:44 AM	#6
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	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!

01-09-2012, 08:22 AM	#10
stefanoberri Member Location: Cambridge area, UK Join Date: Jan 2010 Posts: 35	You now need to use Bioconductor, see the webpage and vignette for instruction. The difference is that, instead of loaing built-in data with Code: data(LS041) you load the data frame with Code: myData <- read.delim("path/to/table.tab") CN <- dataFrame2object(myData)

01-09-2012, 08:41 AM	#11
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	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.

01-09-2012, 10:31 AM	#12
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	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.

01-13-2012, 07:19 AM	#13
stefanoberri Member Location: Cambridge area, UK Join Date: Jan 2010 Posts: 35	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. Last edited by stefanoberri; 01-13-2012 at 07:22 AM. Reason: changing formatting

01-13-2012, 09:34 AM	#14
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	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.

01-17-2012, 12:19 PM	#15
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	Hi stefanoberri, I am having a different issue now. Firstly, this comes up in stderr: Code: 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 Code: 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.

01-25-2012, 08:38 AM	#17
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	Hey Stefano, Do you have time to explain your ways to work with large .BAM files? I am managing to run your program on my large files without any errors. However, when following the manual I end up with ~20% tumour content estimation whereas the oncologist estimated our tumour content to be 96%. I feel I may be doing something wrong in the post processing of the tab file in R. Any suggestions would be appreciated.

01-25-2012, 09:40 AM	#18
stefanoberri Member Location: Cambridge area, UK Join Date: Jan 2010 Posts: 35	Hi. To start, I would suggest to skip the smoothing step. Code: CN <- addSmooth(CN, lambda = 7) If you have very high coverage, the signal should be already well smoothed, and further decresing noise might over-fit mixture models. Lower estimation of tumour content can be due to several reasons such as high tumur heterogeneity (see the example of COLO829 as described in the sup material of the publication) or full genome duplication without homozygous deletions. In both these cases, the peak detection might be inaccurate or too "conservative" (most of the genome still diploid) If you trust the oncologist estimation, you could use that piece of information to obtain a more reliable normalisation. If it is really almost pure, then the main peaks should get a ratio of 0, 0.5, 1, 1.5...

02-02-2012, 08:51 AM	#19
dmacmillan Member Location: British Columbia Join Date: Jan 2012 Posts: 49	Hi Stefano, I had another question. If I wanted to view the ploidy at each region (like in the tab file) how would I go about doing this?

02-02-2012, 09:00 AM	#20
stefanoberri Member Location: Cambridge area, UK Join Date: Jan 2010 Posts: 35	I am not sure I understand what you mean. How about Code: exportTable(CN, file = "CNAnorm_table.tab", show = 'ploidy') It will produce a tab file with chromosome, position, ratio... and segmented normalised copy number (last column). alternatively, you can do Code: ploidy <- 2 * segMean.n(CN) However, it is the second time this question get asked. I'll put in the todo list an easier and documented way to get the ploidy from CN. Also, I found a bug in exportTable. Please download the latest version (1.1.4) from the developers page.