Closed oalavijeh closed 4 years ago
Hi @oalavijeh , I can at least answer your last question. You need the BSgenome.Hsapiens.UCSC.hg19
bioconductor package to make this work, that's what the hg19
refers to. See https://bioconductor.org/packages/release/bioc/manuals/VariantAnnotation/man/VariantAnnotation.pdf for details.
or just use vcf <- VariantAnnotation::readVcf(vcf.file)
without specifying a genome.
Can SVA find those that overlap by a certain percentage?
Yes. This is what findBreakpointOverlaps()
does. Essentially, it does what findOverlaps()
does but it calculates overlap for the two breakends. In practice, this means that for simple events the terminology is different. Since the comparison is between the intervals of uncertainty around the start and end of the SVs, maxgap
means the maximum distance the start|end can be away from each other before they are considered a match, and minoverlap
is kind of useless and just kept for API consistency with findOverlaps()
.
Since SVA can handle arbitrary SVs, it represents them not as {type, chr, start, end}, but as {breakend1, breakend2}. That is, {(chr, lower_pos_bound, upper_pos_bound, orientation), (chr, lower_pos_bound, upper_pos_bound, orientation)}. lower_pos_bound and upper_pos_bound for a breakend can be different either due to the call itself being IMPRECISE
(e.g. RP based callers), or due to a precise variant having a microhomology at the breakpoint position (e.g. flank polyA on both sides of both breakends)
In summary:
maxgap
= max allowable distance between start positions and max allowable distance between end positionsminoverlap
= not useful for SVssizemargin
= multiplier specifying max allowable difference in size. A margin of 0.5 will allow a 100bp event to match something in the [50, 200]bp range. This is need because if you allow maxgap=50
bp of wiggle room around your start and end positions, you still don't want to match a `bp event with a 99bp events.restrictMarginToSizeMultiple
= further restict maxgap
based on event size. A maxgap of 50bp and a sizeMargin of 0 will allow 10bp events that are 40bp apart from each other to match. Setting will dyamically lower maxgap
for small events such that effective maxgap = min(maxgap, event size * restrictMarginToSizeMultiple)
. restrictMarginToSizeMultiple=0.5
will require a 50% overlap between events. For simple SVs, you can think of this as the SV equivalent of minoverlap
.Can SVA then merge these?
No. Merging identical events is trivial but there's too many choices to make when merging SVs that do not match exactly. For example, if you have a 50bp maxgap
you have have event A overlap B, B overlap C, but A not overlap C. There's no right answer about how to merge these events.
Example: A: DEL [100, 150] B: DEL [110, 160] C: DEL [120, 170] D: DEL [130, 180] E: DEL [140, 190] F: DEL [150, 200] G: DEL [160, 210]
With maxgap=10
, each event overlaps the next event, but by the time you get to G, there's no common deleted bases at all. There's no right answer about what to do in such scenario so I have intentionally excluded such functionality from SVA.
Thanks Daniel and Christian for the very helpful answers!
Dear Daniel, I was trying to use SVA and had a few questions (that I'm sure are very simple so apologies):
I have divided Manta calls for a large cohort in vcf format that have been merged and filtered for basic QC. They have then been subdivided by type of call so each vcf represent just one type of SV merge across the cohort
Can SVA find those that overlap by a certain percentage? Can SVA then merge these? I can see once can then plot these positions which will be very useful.
In the example script on your vignette when I run:
suppressPackageStartupMessages(require(StructuralVariantAnnotation)) suppressPackageStartupMessages(require(VariantAnnotation)) vcf.file <- system.file("extdata", "gridss.vcf", package = "StructuralVariantAnnotation") vcf <- VariantAnnotation::readVcf(vcf.file, "hg19") gr <- breakpointRanges(vcf)
at "vcf <- VariantAnnotation::readVcf(vcf.file, "hg19")" I get "file does not exist" message. Is the hg19 a file from the VariantAnnotation package or a different reference file we are meant to provide?
Many thanks for your help