breseq-ext-cnv usage

[gabypetrungaro]

Hi @p3rao ,
I am trying to use breseq-ext-cnv to analyse my coverage tables created with BAM2COV. Basically, I am interested in using the software for detecting amplifications and deletions that do not appear in my Breseq output tables.

I have managed to install the software and run the help command. However, I don't quite understand how to use it, or what the different optional parameters are used for:

1. Could you maybe shortly explain what the python software does? I think with a brief explanation I could interpret the optional parameters better to understand what values of those (e.g. Window and Step size) I should choose for my samples.
2. I failed my first tries to run the software and I don't understand why... I added values for -ori and -ter when I first got this error, but doesn't seem to have changed anything:

$breseq-ext-cnv -i ./input/CP009273\:1-4631469.tab -ori 0 -ter 4631469

Ori has not been set (default value is None)
Ter has not been set (default value is None)
CP009273:1-4631469.tab: Calculating coverage and GC% across sliding window over the genome.
Traceback (most recent call last):
  File "/home/gaby/spyder-env/bin/breseq-ext-cnv", line 8, in <module>
    sys.exit(main())
  File "/home/gaby/spyder-env/lib/python3.10/site-packages/breseq_ext_cnv/get_CNV.py", line 811, in main
    cnv,smpl, ori_win, ter_win = run_HMM(
  File "/home/gaby/spyder-env/lib/python3.10/site-packages/breseq_ext_cnv/get_CNV.py", line 459, in run_HMM
    df, ori_loc, ter_loc = bias_correction(filepath, output, ori, ter, enforce, win, step, frag)
  File "/home/gaby/spyder-env/lib/python3.10/site-packages/breseq_ext_cnv/get_CNV.py", line 544, in bias_correction
    df = preprocess(filepath, win, step, frag)
  File "/home/gaby/spyder-env/lib/python3.10/site-packages/breseq_ext_cnv/get_CNV.py", line 21, in preprocess
    if (step > win) :
TypeError: '>' not supported between instances of 'NoneType' and 'NoneType'

Thank you in advance for your help,
Kind regards,
Gabriela (Senior Postdoc at University of Cologne)

[p3rao]

Hi Gabriela,

Sorry about this issue, I think the run is failing because the current version of this tool is not updated on pypi. Could you try re-installing it using pip install git+https://github.com/barricklab/breseq-ext-cnv.git and rerunning that command.
Regarding the explanation of this pipeline, I will work on a better readme file soon. But here's a short version:

1. The code takes your BAM2COV output (coverage data at all base positions in the genome) and groups the base positions into --window that slide across the genomic positions in steps of length --step-size. So this summarizes sequencing coverage in the local segments of the genome. The step size determines if the subsequent windows will overlap each other and by how many bases.
2. The script then corrects for the GC bias across these genomic segments based on the GC content of the sequence within these segments. The -f --frag_size parameter can be used to specify the fragment sizes during sequencing or PCR steps and thus were subject to GC bias.
3. Then it looks for and corrects coverage bias due to proximity to origin or terminus of DNA replication in circular chromosomes. If the origin or terminus is not known the script will try to detect the bias pattern and infer the approximate genomic segment corresponding to coverage bias related to ori/ter.

All these parameters should be optional where the default values (-w 200 -s 100 -f 500) are used or the values are inferred. So maybe for a first pass you could re-run the code with just the input file with default parameters and then adjust once you have the script running.
This is very much a work in progress so hopefully it will continue to improve.

Best,

-Pranesh

[gabypetrungaro]

Dear Pranesh,
thanks a lot for your answer and explanation. I have now managed to run your code and I get more or less the idea of what it does, I think :) However, I'm afraid my understanding of it is not yet at a point that I can make it work to analyse my data.

1. I would very much appreciate a bit more help with interpreting the csv output.
   I was inspecting mostly the table in the file CP009273:1-4631469_break_pts.csv , and what I can understand is that this table gives me the starting position of the windows where the 'state' changes, right? Like in the following table, I can interpret that I have a segment of size 901 that has state 1 and then another segment with state 2 that starts at 801. However I can see that that is not quite the case right? Because otherwise the second segment should start at Startpos = 901 and this is not what the table says.
   Startpos | State | Segment_Size
   0 | 1 | 901
   801 | 2 | 279000
   279701 | 1 | 1300
   280901 | 2 | 2300
   283101 | 1 | 4700
   287701 | 2 | 91000
   378601 | 1 | 4111500
   4490001 | 2 | 141400
   I get that this is related to the window size and steps, but:
   - could you maybe explain how exactly I should interpret these positions and sizes?
   - Also, what does the value of 'state' tells me? Does 2 mean that the copy number variation is 2 there?
   - Is the code also supposed to find deletions? Like is a copy number equal or close to zero a possibility?

2. I have run the code now for one of my samples where I know that there is a large deletion at 2635210-2636199 (length 989) and an amplification at 3394142-3394240 (length 98). I have tried several window sizes and step values now and I cannot get the code to detect these...
   When I plot lines at the break points in the table on top of the coverage plot, I don't really understand what the code is identifying if I am honest... I attach my plots and the output, in case you have some time to look at them and maybe explain why my output doesn't find the amplification and do find these cuts in weird/wrong regions.. (The last plot has the region 279701 to 287701 approx. and lines in the positions mentioned in the output table lines 3 to 6).

CP009273:1-4631469_CNV.csv
CP009273:1-4631469_break_pts.csv

Thank you in advance for taking the time to read this. I really appreciate your help and I hope at least some of my comments also help you to improve your code and documentation :)

Kind regards,
Gabriela.

[p3rao]

Hi Gabriela,

Apologies again for the delayed response. This project is currently sitting on the backburner due to other priorities. I am trying to work on improving documentation and how to interpret the output files. But to answer your question:

1. The "_break_pts.csv" file outputs the genomic segments that the code finds having unique copy numbers. This pipeline uses a Hidden Markov Model to predict the next "state" or copy number given the change from current coverage. The "Startpos" is the starting coordinate of that current state/copy-number and the length denotes the size of the segment that is part of the said state. The "_CNV.csv" file outputs the copy number prediction as well as the intermediate bias and normalization parameters per sliding "window" coordinate.
2. The coverage plots you attached look like those are real deletions and amplifications. I understand you tried multiple window and step sizes, but looks like in the output file you sent the window size you picked is 200, however the amplification length is 1/2 of the window size. I suggest trying smaller window sizes and step sizes that are 50% of the window size. This ensures region of interest is represented in a window (and not averaged out due to large window size). As a rule set the window size <= the smallest CNV you'd want to detect. Very small window sizes (<10) may cause other issues and predict false positives, so I'd be mindful of that.

Could you re-run this with a window size of 100 and a step size of 50? It should pick both your CNVs with this input setting. I would also recommend re-installing this from github. I added a new feature where you can define a --region (xxxx-yyyy) of interest and the copy number prediction plot will output only your region of interest.