Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add "--start" and "--end" arguments to newreference.py to allow for creating subgenic trees #58

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add "--start" and "--end" arguments to newreference.py to allow for creating subgenic trees #58

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
22 changes: 20 additions & 2 deletions scripts/newreference.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
import argparse
import sys

def new_reference(referencefile, outgenbank, outfasta, gene):
def new_reference(referencefile, outgenbank, outfasta, gene, start, end):
ref = SeqIO.read(referencefile, "genbank")
startofgene = None
endofgene = None
Expand All @@ -25,10 +25,20 @@ def new_reference(referencefile, outgenbank, outfasta, gene):
sys.exit(1)

record = ref[startofgene:endofgene]
# Allows for subgenic regions
if (gene is not None and start is not None and end is not None):
record = record[int(start):int(end)]

source_feature = SeqFeature(FeatureLocation(start=0, end=len(record)), type='source',
qualifiers=ref_source_feature.qualifiers)
record.features.append(source_feature)

# For subgenic regions, adds a CDS feature {gene}_{start}_{end}
if(gene is not None and start is not None and end is not None):
gene = gene + "_" + start + "_" + end
record.features.append(SeqFeature(FeatureLocation(start=0, end=len(record)), type='CDS',
qualifiers={'gene': gene}))

SeqIO.write(record, outgenbank, 'genbank')
SeqIO.write(record, outfasta, "fasta")

Expand All @@ -43,11 +53,19 @@ def new_reference(referencefile, outgenbank, outfasta, gene):
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")
parser.add_argument("--start", help="Start 0-based position relative to the gene (requires --gene argument)")
parser.add_argument("--end", help="End 0-based position relative to the gene (requires --gene argument)")
args = parser.parse_args()

if args.gene=='genome':
# Check if start and end are specified here
if (args.start is not None and args.end is not None):
print(f"ERROR: --start '{args.start}' --end '{args.end}' is not supported for full genome.", file=sys.stderr)
sys.exit(1)


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)
new_reference(args.reference, args.output_genbank, args.output_fasta, args.gene, args.start, args.end)