Publication: https://doi.org/10.1101/2023.11.01.565041
This tool counts the number of specific k-mers within sequence data. The counts can then be compared to other counts to determine and compute the probability that samples are of the same origin to discover incongruent samples or sample swaps. It is intended to be run before any analysis and can provide additional QC information like sequencing error rate.
By default, this tool will only return pairs of samples with the same origin. This tool can theoretically also be used for relatedness inference using the -a parameter, however in this case the PCA-based heuristic should not be used (see below).
General:
- GCC (Tested on 9.3.0)
- zlibdev
- Autotools (if directly cloning from repo)
For generating site fasta files given a VCF file
- Python (Tested with 3.8.5)
- pyfaidx python module
- scikit-learn python module
- Perl (Tested with v5.26.2)
- bwa (Tested with 0.7.17)
ntsm can be installed with conda using the command:
conda install bioconda::ntsmIf cloning directly from the repository run:
./autogen.shCompiling should be as easy as:
./configure && makeTo install in a specified directory:
./configure --prefix=/PATH && make installGiven a VCF file and a reference genome, you can produce fasta files with k-mers that one can use to create fingerprinting. We have provided a set of human data based on similar criteria found in SNP microarrays.
The VCF file in this stage can be a single sample VCF, it just needs the variants.
Example:
ntsmSiteGen generate-sites name=prefix ref=reference.fa vcf=snps.vcfCreates site fasta files referred to as prefix_n{min missing sub k-mers}.fa which is the set of k-mers used by the tool. By default all non C/G <-> A/T conversions are ignored.
Parameters:
w=31 #window size to consider sequences in this region
k=19 #kmer size used in the window region
t=4 #threads for any subprocess or tools
The sites fasta file generated is a set of k-mers for each site. n for each fasta file generated refers to the number of missing k-mers allowed after filtering out repetitive k-mers (e.g. n=0 would mean you need each site to have all possible k-mers). The max value of n is w - k (so with default parameters 12). Small values of n will help you handle sequencing errors, but you need to have enough sites for our statistical analysis to work. In our testing retaining around ~10^5 sites seemed to work well, though if you end up with slightly less or a lot more ntsm should still function well.
If you do not wish to select your own sites, we currently include data/human_sites_n10.fa a fasta file with selected 96287 sites adequate for sample swap detection for human samples in the data folder.
Once fasta files for sites have been created, it is possible to create a PCA rotation matrix for speeding up the analysis. To do so you must supply a multiVCF file from which the PCA will be built. This multi-sample VCF ideally should not contain the same samples as the VCF used in the sample swap detection process. It should be a set of reliable samples on which a PCA and rotational matrix would be based on. We note that the use of a rotational matrix is optional.
Example (can be run in same directory as ntsmSiteGen generate-sites after picking a sites file to use):
ntsmSiteGen generate-pca-rot-mat name=prefix ref=reference.fa multivcf=/mnt/1886F90E86F8ED5EmultiVCF.vcf sites=prefix_n10.fa
Again if you are working with human samples and do not wish to generate your own, we currently include data/human_sites_rotationalMatrix.tsv and human_sites_center.txt to use in our PCA-based heurstic. We based our PCA and rotation matrix on 3202 samples from the 1000 Genomes Project.
Using this set of k-mers we can then count all of these k-mers within a fastq file. Files may be gzipped and multiple threads can be used. Each sample needs a separate run of this command and its own count files.
Example:
ntsmCount -t 2 -s sites.fa sample_part1.fq sample_part2.fq > counts.txtCreates count file using 2 threads. A sliding window using 19-mers is used in this case and the highest count in the window is recorded.
If your files are unsorted and have massive coverage, you may also intentionally run less reads using the -m parameter:
ntsmCount -t 2 -m 10 -s sites.fa sample_part1.fq sample_part2.fq > counts.txtThis will run the file until the average site coverage reaches 10x, which should be adequate for most sequencing data types. Lower coverage is possible if the read error rate is low enough.
Output Example:
#@TK 119443488624
#@KS 19
#locusID countAT countCG sumAT sumCG distinctAT distinctCG
rs1741692 23 22 68 190 3 9
rs6419870 0 43 0 86 1 2
rs3171927 19 20 91 20 5 1
rs12057128 16 0 31 0 2 5
rs11976368 43 0 43 0 1 10
rs4545798 17 13 34 65 2 5
rs10888802 0 37 0 355 9 10
...
Header lines (#@) helps in error rate estimation.
Example command:
ntsmEval HG002_rep1_counts.txt HG002_rep2_counts.txt HG003_counts.txt HG004_counts.txt > summary.tsvor if you wish optionally to speed up the analysis using a PCA rotation matrix:
ntsmEval -a -t 16 -n data/human_sites_center.txt -p data/human_sites_rotationalMatrix.tsv HG002_rep1_counts.txt HG002_rep2_counts.txt HG003_counts.txt HG004_counts.txt > summary.tsvOutput Example (with -a option):
sample1 sample2 score same dist relate ibs0 ibs2 homConcord het1 het2 sharedHet hom1 hom2 sharedHom n cov1 cov2 errorRate1 errorRate2 miss1 miss2 allHom1 allHom2 allHet1 allHet2
HG002_rep1_counts.txt HG002_rep2_counts.txt 0.07988 1 0.004839 0.996827 0 95971 0.998287 33720 33787 33613 62532 62465 62358 96252 37.416162 45.260554 0.003493 0.004301 35 35 62532 62465 33720 33787
HG003_counts.txt HG004_counts.txt 3.430842 0 7.512569 -0.003973 6649 48672 0.355549 33473 33781 13165 62772 62464 35507 96245 44.931787 44.068285 0.005968 0.004208 34 38 62779 62466 33474 33783
HG002_rep1_counts.txt HG003_counts.txt 1.660803 0 4.732675 0.498775 24 62518 0.731288 33720 33474 16744 62528 62774 45774 96248 37.416162 44.931787 0.003493 0.005968 35 34 62532 62779 33720 33474
HG002_rep1_counts.txt HG004_counts.txt 1.653478 0 2.872071 0.500089 19 62525 0.72982 33720 33783 16901 62525 62462 45624 96245 37.416162 44.068285 0.003493 0.004208 35 38 62532 62466 33720 33783
HG002_rep2_counts.txt HG003_counts.txt 1.760081 0 4.707002 0.499821 24 62521 0.73156 33787 33474 16779 62461 62774 45742 96248 45.260554 44.931787 0.004301 0.005968 35 34 62465 62779 33787 33474
HG002_rep2_counts.txt HG004_counts.txt 1.74858 0 2.78644 0.4996 19 62488 0.729034 33787 33783 16916 62458 62462 45572 96245 45.260554 44.068285 0.004301 0.004208 35 38 62465 62466 33787 33783
...
Note that the -a parameter will not output all pairwise comparisons if the PCA heuristic is used and will only consider comparing samples within a certain distance in PCA space. Though related samples are more likely to exist closer in PCA space, if using ntsm for relatedness inference you may miss some related pairs if you use this heuristic. You would not likely want to use -a with the PCA heuristic in most cases and the example above is for more illustrative purposes of the tool's use and output.
Column explanations:
- sampleX: Filename for sample X
- score: Log-likelihood-based score to determine if samples are the same or differ
- same: 1 means the tool thinks the sample is the same and 0 is if the tool thinks they differ
- dist: Distance of sample in PCA space
- relate: Relatedness determined via shared heterozygous sites
- relate: Relatedness determined via shared heterozygous sites
- ibs0: Number of sites with alleles not shared between two samples
- ibs2: Number of sites with the same genotype between two samples
- homConcord: Homozygous concordance determined via shared homozygous sites
- hetX: Number of heterozygous sites for sample X for all sites considered in comparison
- sharedHets: Number of shared heterozygous sites
- homX: Number of homozygous sites for sample X for all sites considered in comparison
- sharedHom: Number of shared homozygous sites
- n: number of unfiltered sites used in comparison
- covX: Coverage of sample X
- errorRateX: Error rate of sample X. May underestimate error caused by long indels.
- missX: Total number of missing sites in sample X
- allHomX: Total number of homozygous sites in sample X
- allHetX: Total number of heterozygous sites in sample X
If run on a single counts file the output can look like this:
sample cov errorRate miss hom het PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19 PC20
HG002_rep1_counts.txt 37.416162 0.003493 35 62532 33720 -14.254352 -23.285693 -0.373179 -7.315873 -1.187992 -5.494577 -0.434657 -0.334589 -0.832297 -1.160507 0.286102 -0.114464 1.013333 0.252766 -0.204102 0.465836 0.694361 0.099620 0.019345 -1.195279
This provides generic QC information useful without other samples (e.g. error rate) and if inclined a means of plotting the sample relative to others on a PCA plot. The number of columns is dependant on the number of principal components used.