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?