SEQanswers

Go Back   SEQanswers > Literature Watch



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

Reply
 
Thread Tools
Old 11-18-2011, 02:10 AM   #1
stefanoberri
Member
 
Location: Cambridge area, UK

Join Date: Jan 2010
Posts: 35
Default 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.
stefanoberri is offline   Reply With Quote
Old 01-06-2012, 08:49 AM   #2
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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?
dmacmillan is offline   Reply With Quote
Old 01-06-2012, 08:58 AM   #3
stefanoberri
Member
 
Location: Cambridge area, UK

Join Date: Jan 2010
Posts: 35
Default

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
stefanoberri is offline   Reply With Quote
Old 01-06-2012, 09:08 AM   #4
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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
dmacmillan is offline   Reply With Quote
Old 01-06-2012, 09:20 AM   #5
stefanoberri
Member
 
Location: Cambridge area, UK

Join Date: Jan 2010
Posts: 35
Default

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.
stefanoberri is offline   Reply With Quote
Old 01-06-2012, 10:44 AM   #6
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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 is offline   Reply With Quote
Old 01-06-2012, 12:03 PM   #7
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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.
dmacmillan is offline   Reply With Quote
Old 01-09-2012, 02:20 AM   #8
stefanoberri
Member
 
Location: Cambridge area, UK

Join Date: Jan 2010
Posts: 35
Default

Quote:
Originally Posted by dmacmillan View Post
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 View Post
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.
stefanoberri is offline   Reply With Quote
Old 01-09-2012, 08:16 AM   #9
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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.
dmacmillan is offline   Reply With Quote
Old 01-09-2012, 08:22 AM   #10
stefanoberri
Member
 
Location: Cambridge area, UK

Join Date: Jan 2010
Posts: 35
Default

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)
stefanoberri is offline   Reply With Quote
Old 01-09-2012, 08:41 AM   #11
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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 is offline   Reply With Quote
Old 01-09-2012, 10:31 AM   #12
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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.
dmacmillan is offline   Reply With Quote
Old 01-13-2012, 07:19 AM   #13
stefanoberri
Member
 
Location: Cambridge area, UK

Join Date: Jan 2010
Posts: 35
Default

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
stefanoberri is offline   Reply With Quote
Old 01-13-2012, 09:34 AM   #14
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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 is offline   Reply With Quote
Old 01-17-2012, 12:19 PM   #15
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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.
dmacmillan is offline   Reply With Quote
Old 01-18-2012, 07:35 AM   #16
stefanoberri
Member
 
Location: Cambridge area, UK

Join Date: Jan 2010
Posts: 35
Default

Quote:
Originally Posted by dmacmillan View Post
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 View Post
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 View Post
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

Last edited by stefanoberri; 01-19-2012 at 12:53 AM.
stefanoberri is offline   Reply With Quote
Old 01-25-2012, 08:38 AM   #17
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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.
dmacmillan is offline   Reply With Quote
Old 01-25-2012, 09:40 AM   #18
stefanoberri
Member
 
Location: Cambridge area, UK

Join Date: Jan 2010
Posts: 35
Default

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...
stefanoberri is offline   Reply With Quote
Old 02-02-2012, 08:51 AM   #19
dmacmillan
Member
 
Location: British Columbia

Join Date: Jan 2012
Posts: 49
Default

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?
dmacmillan is offline   Reply With Quote
Old 02-02-2012, 09:00 AM   #20
stefanoberri
Member
 
Location: Cambridge area, UK

Join Date: Jan 2010
Posts: 35
Default

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.
stefanoberri is offline   Reply With Quote
Reply

Tags
cancer genomics, ngs, normalization, software

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 03:24 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO