Core genome too short

Question

I'm trying to create a snp distance matrix from M tuberculosis isolates to try and infer transmission networks. My plan was to use snippy to make a core genome alignment and then snp-dists to create the matrix but have hit a snag with creating a core genome alignment. When I try it with all 972 of my samples, snippy halts when trying to run snp-sites, saying Warning: No SNPs were detected so there is nothing to output. The issue is that with all the samples included, the core genome size is 0.

When looking at the alignment statistics provided by snippy, it looks like some isolates have very low numbers of aligned base pairs (??due to low coverage or contamination). I have progressively filtered out samples with low numbers of aligned base pairs and can then get snippy/snp-sites to work. As I increase my filter stringency, the resulting core genome is still short. Eg Filtering out those with aligned base pairs of <90% of the reference had a core genome of only 5747bp (from a reference of 4.4 million bp), while excluding 88 (9%) of my samples.

snippy's author recommends using it's output file core.txt to figure out which samples are the "bad" outliers. That file provides 1) the length of the reference, 2) the number of aligned base pairs, 3) the number of unaligned base pairs, 3) the number of variant sites, 4) the number of heterogenous sites, 5) the number of masked sites, and 6) the number of low coverage sites.

#>   ID                        LENGTH ALIGNED UNALIGNED VARIANT   HET MASKED LOWCOV
#>   <chr>                      <dbl>   <dbl>     <dbl>   <dbl> <dbl>  <dbl>  <dbl>
#> 1 R15795_CATCAAGT_S34_L006 4411532 4192432      6818     755   428 209178   2676
#> 2 R15842_GTCTGTCA_S49_L006 4411532 4187344     10304     784   441 209178   4265
#> 3 R15876_CGCTGATC_S36_L006 4411532 4176662     18992    1292   660 209178   6040
#> 4 R15951_ATTGGCTC_S7_L002  4411532 4170649     14732    1281   980 209178  15993
#> 5 R16019_TGGAACAA_S78_L001 4411532 4190733      8132     712   715 209178   2774
#> 6 R16046_GACTAGTA_S6_L002  4411532 4186069      4309    1257  1084 209178  10892

My question is: what heuristic would you use to filter out samples prior to creating a core genome?

And secondarily, what would a reasonable core genome size be for Mtb?

Answer 1

Targeting Illumina reads with >70X depth-of-coverage is a good place to start. Once a robust database is established lower coverage isolates can be used effectively too. More specific to individual SNPs... filtering on VCF QUAL value, Map Quality and/or AC value have shown to be useful.

Core genome will likely be around 1,000 SNPs, but it will depend on the lineage the reference belongs to and the aligned sample lineage. If both are from the same lineage you're likely to be less than 700 SNPs. When from different lineages SNP counts will obviously be higher in addition to the alignment error seen.

Pertaining to filtering individual SNPs...

vSNP is a tool developed for SNP analysis. It can be applied to different specie types but was originally developed for use with TB complex isolates. Code detail is here.

vSNP is a 2-step process. First, generating VCF files. Second, using those VCF files to output SNP tables and trees. The goal of the script is to achieve a high resolution genotyping analysis that can be validated with minimal effort. Some steps are included to help achieve this such as reference selection, automatic grouping and position filtering.

Pseudocode:

For set of VCF files created from the same reference...

def find_positions()
    AC=2
    QUAL > 300 

For each found position...

def get_snps()
    AC=2 (solid call)
    QUAL > 50
        THEN --> SNP
    AC=1 (mix call)
    QUAL > 150
        THEN --> AMBIGIOUS
    AC=2
    QUAL >= 50
       THEN --> REF
    QUAL < 50
       THEN --> N
    ALT == "None"
       THEN --> -

First find all possible positions that could be informative in all samples/VCF files: find_positions(), then define each position to be included in the alignment: get_snps().

Then it can be convenient to filter consistently poor alignment regions such as PPE/PGRS.