-
Notifications
You must be signed in to change notification settings - Fork 11
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
13 changed files
with
165 additions
and
46 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Empty file.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,53 @@ | ||
| from Bio import SeqIO | ||
| from Bio.SeqRecord import SeqRecord | ||
| from Bio.SeqFeature import SeqFeature, FeatureLocation, Seq | ||
| import shutil | ||
| import argparse | ||
| import sys | ||
|
|
||
| def new_reference(referencefile, outgenbank, outfasta, gene): | ||
| ref = SeqIO.read(referencefile, "genbank") | ||
| startofgene = None | ||
| endofgene = None | ||
| for feature in ref.features: | ||
| if feature.type == 'source': | ||
| ref_source_feature = feature | ||
| if feature.type =='gene' or feature.type == 'CDS': | ||
| a = list(feature.qualifiers.items())[0][-1][0] | ||
| if a == gene: | ||
| startofgene = int(list(feature.location)[0]) | ||
| endofgene = int(list(feature.location)[-1])+1 | ||
|
|
||
| # If user provides a --gene 'some name' that is not found, error out as this may indicate that | ||
| # the gene name is misspelled or the user may be using the wrong GenBank file. | ||
| if(gene is not None and startofgene is None and endofgene is None): | ||
| print(f"ERROR: No '{gene}' was found under 'gene' or 'CDS' features in the GenBank file.", file=sys.stderr) | ||
| sys.exit(1) | ||
|
|
||
| record = ref[startofgene:endofgene] | ||
| source_feature = SeqFeature(FeatureLocation(start=0, end=len(record)), type='source', | ||
| qualifiers=ref_source_feature.qualifiers) | ||
| record.features.append(source_feature) | ||
|
|
||
| SeqIO.write(record, outgenbank, 'genbank') | ||
| SeqIO.write(record, outfasta, "fasta") | ||
|
|
||
|
|
||
|
|
||
| if __name__ == '__main__': | ||
| parser = argparse.ArgumentParser( | ||
| description="make new reference depending on whether the entire genome or only part is to be used for the tree", | ||
| formatter_class=argparse.ArgumentDefaultsHelpFormatter | ||
| ) | ||
| parser.add_argument("--reference", required=True, help="GenBank file with reference sequences") | ||
| parser.add_argument("--output-fasta", required=True, help="GenBank new reference file") | ||
| parser.add_argument("--output-genbank", required=True, help="GenBank new reference file") | ||
| parser.add_argument("--gene", help="gene name or genome for entire genome") | ||
| args = parser.parse_args() | ||
|
|
||
| if args.gene=='genome': | ||
| shutil.copy(args.reference, args.output_genbank) | ||
| SeqIO.write(SeqIO.read(args.reference, 'genbank'), args.output_fasta, 'fasta') | ||
| else: | ||
| new_reference(args.reference, args.output_genbank, args.output_fasta, args.gene) | ||
|
|
File renamed without changes.
File renamed without changes.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
File renamed without changes.
File renamed without changes.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,60 @@ | ||
| """ | ||
| This part of the workflow prepares reference files and sequences for constructing the gene phylogenetic trees. | ||
| REQUIRED INPUTS: | ||
| reference = path to reference sequence or genbank | ||
| sequences = path to all sequences from which gene sequences can be extracted | ||
| OUTPUTS: | ||
| gene_fasta = reference fasta for the gene (e.g. E gene) | ||
| gene_genbank = reference genbank for the gene (e.g. E gene) | ||
| sequences = sequences with gene sequences extracted and aligned to the reference gene sequence | ||
| This part of the workflow usually includes the following steps: | ||
| - newreference.py: Creates new gene genbank and gene reference FASTA from the whole genome reference genbank | ||
| - nextclade: Aligns sequences to the reference gene sequence and extracts the gene sequences to ensure the reference files are valid | ||
| See Nextclade or script usage docs for these commands for more details. | ||
| """ | ||
|
|
||
| ruleorder: align_and_extract_E > decompress | ||
|
|
||
| rule generate_E_reference_files: | ||
| """ | ||
| Generating reference files for the E gene | ||
| """ | ||
| input: | ||
| reference = "config/reference_{serotype}_genome.gb", | ||
| output: | ||
| fasta = "results/config/reference_{serotype}_E.fasta", | ||
| genbank = "results/config/reference_{serotype}_E.gb", | ||
| params: | ||
| gene = "E", | ||
| shell: | ||
| """ | ||
| python3 bin/newreference.py \ | ||
| --reference {input.reference} \ | ||
| --output-fasta {output.fasta} \ | ||
| --output-genbank {output.genbank} \ | ||
| --gene {params.gene} | ||
| """ | ||
|
|
||
| rule align_and_extract_E: | ||
| """ | ||
| Cutting sequences to the length of the E gene reference sequence | ||
| """ | ||
| input: | ||
| sequences = "data/sequences_{serotype}.fasta", | ||
| reference = "results/config/reference_{serotype}_E.fasta" | ||
| output: | ||
| sequences = "results/E/sequences_{serotype}.fasta" | ||
| params: | ||
| min_length = 1000, | ||
| shell: | ||
| """ | ||
| nextclade run \ | ||
| -j 1 \ | ||
| --input-ref {input.reference} \ | ||
| --output-fasta {output.sequences} \ | ||
| --min-seed-cover 0.01 \ | ||
| --min-length {params.min_length} \ | ||
| --silent \ | ||
| {input.sequences} | ||
| """ |