Mapping using Tabix and VCF files

Mapping using Tabix and VCF files#

We have seen how Ensembl can be used to localise variants in the genome. This works reasonably well for low throughput applications, however, it is not scalable. If you want to validate thousands or millions of variants then you will need to use locally stored files. This example demonstrates how to use a tabix indexed VCF file to map variants. The mapping file used in this example is a very small test one used in the package.

Having said that VCF files are locally stored, this does not have to be the case as all the files are read with pysam, which is a Python interface to the htslib library that allows remove access to VCF files. However, I do not think that Ensembl or the NCBI would thank you for hitting their FTP servers with long running variant mapping tasks.

If you have not done so already, It is recommended that you work your way through the Ensembl mapping tutorial first, as this contains some more background than this one.

[11]:

import csv
import gzip
import pandas as pd
from variant_mapper import mapper, examples, constants as con

We will download the test mapper and some test variants to map against it.

[2]:

outdir = "/data/test-mapper"

[3]:

files = examples.download_test_mapper(outdir)

The files are listed below. There is a primary and secondary mapper, in this example we will just use one of them. When we come to the scan mapper example in the next workbook, we will use both. There is also a reference genome sequence small-mepper-genome.faa.gz. In addition, there are indexes for all of these. Finally, there are some test variants to run the mapping on.

[4]:

for i in files:
    print(i)

/data/test-mapper/small-mapper-genome.faa.gz
/data/test-mapper/small-mapper-genome.faa.gz.fai
/data/test-mapper/small-mapper-genome.faa.gz.gzi
/data/test-mapper/small-mapper-primary-mapper.vcf.gz
/data/test-mapper/small-mapper-primary-mapper.vcf.gz.tbi
/data/test-mapper/small-mapper-secondary-mapper.vcf.gz
/data/test-mapper/small-mapper-secondary-mapper.vcf.gz.tbi
/data/test-mapper/test-variants.txt.gz

[5]:

with gzip.open('/data/test-mapper/test-variants.txt.gz', 'rt') as infile:
    reader = csv.DictReader(infile, delimiter="\t")
    rows = [i for i in reader]
pd.DataFrame(rows)

[5]:
chr pos rsid ref alt notes
0 1 259 . T G all-match
1 1 3368 . C T ref-flip
2 1 5957 . G GGGTTGTCGC multi-site
3 1 6500 . A T missing
4 1 8366 . T A all-match
5 2 1613 . TCCTG T all-match
6 2 1807 . CTG C all-match
7 2 3136 . CCCATTC C all-match
8 2 8190 rs1941537900 T C all-match
9 2 8647 . G A missing
10 3 2984 rs1838945919 TTACGGACT T all-match
11 3 3621 . G T all-match
12 3 3646 . A G all-match
13 3 8244 . T G all-match
14 4 217 rs1146836854 A T all-match
15 4 3402 rs1684565759 C A all-match

Now we will use the primary mapper to run these, using Tabix as the method to localise the variant sites to the genome.

[12]:

mapper_file = "/data/test-mapper/small-mapper-primary-mapper.vcf.gz"

# Pass the mapping VCF file to the tabix mapper
with mapper.TabixVcfVariantMapper(mapper_file) as vmap:
    for i in rows:
        chr_name, start_pos, ref_allele, alt_allele = i['chr'], int(i['pos']), i['ref'], i['alt']

        # Map the variant
        mapped = vmap.map_variant(chr_name, start_pos, ref_allele, alt_allele=alt_allele)

        # Extract the mapped data
        i['map_bits'] = mapped.map_bits
        i['nsites'] = mapped.nsites
        try:
            # Get the mapping alleles
            i['mapped_ref'] = mapped.mapping_coords.ref_allele
            i['mapped_alt'] = mapped.mapping_coords.alt_allele
        except AttributeError:
            # When the variants can't be mapped, mapped.mapping_coords will be None
            i['mapped_ref'] = None
            i['mapped_alt'] = None
        # Decode the mapping bits into a human readable list
        i['decoded_map_bits'] = con.decode_mapping_flags(mapped.map_bits)

The results can be seen below, we have added the following columns:

  • map_bits - The mapping bits, the fields that have aligned with the site in the mapping file.

  • nsites - The number of other variant sites that overlap the mapping site.

  • mapped_ref - The reference allele that has been mapped.

  • mapped_alt - The alternate allele that has been mapped

  • decoded_map_bits - A human readable form of the mapping bits.

We can see that rows 3 and 9 have not mapped as all, as they are not present in the mapping file.

[13]:

pd.DataFrame(rows)

[13]:
chr pos rsid ref alt notes map_bits nsites mapped_ref mapped_alt decoded_map_bits
0 1 259 . T G all-match 28160 1 T G [CHR, START, STRAND, REF, ALT]
1 1 3368 . C T ref-flip 28224 1 T C [CHR, START, STRAND, REF, ALT, REF_FLIP]
2 1 5957 . G GGGTTGTCGC multi-site 28160 3 G GGGTTGTCGC [CHR, START, STRAND, REF, ALT]
3 1 6500 . A T missing 0 0 None None [NO_DATA]
4 1 8366 . T A all-match 28160 1 T A [CHR, START, STRAND, REF, ALT]
5 2 1613 . TCCTG T all-match 28160 2 TCCTG T [CHR, START, STRAND, REF, ALT]
6 2 1807 . CTG C all-match 28160 1 CTG C [CHR, START, STRAND, REF, ALT]
7 2 3136 . CCCATTC C all-match 28160 1 CCCATTC C [CHR, START, STRAND, REF, ALT]
8 2 8190 rs1941537900 T C all-match 28160 1 T C [CHR, START, STRAND, REF, ALT]
9 2 8647 . G A missing 0 0 None None [NO_DATA]
10 3 2984 rs1838945919 TTACGGACT T all-match 28160 1 TTACGGACT T [CHR, START, STRAND, REF, ALT]
11 3 3621 . G T all-match 28160 1 G T [CHR, START, STRAND, REF, ALT]
12 3 3646 . A G all-match 28160 2 A G [CHR, START, STRAND, REF, ALT]
13 3 8244 . T G all-match 28160 1 T G [CHR, START, STRAND, REF, ALT]
14 4 217 rs1146836854 A T all-match 28160 1 A T [CHR, START, STRAND, REF, ALT]
15 4 3402 rs1684565759 C A all-match 28160 2 C A [CHR, START, STRAND, REF, ALT]

Now we will perform the same operation but passing the reference genome assembly for the mapping file. This time, if a variant does not map then it is checked for a reference allele match in the reference genome assembly.

[21]:

# Pass the mapping VCF file to the tabix mapper
with mapper.TabixVcfVariantMapper(mapper_file, ref_genome="/data/test-mapper/small-mapper-genome.faa.gz") as vmap:
    for i in rows:
        chr_name, start_pos, ref_allele, alt_allele = i['chr'], int(i['pos']), i['ref'], i['alt']

        # Map the variant
        mapped = vmap.map_variant(chr_name, start_pos, ref_allele, alt_allele=alt_allele)

        # Extract the mapped data
        i['map_bits'] = mapped.map_bits
        i['nsites'] = mapped.nsites
        try:
            # Get the mapping alleles
            i['mapped_ref'] = mapped.mapping_coords.ref_allele
            i['mapped_alt'] = mapped.mapping_coords.alt_allele
        except AttributeError:
            # When the variants can't be mapped, mapped.mapping_coords will be None
            i['mapped_ref'] = None
            i['mapped_alt'] = None
        # Decode the mapping bits into a human readable list
        i['decoded_map_bits'] = con.decode_mapping_flags(mapped.map_bits)
        # print(mapped)

Here we can see that the mapping bits for the two rows that do not match are errors, as they do not even map to the reference assembly.

[15]:

pd.DataFrame(rows)

[15]:
chr pos rsid ref alt notes map_bits nsites mapped_ref mapped_alt decoded_map_bits
0 1 259 . T G all-match 28160 1 T G [CHR, START, STRAND, REF, ALT]
1 1 3368 . C T ref-flip 28224 1 T C [CHR, START, STRAND, REF, ALT, REF_FLIP]
2 1 5957 . G GGGTTGTCGC multi-site 28160 3 G GGGTTGTCGC [CHR, START, STRAND, REF, ALT]
3 1 6500 . A T missing 1 0 None None [ERROR]
4 1 8366 . T A all-match 28160 1 T A [CHR, START, STRAND, REF, ALT]
5 2 1613 . TCCTG T all-match 28160 2 TCCTG T [CHR, START, STRAND, REF, ALT]
6 2 1807 . CTG C all-match 28160 1 CTG C [CHR, START, STRAND, REF, ALT]
7 2 3136 . CCCATTC C all-match 28160 1 CCCATTC C [CHR, START, STRAND, REF, ALT]
8 2 8190 rs1941537900 T C all-match 28160 1 T C [CHR, START, STRAND, REF, ALT]
9 2 8647 . G A missing 1 0 None None [ERROR]
10 3 2984 rs1838945919 TTACGGACT T all-match 28160 1 TTACGGACT T [CHR, START, STRAND, REF, ALT]
11 3 3621 . G T all-match 28160 1 G T [CHR, START, STRAND, REF, ALT]
12 3 3646 . A G all-match 28160 2 A G [CHR, START, STRAND, REF, ALT]
13 3 8244 . T G all-match 28160 1 T G [CHR, START, STRAND, REF, ALT]
14 4 217 rs1146836854 A T all-match 28160 1 A T [CHR, START, STRAND, REF, ALT]
15 4 3402 rs1684565759 C A all-match 28160 2 C A [CHR, START, STRAND, REF, ALT]

Finally, lets try to map a variant that is not in the mappping file but is in the reference assembly to see what happens.

[19]:

# Pass the mapping VCF file to the tabix mapper
with mapper.TabixVcfVariantMapper(mapper_file, ref_genome="/data/test-mapper/small-mapper-genome.faa.gz") as vmap:
    # Map the variant
    mapped = vmap.map_variant("1", 259, "T", alt_allele="A")
    print(mapped)

MappingResult(source_coords=MapCoord(chr_name='1', start_pos=259, strand=1, ref_allele='T', alt_allele='A'), mapping_coords=MapCoord(chr_name='1', start_pos=259, strand=1, ref_allele='T', alt_allele='A'), map_bits=32768, source_row=None, map_row=None, errors=None, nsites=1, resolver=None)

We can see that while there is no mapping variant the mapping bits now indicate that the reference assembly does cantain the reference allele.

[20]:

con.decode_mapping_flags(mapped.map_bits)

[20]:

['REF_GENOME_MATCH']

Summary#

The above is the very basic usage of the tabix mapper. we have not covered using resolvers to extract data from the mapping file. will eventually update this notebook to give some examples.

Contents