A tool to filter sites in a FASTA-format whole-genome pseudo-alignment

Overview

Core-SNP-filter

License GPL v3

This is a tool to filter sites (i.e. columns) in a FASTA-format whole-genome pseudo-alignment based on:

  • Whether the site contains variation or not.
  • How conserved the site is, i.e. contains an unambiguous base in a sufficient fraction of the sequences.

I wrote Core-SNP-filter because I was using Snippy, and the snippy-core command produces a core.full.aln file (contains all sites regardless of variation and conservation) and core.aln (only contains invariant sites with 100% conservation). I wanted a tool that could produce a core SNP alignment, but with more flexibility, e.g. including sites with ≥95% conservation.

Core-SNP-filter is efficient. On a small input alignment (2 Mbp in length, 100 sequences), it runs in seconds. On a large input alignment (5 Mbp in length, 5000 sequences, 25 GB file size), it takes less than 10 minutes and only uses ~50 MB of RAM.

Important caveat: Core-SNP-filter is only appropriate for DNA alignments, not protein alignments.

Usage

The executable named coresnpfilter takes a FASTA file as input. This must be an aligned FASTA file, i.e. all sequences must be the same length. The characters in the FASTA sequences can be bases (e.g. A or c), gaps (-) or any other ASCII character (e.g. N for ambiguous bases or X for masked bases). The input FASTA can be gzipped, and line breaks (multiple lines per sequence) are okay.

There are two main options:

  • -e/--exclude_invariant: if used, all invariant sites in the alignment are removed. A site counts as invariant if the number of unique unambiguous bases (A, C, G or T) at that site is one or zero. For example, a site with only A is invariant, but a site with both A and C is not invariant. Gaps and other characters do not count, e.g. a site with only A, N and - is invariant. Case does not matter, e.g. a site with only A and a is invariant.
  • -c/--core: at least this fraction of the sequences must contain an unambiguous base (A, C, G or T) at a site for the site to be included. The default is 0.0, i.e. sites are not filtered based on core fraction. If 1.0 is given, all sites with gaps or other characters will be removed, leaving an alignment containing only unambiguous bases. A more relaxed value of 0.95 will ensure that each site contains mostly unambiguous bases, but up to 5% of the bases can be gaps or other characters.

Core-SNP-filter outputs a FASTA alignment to stdout. The output will have the same number of sequences as the input, but (depending on the options used) the length of the sequences will likely be shorter. The header lines (names and descriptions) of the output will be the same as the input, and there will be no line breaks in the sequences (each sequence gets one line). Some basic information (input file, input sequence length, number of sequences and output sequence length) is printed to stderr.

Some example commands:

# Exclude invariant sites:
coresnpfilter -e core.full.aln > filtered.aln

# With a strict core threshold (same as Snippy's core.aln):
coresnpfilter -e -c 1.0 core.full.aln > filtered.aln

# With a slightly more relaxed core threshold:
coresnpfilter -e -c 0.95 core.full.aln > filtered.aln

# Use gzipped files to save disk space:
coresnpfilter -e -c 0.95 core.full.aln.gz | gzip > filtered.aln.gz

# Running without any options will work, but the output will be the same as the input:
coresnpfilter core.full.aln > filtered.aln

Full help text:

Core SNP filter

Usage: coresnpfilter [OPTIONS] 

Arguments:
    Input alignment

Options:
  -c, --core         Restrict to core genome (0.0 to 1.0, default = 0.0) [default: 0.0]
  -e, --exclude_invariant  Exclude invariant sites
      --verbose            Verbose output
  -h, --help               Print help
  -V, --version            Print version

Installation from pre-built binaries

Core-SNP-filter compiles to a single executable binary (coresnpfilter), which makes installation easy!

You can find pre-built binaries for common OSs/CPUs on the releases page. If you use one of these OSs/CPUs, download the appropriate binary for your system and put the coresnpfilter file in a directory that's in your PATH variable, e.g. /usr/local/bin/ or ~/.local/bin/.

Alternatively, you don't need to install Core-SNP-filter at all. Instead, just run it from wherever the coresnpfilter executable happens to be, like this: /some/path/to/coresnpfilter --help.

Installation from source

If you are using incompatible hardware or a different OS, then you'll have to build Core-SNP-filter from source. Install Rust if you don't already have it. Then clone and build Core-SNP-filter like this:

git clone https://github.com/rrwick/Core-SNP-filter.git
cd Core-SNP-filter
cargo build --release

You'll find the freshly built executable in target/release/coresnpfilter, which you can then move to an appropriate location that's in your PATH variable.

Demo dataset

This repo's demo.fasta.gz file is a pseudo-alignment made from 40 Klebsiella samples (the original was ~5 Mbp long but I subsetted it down to 100 kbp to save space). It has many gaps, invariant sites and Ns, so Core-SNP-filter can make it a lot smaller.

For example, you can use Core-SNP-filter to create an invariant-free 95%-core alignment:

coresnpfilter -e -c 0.95 demo.fasta.gz > demo_core.fasta

The stderr output will look like this:

Core-SNP-filter
  input file:                    demo.fasta.gz
  number of sequences:           40
  input sequence length:         10000
  invariant-A sites removed:     1394
  invariant-C sites removed:     1763
  invariant-G sites removed:     1849
  invariant-T sites removed:     1378
  other invariant sites removed: 322
  non-core sites removed:        2143
  output sequence length:        1151

You can then build a tree with a program such as IQ-TREE:

iqtree2 -s demo_core.fasta -T 4

Verbose output

Using the --verbose option will make Core-SNP-filter print a per-site table to stderr. Save it to file like this:

coresnpfilter -e -c 0.95 --verbose core.full.aln 1> filtered.aln 2> core_snp_table.tsv

This is mainly for debugging purposes, so you probably don't want to use it. But if you do, the columns are:

  1. pos: 1-based index of the input alignment site
  2. a: whether any sequence at this site contains A or a
  3. c: whether any sequence at this site contains C or c
  4. g: whether any sequence at this site contains G or g
  5. t: whether any sequence at this site contains T or t
  6. count: the number of sequences at this site which contain an unambiguous base
  7. frac: the fraction of sequences at this site which contain an unambiguous base
  8. var: whether there is any variation at this site (i.e. two or more of the a/c/g/t columns are true)
  9. keep: whether the site passed the filter and is included in the output

Boolean columns use 0 for no and 1 for yes.

For example, you can use this table to see which sites in your input alignment have made it into the output alignment:

awk '{if ($9==1) print $1;}' core_snp_table.tsv

License

GNU General Public License, version 3

You might also like...
Filter, Sort & Delete Duplicate Files Recursively
Filter, Sort & Delete Duplicate Files Recursively

Deduplicator Find, Sort, Filter & Delete duplicate files Usage Usage: deduplicator [OPTIONS] [scan_dir_path] Arguments: [scan_dir_path] Run Dedupl

A fast, simple and lightweight Bloom filter library for Python, fully implemented in Rust.

rBloom A fast, simple and lightweight Bloom filter library for Python, fully implemented in Rust. It's designed to be as pythonic as possible, mimicki

Rust crate for interacting with the Windows Packet Filter driver.

NDISAPI-RS NDISAPI-RS is a Rust crate for interacting with the Windows Packet Filter driver. It provides an easy-to-use, safe, and efficient interface

The fastest bloom filter in Rust. No accuracy compromises. Use any hasher.
The fastest bloom filter in Rust. No accuracy compromises. Use any hasher.

b100m-filter The fastest bloom filter in Rust. No accuracy compromises. Use any hasher. Usage # Cargo.toml [dependencies] b100m-filter = "0.3.0" use b

CLI tool that make it easier to perform multiple lighthouse runs towards a single target and output the result in a plotable format.

Lighthouse Aggregator CLI tool that make it easier to perform multiple lighthouse runs towards a single target and output the result in a "plotable" f

UniSBOM is a tool to build a software bill of materials on any platform with a unified data format.

UniSBOM is a tool to build a software bill of materials on any platform with a unified data format. Work in progress Support MacOS Uses system_profile

a command-line tool that transforms a Git repository into a minimal format for ChatGPT queries
a command-line tool that transforms a Git repository into a minimal format for ChatGPT queries

gprepo /dʒiːpiːˈɹi:pi:oʊ/ a command-line tool that transforms a Git repository into a minimal format for ChatGPT queries. Features Excludes LICENSE an

A tool to format codeblocks inside markdown and org documents.
A tool to format codeblocks inside markdown and org documents.

cbfmt (codeblock format) A tool to format codeblocks inside markdown, org, and restructuredtext documents. It iterates over all codeblocks, and format

⚗️ Superfast CLI interface for the conventional commits commit format
⚗️ Superfast CLI interface for the conventional commits commit format

resin ⚗️ Superfast CLI interface for the conventional commits commit format ❓ What is resin? resin is a CLI (command-line interface) tool that makes i

Releases(v0.1.1)
Owner
Ryan Wick
I'm a bioinformatician with a focus on bacterial genome assembly.
Ryan Wick
Rust crate `needleman_wunsch` of the `fasebare` package: reading FASTA sequences, Needleman-Wunsch alignment

fasebare Rust crate needleman_wunsch of the fasebare package: reading FASTA sequences, Needleman-Wunsch alignment. Synopsis The crate needleman_wunsch

Laurent Bloch 2 Nov 19, 2021
Parallel iteration of FASTA/FASTQ files, for when sequence order doesn't matter but speed does

Rust-parallelfastx A truly parallel parser for FASTA/FASTQ files. Principle The input file is memory-mapped then virtually split into N chunks. Each c

Rayan Chikhi 8 Oct 24, 2022
A command-line downloader for sites archived on the Wayback Machine

This is a small command-line utility I wrote to help with browsing archived websites from the Wayback Machine, which can sometimes be pretty slow.

Jonas Schievink 7 Oct 18, 2022
Mini Rust CLI to deploy sites to Netlify using their API

This is a Rust CLI that uses the Netlify API to deploy sites.

Benjamin Lannon 10 May 12, 2022
Just a collection of tiny Rust projects I've did. None warrant a whole repo rn

Daily Rust I try to write some tiny programs daily to gradually improve my Rust skills! Current Program Descriptions first_prog.rs: My first program a

null 3 Nov 5, 2022
Sero is a web server that allows you to easily host your static sites without pain. The idea was inspired by surge.sh but gives you full control.

sero Lightning-fast, static web publishing with zero configuration and full control ?? Table Of Contents ?? Table Of Contents ?? Tools ❓ About The Pro

Dmitry Miasnenko 6 Nov 13, 2023
Simple terminal alignment viewer

Alen Simple terminal sequence alignment viewer. What is Alen? It's a command-like program to view DNA or protein alignments in FASTA formats. Alen is

Jakob Nybo Nissen 51 Dec 19, 2022
Tight Model format is a lossy 3D model format focused on reducing file size as much as posible without decreasing visual quality of the viewed model or read speeds.

What is Tight Model Format The main goal of the tmf project is to provide a way to save 3D game assets compressed in such a way, that there are no not

null 59 Mar 6, 2023
Alexander Mongus is a state-of-the-art filter to sneak amogus characters in pictures

A. Mongus Go to: http://www.lortex.org/amogu/ ??? This is a client-side, Webassembly-based filter to hide amongus characters in your images. Example:

Lucas Pluvinage 3 Apr 16, 2022
Yet Another Kalman Filter Implementation. As well as Lie Theory (Lie group and algebra) on SE(3). [no_std] is supported by default.

yakf - Yet Another Kalman Filter Yet Another Kalman Filter Implementation, as well as, Lie Theory (Lie group, algebra, vector) on SO(3), SE(3), SO(2),

null 7 Dec 1, 2022