First look at the gnomAD SV catalog

TL;DR

• The vast majority of SVs look good and there is a lot of information in the VCF.
• Some insertions and complex variants have contradicting POS/END fields (not sure what to do with those).
• Thousands of SVs are so similar that they could be the same SVs that weren’t merged.
• In a few cases, duplicates might lead to under-estimating the allele frequency.
• Take-home message: don’t annotate SVs using only the best matching SV in gnomAD.
• Most likely, the majority of variants are not really sequence-resolved (like in most SV catalogs).

SV types

Detection methods

How many variants are supported by each method, for each SV type:

Note: the same variant is counted several times in this table if it is supported by multiple methods.

Interesting to see that MELT supports a few duplications (are they MEI close to an existing TE?) and that read depth changes support a few insertions (maybe some very recent MEI?).

Allele frequency

• AF field: sometimes one value (when biallelic), sometimes multiple values.
• Only MCNVs contain multiple AF values (one for each copy number).
• All MCNVs have a CN=2 allele.

I want one value of allele frequency per variant so I sum up the AF values across alleles excluding the CN=2 allele for MCNVs.

Size distribution

The peak at 5 kbp is a bit puzzling. Maybe it’s a technical artifact, e.g. due to a method that detects DEL/DUP only down to 5 kbp. Looking at the variants with absolute SVLEN of exactly 5 kbp:

It’s most likely an artifact of the depth algorithm. Does it affect the allele frequency distribution?

The tail of the distribution is longer in the 5-6 kbp class.

SV coordinates confusion

It seems like the END value is not always larger than the position POS in the VCF. Most of it comes from CTX (reciprocal chromosomal translocation) variants or translocations which makes sense: the END might relate to the second breakpoint, potentially anywhere in the genome.

Some are insertions and I’m not sure how to interpret them. The vast majority of these variants have a PASS filter.

For most of these insertions the END is not that far from the POS (<100bp) but for some the END and POS are hundreds of kbp apart. Hum. Let’s have a look at some of those:

## GRanges object with 3 ranges and 13 metadata columns:
##                         seqnames                 ranges strand |
##                            <Rle>              <IRanges>  <Rle> |
##   gnomAD_v2_INS_4_34681        4 [116015282, 116015282]      * |
##   gnomAD_v2_INS_2_13594        2 [ 85885163,  85885163]      * |
##   gnomAD_v2_INS_9_69619        9 [ 28199346,  28199346]      * |
##                         paramRangeID            REF             ALT
##                             <factor> <DNAStringSet> <CharacterList>
##   gnomAD_v2_INS_4_34681         <NA>              N    <INS:ME:ALU>
##   gnomAD_v2_INS_2_13594         <NA>              N    <INS:ME:ALU>
##   gnomAD_v2_INS_9_69619         <NA>              N           <INS>
##                              QUAL      FILTER      SVTYPE     SVLEN
##                         <numeric> <character> <character> <integer>
##   gnomAD_v2_INS_4_34681       646        PASS         INS       281
##   gnomAD_v2_INS_2_13594       652        PASS         INS       279
##   gnomAD_v2_INS_9_69619       407        PASS         INS       136
##                               END        CHR2    CPX_TYPE          methods
##                         <integer> <character> <character>      <character>
##   gnomAD_v2_INS_4_34681  44565090           4        <NA> delly_manta_melt
##   gnomAD_v2_INS_2_13594   8372664           2        <NA> delly_manta_melt
##   gnomAD_v2_INS_9_69619  17898584           9        <NA>      delly_manta
##                                            ID        AF
##                                   <character> <numeric>
##   gnomAD_v2_INS_4_34681 gnomAD_v2_INS_4_34681  0.001490
##   gnomAD_v2_INS_2_13594 gnomAD_v2_INS_2_13594  0.005355
##   gnomAD_v2_INS_9_69619 gnomAD_v2_INS_9_69619  0.057464
##   -------
##   seqinfo: 24 sequences from an unspecified genome

I don’t see anything else suspicious.

For complex variants, the POS and END are always less than 100 bp apart and all for one type of complex SV: “dispersed duplications”. If I understand correctly, this type of complex SV is pretty much an insertion of a sequence already in the genome (but not a transposable element I guess?).

I’m not sure how to interpret them either.

For now, I’ll use the POS field for both insertions and dispersed duplications when manipulating the variants as genomic intervals. So I’ll assume that they are both simple insertions (which I think is fair enough).

Potential duplicates

Are there any pairs of variants with a reciprocal overlap higher than 90%? For insertions, I match them if they are located at less than 30 bp from each other and their size is 90% similar.

Note: Here I don’t consider BND or CTX variants, and I count each pair of variants only once.

There are thousands of deletions, insertions and duplications with high reciprocal overlap. Also many complex variants but I can imagine that it’s difficult to merge those. For the vast majority of pairs, both variants have a PASS filter. However, for most of the complex variants duplicates, one of the variant doesn’t have a PASS filter. Using FILTER==PASS would remove almost all the complex SV duplicates.

Let’s look at the size distribution of these potential duplicates.

Wow, very “peaky”. Let’s zoom in with no log-scale:

• The “Alu” peak around 300 bp is made of around 200 insertion pairs.
• The “L1” peak around 6 kbp is made of around 30 insertion pairs and almost as many duplication pairs.
• Again the technical peak at 5 kbp, here mostly pairs of deletions.
• The complex variants that are potentially duplicated form two size clusters.

GiaB comparison to investigate sequence resolution

The Genome in a Bottle consortium combined many types of techonologies (ultra deep short-read, linked-reads, long reads, optical mapping) to generate a high-quality SV catalog. Maybe not all variants are correct but it’s one of the best quality catalog we have. I want to quickly compare these SVs with the ones in the gnomAD catalog. If many variants look exactly the same in both datasets, it would be good evidence that they are both sequence-resolved.

I downloaded the NIST SVs Integration v0.6. This catalog also uses GRCh37. It includes short indels that I decided to filter out, keeping only variants of size 50 bp or more.

I’m curious to compare the exact breakpoint locations of SVs shared by both catalogs. Do they have exactly the same breakpoints/sequence? To make things easier let’s just look at deletions.

Many deletions in the gnomAD catalog are extremely similar to the GiaB dataset but not exactly the same. Some examples (randomly selected):

This quick comparison suggests that the variants are not really sequence-resolved. In gnomAD, the GiaB dataset, or both, the breakpoints might have a few errors. That will be important when genotyping these SVs or injecting them in genome graphs.

Of note, I’m not sure if the variants were left-aligned. In this case it wouldn’t make a big difference because we see that the deletions are discordant in term of size also. Still, ideally both VCF should be normalized first to better estimate how many deletions are exactly matched.

CHR2 and END

CHR2 is supposed to be the “Chromosome for END coordinate”. However I noticed that it might not be the case sometimes. For example, when CHR2 is different from #CHROM but the END is very close to the variant’s position. In those cases it might be that the END was set to the position in the #CHROM/POS side, like for insertions for example.

It happens a lot (the y-axis is log-scale) but mostly in variants without a FILTER=PASS. Just a few dozens variants with FILTER=PASS may have a problem with their CHR2 or END.

Some examples:

##   seqnames    start      end width strand paramRangeID REF   ALT QUAL
## 1        1  8620717  8620717     1      *         <NA>   N <INS>  999
## 2        1  9121444  9121444     1      *         <NA>   N <CPX>  331
## 3        1 23636591 23636591     1      *         <NA>   N <CPX>  330
## 4        1 30893626 30893626     1      *         <NA>   N <INS>  534
## 5        1 59880893 59880893     1      *         <NA>   N <INS>  351
## 6        1 66188203 66188203     1      *         <NA>   N <INS>  175
##   FILTER SVTYPE SVLEN      END CHR2 CPX_TYPE methods                   ID
## 1   PASS    INS   277  8620717   22     <NA>   manta  gnomAD_v2_INS_1_295
## 2   PASS    CPX   718  9121444   14     dDUP   manta   gnomAD_v2_CPX_1_20
## 3   PASS    CPX   270 23636591    2     dDUP   manta   gnomAD_v2_CPX_1_57
## 4   PASS    INS   172 30893626    X     <NA>   manta gnomAD_v2_INS_1_1114
## 5   PASS    INS   340 59880893    2     <NA>   manta gnomAD_v2_INS_1_2293
## 6   PASS    INS   344 66188203    7     <NA>   manta gnomAD_v2_INS_1_2601
##         AF
## 1 0.000047
## 2 0.225787
## 3 0.000047
## 4 0.001762
## 5 0.001261
## 6 0.000047