oprn_consv_calc
Overview
The operon_conserve_detect.py script takes a reference operon as a list of genome accessions and will determine its structural conservation (i.e. gene order) within a specific taxonomic clade.
Pipeline
Object Definitions
-
GenomeFeature (feature.py)
- Holds information about a feature
- Member variables:
- Genome accession
- Coding start and end
- Strand (+ or -)
- Protein accession
- Locus tag
- Amino acid sequence
- Also requires the following parameters to be passed in when initializing the object:
- Request limit for the Entrez searches
- Sleep time between consecutive Entrez searches
- Main Functions:
- get_intergenic_distance(self, other): Determines the distance to another GenomeFeature
- __str__() and __eq__() are overridden string and equivalence functions
-
AnnotatedHit (feature.py)
- A subclass of GenomeFeature. Contains additional information pertaining to the BLAST hit results
- Member variables:
- The query accession for the BLAST search that resulted in this hit
- The alignment start and end positions
- The percent amino acid identity
- The alignment sequence
- Main Functions:
- fetch_feature(self, record, margin_limit, max_attempts, mult_factor): Determines the feature that corresponds to this hit. It works by finding a feature in the full genome accession record that contains the alignment start and end positions in bounds with a set margin.
- __str__() and __eq__() are overriden string and equivalence functions
-
Operon (operon.py)
- Holds a set of GenomeFeature (and/or by extension AnnotatedHit) objects that have been assigned to the same operon.
- Member variables:
- features[]: A list of GenomeFeature (and/or by extension AnnotatedHit) objects
- genome_accession: The nucleotide accession for the GenomeFragment this operon belongs to
- genome_features[]: A list of all the features in the GenomeFragment this operon belongs to.
- strand: Whether the operon is on the plus or minus strand
- Main Functions:
- add_feature(self, feature): Checks if the feature belongs to the GenomeFragment this operon belongs to, and adds it to the list of features. It then sorts the current list of features from 5' to 3'
- __str__() is a overridden string function
-
GenomeFragment (genome_fragment.py)
- Represents a single nucleotide accession (i.e. plasmid, chromosome, contig, ...etc.)
- Member variables:
- hits[]: The AnnotatedHits associated with this GenomeFragment
- all_features[]: All the GenomeFeatures that are parsed out of the full record for this GenomeFragment
- operons[]: Operon objects associated with this GenomeFragment
- name: The name of this GenomeFragment
- genome_accession: The genome accession for this GenomeFragment
- assembly_accession: The accession for the assembly this nucleotide record belongs to
- taxid: the taxonmic id
- full_record: The entire nucleotide record. Downloaded locally in the cache directory.
- Main Functions:
- fetch_features(self): Obtains all coding features in the genome and saves them as a list of GenomeFeature obejcts
- fetch_record(self): Obtains the full genome record for this nucleotide accession. Reads in from a file or downloads it.
- fetch_hit_features(self, margin_limit, max_attempts, mult_factor): Fetches the features associated with each of the hits associated with this GenomeFragment
- purge_hits(self): Removes any duplicate hits
- add_hit(self, a_hit): Adds an AnnotatedHit object to the current list of hits.
- assemble_operons(self, feature_limit, intergenic_limit): Takes the list of hits and organizes them into putative operons (Operon objects)
-
Species (species.py)
- Holds all the GenomeFragment objects that belong to the same species
- Member variables:
- assembly_accession: The assembly accession for this Species
- species_name: The name of the assembly
- genome_fragments[]: A list of GenomeFragment objects associated with this Species
- genome_fragment_accessions[]: A list of the accessions for all the GenomeFragments in genome_fragments[]
- sim_score: The structural similarity score for this Species relative to the reference
- query_percent_ids[]: The amino acid percent identity for each of the reference genes.
- taxid: The taxonomic id for this assembly
- Main Functions
- add_genome_fragment(self, genome_fragment): Adds a GenomeFragment object to the species
- measure_sim(self, reference_operon): Calculates the structural similarity of the Species relative to the reference operon.
- extract_taxids(self), extract_genome_accessions(self), extract_names(self): These functions set the respective member variables so that the information is stored when the Species is cleaned from memory.
- draw_figure(self, color_code, output_dir): Draw a genome diagram of each of the Operons to SVG files.
- clean(self): Deletes all the GenomeFragments to clear up memory.
Sample Input
See the template input JSON file in /input/.
Parameter | Description |
---|---|
entrez Parameters | |
request_limit |
The max number of attempts that should be made to complete an Entrez request |
sleep_time |
The amount of time that should be waited between consective attempts at an Entrez request (ms) |
email |
Email to use for the Entrez api |
api_key |
Entrez API key |
--- | --- |
blast Parameters | |
tax_include and tax_exclude |
These will restrict the BLAST to include (tax_include ) and exclude (tax_exclude ) specific taxa. |
database |
The database to conduct the BLAST search against. |
blast_type |
Can be set to remote or local . If local , use local_db_path to set the path. |
e-val |
The maximum e-value threshold for the BLAST searches. |
coverage_min |
The coverage threshold to restrict the BLAST results. |
max_hits |
The max hits to recieve from a BLAST search. |
max_attempts |
The maximum number of attempts that should be made to pull all the BLAST hits within the e-value cutoff. This is only applicable is the extensive_search option is set to true . |
search_mult_factor |
The max_hits is increased by this factor through each iteration if the extensive_search option is set to true . |
annotate |
Can be true or false . Determines if the BLAST search will return AnnotatedHit objects. Current functionality requires that it MUST BE true . |
extensive_search |
If true , multiple BLAST searches will be conducted (if neccessary) to pull all hits within the e-value limit. If false , a single BLAST search will be conducted regardless of whether all hits witin the e-value threshold are returned. |
reverse_blast |
If true , all return hits will be tested for true homology by BLAST searching the reference genome with the hit. |
--- | --- |
hit_feature_detection | |
margin_limit |
The acceptable margin between the alignment positions and the positions of the annotated feature. |
max_attempts |
The maximum number of attempts that will be made to identify a feature for a hit. |
mult_factor |
The factor by which the margin_limit will increase for every attempt at identifying the feature for a hit. |
--- | --- |
operon_assembly | |
feature_limit |
Maximum number of features allowed between two hits on an operon. |
intergenic_limit |
Maximum distance in bp allowed between two genes in an operon. |
use_ref_limit and ref_limit_margin |
If true , the intergenic_limit will be determined by the max intergenic distance in the refernce operon multiplied by the margin. |
--- | --- |
thread_limit |
A multithreaded approach is implemented when processing each of the Species objects. This is the maximum number of threads. |
species_percent_id_limit |
For a species to be outputted, it must have at least one hit with an amino acid percent identity above this limit. |
input_records |
The protein accessions for the genes in an operon. The order the genes are representative of the order of the genes in the reference genome. |
reference_genome_accession |
The nucleotide accession of the nucleotide record from the reference speices. This nucleotide record contains the reference operon. |
reference_genome_assembly |
The genome assembly accession from which the reference operon comes from. Needed to conduct the reverse BLAST locally. |
reference_genome_name |
The name of the reference assembly. |
cache_dir |
The cache directory where all nucleotide records will be downloaded locally. |
color_code |
A dictionary that defines the colors for each of the genes in the reference operon. These colors are used to draw the opeon diagram for the putative operons. intergenic refers to features that are inserted into the operon, but were not a hit to any of the reference genes. |
Output
Results are outputted to a CSV and each operon is drawn to a SVG:
Usage
Once the conda
environment is set up, the script can be ran via command line:
python operon_conserve_detect.py [input file name].json
where the input file name is the name of the input file, not including the directory. The input file is expected to be in /input/
.