remoVecSec is a Python toolkit for detecting and removing contamination in assembled genomes, designed to facilitate pre-submission quality control for NCBI WGS and similar repositories. It provides command-line tools and Python modules to automate the screening and cleaning of genome assemblies from common sources of contamination (vectors, adaptors, E. coli, phage, mitochondrial, and more).
- Automated detection and removal of vector, adaptor, and mitochondrial contamination
- BLAST-based screening using NCBI-recommended databases
- Modular Python code for integration into custom pipelines
- Command-line interface for batch processing
- Output of cleaned genome sequences and summary reports
Clone this repository and install dependencies (requires Python 3.7+ and Biopython):
git clone https://github.com/htafer/remoVecSec.git
cd remoVecSec
pip install -r requirements.txtYou will also need BLAST+ and VecScreen installed and accessible in your PATH. See NCBI BLAST+ download and VecScreen.
The main entry point is the remower.py script:
python3 bin/remower.py -g genome.fa -c contam_db -v vec_db -m mito_dbArguments:
--genomefile,-g: Input genome FASTA file (required)--dbvec,-v: VecScreen database (optional)--dbmito,-m: Organelle (mitochondrial) BLAST database (optional)--dbcont,-c: Contaminant BLAST database (optional)--dist,-d: Maximal distance for merging two intervals (default: 50)
Example:
python3 bin/remower.py -g my_genome.fa -c contam_in_euks.fa -v adaptors_for_screening_euks.fa -m mito.ntremoVecSec/removeVec.py: Vector/adaptor detection and removalremoVecSec/removeMito.py: Mitochondrial contamination detectionremoVecSec/removeContaminant.py: General contaminant detection (E. coli, phage, etc.)remoVecSec/removeUtils.py: Utility functions for interval merging and FASTA correction
Each module can be used independently or as part of the main workflow.
- Contaminant Screening:
- Uses BLAST or VecScreen to identify contaminant regions in the genome.
- Interval Merging:
- Overlapping or adjacent contaminant intervals are merged for robust trimming.
- Sequence Correction:
- Contaminated or flagged regions are trimmed or removed from the genome FASTA.
- Reporting:
- Outputs cleaned sequences and a summary of modifications.
- Python 3.7+
- Biopython
- NCBI BLAST+ suite (blastn, makeblastdb)
- VecScreen (for adaptor screening)
This project is licensed under the MIT License. See the LICENSE file for details.
The following sections summarize the NCBI protocols and command-line examples for contaminant, adaptor, mitochondrial, rRNA, and foreign chromosome screening. For full details, see the original NCBI documentation and the remainder of this README.
-
File to screen for the common contaminants in eukaryotic sequences:
ftp://ftp.ncbi.nlm.nih.gov/pub/kitts/contam_in_euks.fa.gz
Contains the cloning artifacts that are likely to show up as contaminants across all eukaryotic species: vector sequences, E.coli genome, phage genomes, bacterial Insertion Sequences and transposons.
-
File to screen for the common contaminants in prokaryotic sequences:
ftp://ftp.ncbi.nlm.nih.gov/pub/kitts/contam_in_prok.fa
Contains phiX174.
These files need to be unzipped and the resulting FASTA sequence files formatted as BLAST databases using the makeblastdb program.
blastn and makeblastdb are contained in the blast+ package which can be installed following the instruction in the BLAST help documents.
"BLAST Command Line Applications User Manual":
https://www.ncbi.nlm.nih.gov/books/NBK279671/
"Standalone BLAST Setup for Windows PC":
https://www.ncbi.nlm.nih.gov/books/NBK52637/
"Standalone BLAST Setup for Unix":
https://www.ncbi.nlm.nih.gov/books/NBK52640/
A BLAST search is run against either the contam_in_euks or contam_in_prok database, depending on the origin of the input sequences. The common contaminant BLAST results are filtered for hits over various length and percent identity cut-offs.
Command line:
-
for screening eukaryotic sequences:
blastn -query _input_fasta_sequences_ -db contam_in_euks -task megablast -word_size 28 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -dust yes -evalue 0.0001 -perc_identity 90.0 -outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" | awk '($3>=98.0 && $4>=50)||($3>=94.0 && $4>=100)||($3>=90.0 && $4>=200)'
OR with an intermediate file, these 2 commands:
blastn -query _input_fasta_sequences_ -db contam_in_euks -task megablast -word_size 28 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -dust yes -evalue 0.0001 -perc_identity 90.0 -outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -out _out_file_
awk '($3>=98.0 && $4>=50)||($3>=94.0 && $4>=100)||($3>=90.0 && $4>=200)' _out_file_
-
for screening prokaryotic sequences:
blastn -query _input_fasta_sequences_ -db contam_in_prok -task megablast -word_size 28 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -dust yes -evalue 0.0001 -perc_identity 90.0 -outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" | awk '($3>=98.0 && $4>=50)||($3>=94.0 && $4>=100)||($3>=90.0 && $4>=200)'
OR with an intermediate file, these 2 commands:
blastn -query _input_fasta_sequences_ -db contam_in_prok -task megablast -word_size 28 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -dust yes -evalue 0.0001 -perc_identity 90.0 -outfmt "7 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -out _out_file_
awk '($3>=98.0 && $4>=50)||($3>=94.0 && $4>=100)||($3>=90.0 && $4>=200)' _out_file_
VecScreen (https://www.ncbi.nlm.nih.gov/tools/vecscreen/) is run against either the adaptors_for_screening_euks.fa database or adaptors_for_screening_proks.fa database, depending on the origin of the input sequences. Hits are filtered to retain only those matches that VecScreen classifies as "Strong" or "Moderate" (see: https://www.ncbi.nlm.nih.gov/tools/vecscreen/about/#Categories).
The adaptors_for_screening databases are available here:
ftp://ftp.ncbi.nlm.nih.gov/pub/kitts/adaptors_for_screening_euks.fa
ftp://ftp.ncbi.nlm.nih.gov/pub/kitts/adaptors_for_screening_proks.fa
These FASTA sequence files need to be formatted as BLAST databases using the makeblastdb program.
The VecScreen standalone program is available here:
ftp://ftp.ncbi.nlm.nih.gov/blast/demo/vecscreen
The script to filter the VecScreen results is here:
ftp://ftp.ncbi.nlm.nih.gov/pub/kitts/VSlistTo1HitPerLine.awk
Command line:
-
for screening eukaryotic sequences:
vecscreen -d adaptors_for_screening_euks.fa -f3 -i _input_fasta_sequences_ -o _vs_output_file_ -
for screening prokaryotic sequences:
vecscreen -d adaptors_for_screening_proks.fa -f3 -i _input_fasta_sequences_ -o _vs_output_file_
Filter out the "Weak" and "Suspect Origin" hits:
VSlistTo1HitPerLine.awk suspect=0 weak=0 _vs_output_file_ > _filtered_vs_output_file_
BLAST is used to screen the input sequences against a database of the mitochondrial genome sequences in the NCBI Reference Sequences (RefSeq) collection.
ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/mito.nt.gz
This file needs to be unzipped and the resulting FASTA sequence file formatted as a BLAST database using the makeblastdb program.
blastn and makeblastdb are contained in the blast+ package (see above).
The BLAST hits to mitochondrial genomes are filtered for hits over 98.6% identity and at least 120 bases long.
blastn -query _input_fasta_sequences -db mito.nt -out % -task megablast -word_size 28 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -dust yes -evalue 0.0001 -perc_identity 98.6 -soft_masking true -outfmt 7 | awk '$4>=120' > _filtered_mito_output_file_
Ribosomal RNA genes are the cause of many false positives because the include some segments that align to distantly related organisms. Segments that match rRNA genes are identified so that such segments are not reported as being foreign.
BLAST is used to screen the input sequences against a database of the rRNA gene sequences .
ftp://ftp.ncbi.nlm.nih.gov/pub/kitts/rrna.gz
This file needs to be unzipped and the resulting FASTA sequence file formatted as a BLAST database using the makeblastdb program.
blastn and makeblastdb are contained in the blast+ package (see above).
The BLAST hits to rRNA genes are filtered for hits over 95% identity and at least 100 bases long.
blastn -query _input_fasta_sequences_ -db rrna -task megablast -template_length 18 -template_type coding -window_size 120 -word_size 12 -xdrop_gap 20 -no_greedy -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -dust yes -evalue 1E-9 -gapextend 2 -gapopen 4 -penalty -4 -perc_identity 95 -reward 3 -soft_masking true -outfmt 7 | awk '$4>=100' > _filtered_mito_output_file_
Screens for matches to chromosome sequences from foreign organisms. Foreign organisms are those that belong to a different taxonomic group compared to the organism whose sequences are being screened. The taxonomic groups are:
arthropoda, chordata, other_metazoa,
viridiplantae, fungi, other_eukaryota,
bacteria, archaea, viruses_and_viroids
Our databases to detect cross-contamination detection are limited to assemblies that have been publicly released in GenBank/ENA/DDBJ and subsequently picked up by RefSeq. Genome centers can do better by augmenting these databases with additional genomes that they have sequenced but which are not yet represented in the RefSeq collection.
- archaea
Query in Nucleotide :
archaea[porgn] AND srcdb_refseq[prop] AND biomol_genomic[prop] AND complete[prop]
- bacteria
Query in Nucleotide :
bacteria[porgn] AND srcdb_refseq[prop] AND biomol_genomic[prop] AND complete[prop]
- fungi
Query in Nucleotide :
fungi[porgn] AND srcdb_refseq[prop] AND biomol_genomic[prop] AND (NC_000000:NC_999999[pacc] OR AC_000000:AC_999999[pacc] OR (NT_000001:NT_999999999[pacc] AND ("chromosome 2L" OR "chromosome 2R" OR "chromosome 3L" OR "chromosome 3R")))
- arthropoda
Query in Nucleotide :
arthropoda[porgn] AND srcdb_refseq[prop] AND biomol_genomic[prop] AND (NC_000000:NC_999999[pacc] OR AC_000000:AC_999999[pacc] OR (NT_000001:NT_999999999[pacc] AND ("chromosome 2L" OR "chromosome 2R" OR "chromosome 3L" OR "chromosome 3R")))
- chordata
Query in Nucleotide :
chordata[porgn] AND srcdb_refseq[prop] AND biomol_genomic[prop] AND (NC_000000:NC_999999[pacc] OR AC_000000:AC_999999[pacc] OR (NT_000001:NT_999999999[pacc] AND ("chromosome 2L" OR "chromosome 2R" OR "chromosome 3L" OR "chromosome 3R")))
- other_metazoa
Query in Nucleotide :
metazoa[porgn] NOT (arthropoda[porgn] OR chordata[porgn]) AND srcdb_refseq[prop] AND biomol_genomic[prop] AND (NC_000000:NC_999999[pacc] OR AC_000000:AC_999999[pacc] OR (NT_000001:NT_999999999[pacc] AND ("chromosome 2L" OR "chromosome 2R" OR "chromosome 3L" OR "chromosome 3R")))
- viridiplantae
Query in Nucleotide :
viridiplantae[porgn] AND srcdb_refseq[prop] AND biomol_genomic[prop] AND (NC_000000:NC_999999[pacc] OR AC_000000:AC_999999[pacc] OR (NT_000001:NT_999999999[pacc] AND ("chromosome 2L" OR "chromosome 2R" OR "chromosome 3L" OR "chromosome 3R")))
- other_eukaryota
Query in Nucleotide :
eukaryota[porgn] NOT (metazoa[porgn] OR fungi[porgn] OR viridiplantae[porgn]) AND srcdb_refseq[prop] AND biomol_genomic[prop] AND (NC_000000:NC_999999[pacc] OR AC_000000:AC_999999[pacc] OR (NT_000001:NT_999999999[pacc] AND ("chromosome 2L" OR "chromosome 2R" OR "chromosome 3L" OR "chromosome 3R")))
- viruses_and_viroids
Query in Nucleotide :
(viruses[porgn] OR viroids[porgn]) AND srcdb_refseq[prop] AND biomol_genomic[prop] AND (NC_000000:NC_999999[pacc] OR AC_000000:AC_999999[pacc] OR (NT_000001:NT_999999999[pacc] AND ("chromosome 2L" OR "chromosome 2R" OR "chromosome 3L" OR "chromosome 3R")))
The FASTA sequence files resulting from these queries are formatted as nine BLAST databases using the makeblastdb program.
Repeats in the input FASTA sequences are soft-masked to lowercase using WindowMasker. Then BLAST hits over 98% identity are generated to the databases for the 8 taxonomic groups to which the organism being screened does not belong.
blastn -query _input_fasta_sequences_ -db _distant_organism_dbs_ -task megablast -word_size 28 -best_hit_overhang 0.1 -best_hit_score_edge 0.1 -dust yes -evalue 0.0001 -min_raw_gapped_score 100 -penalty -5 -perc_identity 98.0 -soft_masking true
The following heuristic rules help to get rid of most false matches.
Contaminant matches from (1) are merged if they are from the same class of sequence (VECTOR, E.coli, IS, PHG) and they overlap or are separated by 50 bases or less.
If the total coverage of contaminant matches from (1) is >75% of the sequence length then flag the sequence as a contaminant to be excluded.
If the contaminant is classed as VECTOR, E.coli, IS:./, PERM:./ or PHG:* and the contaminant location is within 100 bases of the the start or end of the sequence (or gap is the sequence is not contiguous), or within 100 bases of another contaminant match that is at an end, flag the contaminant span for trimming.
If the contaminant is one of the above, and the match is longer than 700 bases flag the contaminant span for trimming.
Other matches may be false alarms. Treat them as suspect spans and reBLAST the hit span plus 10 Kbp of flanking sequence on each side against nr, HTGS, related and unrelated chromosomes (as described below).
Flag all adaptor spans for trimming.
If the total coverage of mitochondrial matches from (3) is >75% of the sequence length then flag the sequence as being mitochindrial sequence to be excluded.
Ignore any matches to chromosomes from unrelated organisms that lie with a region identified as being rRNA genes from (4) (the spans matched in 4 plus 100 bases on both sides). These are likely to be false matches.
Treat other matched spans as suspect and reBLAST the hit span plus 10 Kbp of flanking sequence on each side against nr, HTGS, related and unrelated chromosomes
Spans identified a contamination suspects in the first pass, plus 10 Kbp of flanking sequence on each side (up to the end of the contig), are BLASTed against nr, HTGS, related and unrelated chromosomes to generate additional data for calling contaminants to be excluded or trimmed.
chromosome databases (a) to (i) from (5) above.
ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.*.tar.gz
ftp://ftp.ncbi.nlm.nih.gov/blast/db/htgs.*.tar.gz
The suspect spans are BLASTed against each of the 10 databases.
blastn -query _suspect_spans_plus_flanks_ -db _reblast_db_ -task megablast -dust yes -evalue 1E-9 -searchsp 1000000000 -perc_identity 98.0 -soft_masking true
Automatically exclude sequence contigs that meet all the following criteria:
>60% of length covered with foreign hits, or less than 200 bp that are NOT covered
Each contributing hits must be 100 bp or longer with identity >= 98%
The best match to chromosomes from unrelated organisms is longer than the best match to chromosomes from the related organism group
Some of the other hits may be reviewed manually.
usage: remower.py [-h] --genomefile GENOMEFILE [--dbvec DBVEC] [--dbmito DBMITO] [--dbcont DBCONT] [--dist DIST]
remower is a script that allows to remove contamination in assembled genome. It takes as input a genome file, contamination databases returns on stdout the corrected genome and on stderr warnings regarding vector sequences not removed.
optional arguments: -h, --help show this help message and exit --genomefile GENOMEFILE, -g GENOMEFILE Genome file --dbvec DBVEC, -v DBVEC The vecscreen database --dbmito DBMITO, -m DBMITO The organelle database --dbcont DBCONT, -c DBCONT The contaminant database --dist DIST, -d DIST Maximal distance for merging two intervals
python3 ./remower.py -c contaminationdb -v vectordb -m mitodb
--genomefile genome.fa
#+ENDSRC