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

Allows for CDS (as well as gene) features to generate a new gene reference #55

Merged
merged 2 commits into from
Mar 12, 2024
Merged

Allows for CDS (as well as gene) features to generate a new gene reference #55

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
570daa6
Allows for CDS (as well as gene) features to generate a new gene refe…
j23414 Mar 8, 2024
8cd6a13
newreference: If gene is given but not found in genbank, error out
j23414 Mar 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 10 additions & 1 deletion scripts/newreference.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,18 +3,27 @@
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
huddlej marked this conversation as resolved.
Show resolved Hide resolved
for feature in ref.features:
if feature.type == 'source':
ref_source_feature = feature
if feature.type =='gene':
if feature.type =='gene' or feature.type == 'CDS':
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ideally we'd use a centralised GenBank parser, such as augur.utils.load_features. It seem like our parser only considers CDS feature types at the moment, but that could be easily extended.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ohh, I'm open to using augur.utils.load_features. Ah, is there an example script that can be easily referenced?

For extending to "gene feature types" is this basically adding a duplicate code block like the following?

if feat.type=='gene':
    # ... copy internal code from "if feat.type=='CDS'"

I assume if both a "gene" and a "CDS" feature are named "E", one won't overwrite the nucleotide coordinates of the other (although hopefully they match)?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for nextstrain/augur#1435! Upon review, and in order to not hold up progress here (and measles), I'd suggest using the code as it stands and then we can move to using a util function once we're happy with the changes there.

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' is not found, print a warning and use the entire genome.
# Otherwise do not print a warning.
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)
Expand Down