<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Sequence records, part 1
Instructions
Read the chapter 4 and the sections 5.1, 9.3, 9.6, 9.15, and 9.16 of the Biopython tutorial. It is also recommended to read the sections 5.4, 5.5, 9.10, and 9.12, which are touched by this notebook.
Objectives
create and manipulate SeqRecord and SeqFeature objects
search and fetch sequence records from Entrez databases
read sequence records from files
Summary
Sequence records are modelled by the SeqRecord class in the Bio.SeqRecord module. Sequence records contain a sequence and its associated metadata (such as annotations).
Step1:  SeqRecord objects contain metadata
SeqRecord objects contain the following attributes
Step2:  The attributes can be modified as needed.
Step3:  Unstructured annotations are organised into a dictionary by their name (dictionary key) and can contain basically anything.
Step4:  Per-letter annotations are also organised by their name into a dictionary. The annotations must be iterable objects of the same length at the sequence.
Step5:  Features are a list of SeqFeature objects (see below).
Step6:  Cross-references are a list of strings.
Step7:  SeqRecord can be modified (slicing, concatenation, etc.) to produce new records. In general, the letter annotations and features are preserved but annotations and cross-references are not. You must add them manually if needed. Remember to change the id, name, and description of the new record, too.
See the Biopython tutorial for details.
Features are SeqFeature objects
Higher-level information about a sequence is modelled by the SeqFeature class in the Bio.SeqFeature module. A SeqFeature object specifies a region of the sequence and associates some information to the region. The objects contain the following attributes
Step8:  Locations can be specified in various ways
A single position of a sequence can be described by ExactPosition objects (but there are also fuzzy variants). Continuous and discontinuous ranges can be described by FeatureLocation and CompoundLocation objects, respectively. See the Bio.SeqFeature module for details.
Step9:  Locations can be queried for their details.
Step10:  Specified region can be extracted from full sequence
Step11:  Features (or locations, in fact) are also aware of the strand (if applicable).
Step12:  Sequence records can be obtained from NCBI's Entrez databases
The Bio.Entrez contains tools to query Entrez databases and fetch sequence records. They make use of the Entrez Programming Utilities (also known as EUtils), described in detail here. The output of these tools is usually in the XML format, which can be parsed by the utility functions of the Bio.Entrez module.
Important
Step13:  ESearch for searching databases
The esearch function can be used to search and retrieve primary IDs of records. It requires the database and the query as input and accepts some optional arguments as well. See the ESearch documentation for details.
Step14:  The output can be stored to a file for later use, which is a cost-effective approach for sequences in particular.
Step15:  The read function in Bio.Entrez can parse the output of EUtils tools into a dictionary. It expects input that contains exactly one record.
In this case, the output is the query result that contains the retrieved IDs (and some other information).
Step16:  EFetch for retrieving records
The efetch function can be used to retrieve specific records. It requires the database and the list of IDs as input and accepts some optional arguments as well. See the EFetch documentation for details.
The output format can be specified. The available formats depend on the database.
Step17:  Full records can be obtained in the GenBank format in plain text (and in XML).
Step18:  SeqIO produces SeqRecord objects
The Bio.SeqIO module contains functions to parse the retrieved sequence records into SeqRecord objects. The read function parses an input that contains a single record, whereas the parse function returns a generator that produces all records from the input sequentially. They require the source (filename or file-like object) and the format as input.
(Bio.Entrez can parse the XML output of efetch but the result is a dictionary rather than SeqRecord because Bio.Entrez is intended for all Entrez databases and is not specific to sequence data.)
Step19:  Since parse returns a generator, you can process each record at the time without storing all records in memory simultaneously.
Step20:  The amount of information varies between formats. You should use the format that is the most appropriate for your application.
Step21:  The SeqRecord objects can also be written to a file by Bio.SeqIO.
Step22:  Large queries with history server and batch retrieve
NCBI expects that large queries are made by using their support for query history because it saves their computational resources. To use the history, add the usehistory="y" argument to your function call.
Step23:  The result will contain WebEnv and QueryKey information, which replace the list of IDs in efetch. These pieces of information identify the result set of your search.
Step24:  A large number of records should be retrieved in batches, which can be achieved with the retstart and retmax arguments. The former indicates how many records to skip from the beginning and the latter how many records to retrieve.
Note that the example below parses the records directly but, also in this case, the retrieved records could be stored to a file before parsing.
<ASSISTANT_TASK:>
Python Code:
import Bio.SeqRecord as BSR
Explanation: Sequence records, part 1
Instructions
Read the chapter 4 and the sections 5.1, 9.3, 9.6, 9.15, and 9.16 of the Biopython tutorial. It is also recommended to read the sections 5.4, 5.5, 9.10, and 9.12, which are touched by this notebook.
Objectives
create and manipulate SeqRecord and SeqFeature objects
search and fetch sequence records from Entrez databases
read sequence records from files
Summary
Sequence records are modelled by the SeqRecord class in the Bio.SeqRecord module. Sequence records contain a sequence and its associated metadata (such as annotations).
End of explanation
import Bio.Seq as BS
import Bio.Alphabet as BA
# sequence
seq = BS.Seq('MDGEDVQALVIDNGSGMCKA', BA.generic_protein)
# sequence record
record = BSR.SeqRecord(seq)
print(record)
# get sequence from record
print(record.seq)
Explanation: SeqRecord objects contain metadata
SeqRecord objects contain the following attributes:
seq : sequence, typically a Seq object
id : primary identifier, string
name : common name, string
description : short human-readable description, string
letter_annotations : metadata pertaining to single letters
features : structured metadata as SeqFeature objects
annotations : unstructured metadata
dbxrefs : cross-references
SeqRecord objects can be easily created.
End of explanation
# add identifier
record.id = "AC500001"
# add name
record.name = "DUMMY"
# add description
record.description = "Dummy sequence"
print(record)
Explanation: The attributes can be modified as needed.
End of explanation
# add annotation named 'source'
record.annotations['source'] = "Imaginary organism"
print(record)
print(record.annotations['source'])
Explanation: Unstructured annotations are organised into a dictionary by their name (dictionary key) and can contain basically anything.
End of explanation
# add letter annotation named 'structure'
record.letter_annotations['structure'] = '---EEE--EEEEE-EEEEEE'
print(record)
print(record.letter_annotations['structure'])
Explanation: Per-letter annotations are also organised by their name into a dictionary. The annotations must be iterable objects of the same length at the sequence.
End of explanation
# no features yet
print(record.features)
Explanation: Features are a list of SeqFeature objects (see below).
End of explanation
# add cross-references
record.dbxrefs = ['DB:dummy', 'Project:imaginary']
print(record)
Explanation: Cross-references are a list of strings.
End of explanation
import Bio.SeqFeature as BSF
feature = BSF.SeqFeature(BSF.FeatureLocation(0,10), type='domain')
feature.qualifiers['evidence'] = 'experimental'
print(feature)
# type
print(feature.type)
# location
print(feature.location)
# qualifiers
for key, value in feature.qualifiers.items():
    print(key, ':', value)
Explanation: SeqRecord can be modified (slicing, concatenation, etc.) to produce new records. In general, the letter annotations and features are preserved but annotations and cross-references are not. You must add them manually if needed. Remember to change the id, name, and description of the new record, too.
See the Biopython tutorial for details.
Features are SeqFeature objects
Higher-level information about a sequence is modelled by the SeqFeature class in the Bio.SeqFeature module. A SeqFeature object specifies a region of the sequence and associates some information to the region. The objects contain the following attributes:
type: feature type as string
location: sequence region
qualifiers: dictionary of unstructured information
End of explanation
# position 5 of sequence (zero-based indexing)
pos = BSF.ExactPosition(5)
print(repr(pos))
# positions can treated as integers
print(repr(pos + 1))
# region 6-10 of sequence (one-based, inclusive)
# which is 5:10 in Python-style indexing (zero-based, end exclusive)
loc1 = BSF.FeatureLocation(5, 10)
print(loc1)
# region 21-25 (inclusive)
loc2 = BSF.FeatureLocation(20, 25)
print(loc2)
# join of two regions (such as gene with two exons)
loc3 = BSF.CompoundLocation([loc1, loc2])
print(loc3)
Explanation: Locations can be specified in various ways
A single position of a sequence can be described by ExactPosition objects (but there are also fuzzy variants). Continuous and discontinuous ranges can be described by FeatureLocation and CompoundLocation objects, respectively. See the Bio.SeqFeature module for details.
End of explanation
# start position
print(loc3.start)
# end position
print(loc3.end)
# continuous sub-regions
for p in loc3.parts:
    print(p)
# is the position 4 in the region?
print(4 in loc3)
# is the position 8 in the region?
print(8 in loc3)
# is one region within another? (not supported)
try:
    print(loc1 in loc3)
except ValueError as e:
    print(e)
Explanation: Locations can be queried for their details.
End of explanation
# full sequence
seq = BS.Seq('GGACTCTTAGCGGCTCACGCACTTTCTTCCGAAGACGGAACCCG', BA.generic_dna)
# two-exon gene within sequence
feature = BSF.SeqFeature(loc3, type='gene', strand=1)
# the nucleotides specified by the feature
print(repr(feature.extract(seq)))
Explanation: Specified region can be extracted from full sequence
End of explanation
# feature in reverse strand
loc4 = BSF.FeatureLocation(0, 5, strand=-1)
print(loc4)
# strand (1 == forward, -1 == reverse)
print(loc4.strand)
feature = BSF.SeqFeature(loc4, type='gene')
# the five nucleotides specified by the feature
# (note that the output is a reverse complement)
print(repr(feature.extract(seq)))
Explanation: Features (or locations, in fact) are also aware of the strand (if applicable).
End of explanation
import Bio.Entrez as BE
# set email globally
BE.email = 'your.name@example.com'
# remember to set your real email address before continuing to the examples below
SET-YOUR-REAL-EMAIL-ADDRESS-BEFORE-COMMENTING-OUT-THIS-LINE
Explanation: Sequence records can be obtained from NCBI's Entrez databases
The Bio.Entrez contains tools to query Entrez databases and fetch sequence records. They make use of the Entrez Programming Utilities (also known as EUtils), described in detail here. The output of these tools is usually in the XML format, which can be parsed by the utility functions of the Bio.Entrez module.
Important: Be aware that the Entrez utilities impose usage limits and that you can get yourself blacklisted if you violate them. Biopython takes care of keeping an eye on those for you. You should also fill in a non-fake email address. See the Frequency, Timing and Registration of E-utility URL Requests section of the manual.
End of explanation
# search "nucleotide" database for the "accD" gene of "opuntia" organism
handle = BE.esearch(db="nucleotide",
                    term="opuntia[ORGN] accD[gene]",
                    # return ID list in XML format
                    rettype="uilist",
                    retmode="xml")
Explanation: ESearch for searching databases
The esearch function can be used to search and retrieve primary IDs of records. It requires the database and the query as input and accepts some optional arguments as well. See the ESearch documentation for details.
End of explanation
# save to file
with open('results.xml', 'w') as f:
    f.write(handle.read())
Explanation: The output can be stored to a file for later use, which is a cost-effective approach for sequences in particular.
End of explanation
# load from file and parse
with open('results.xml') as f:
    results = BE.read(f)
# (you could also parse 'handle' directly)
# retrieved IDs
print(results['IdList'])
Explanation: The read function in Bio.Entrez can parse the output of EUtils tools into a dictionary. It expects input that contains exactly one record.
In this case, the output is the query result that contains the retrieved IDs (and some other information).
End of explanation
# the ids of three records from the previous search
ids = results['IdList']
# fetch records from "nucleotide" database
handle = BE.efetch(db="nucleotide",
                   id=ids,
                   # return sequences in FASTA format
                   rettype="fasta",
                   retmode="text")
# save to file
with open('sequences.fasta', 'w') as f:
    f.write(handle.read())
# show file content
with open('sequences.fasta') as f:
    print(f.read())
Explanation: EFetch for retrieving records
The efetch function can be used to retrieve specific records. It requires the database and the list of IDs as input and accepts some optional arguments as well. See the EFetch documentation for details.
The output format can be specified. The available formats depend on the database.
End of explanation
# fetch records from "nucleotide" database
handle = BE.efetch(db="nucleotide",
                   id=ids,
                   # return sequences in GenBank plain text format
                   rettype="gb",
                   retmode="text")
# save to file
with open('sequences.gb', 'w') as f:
    f.write(handle.read())
# show file content
with open('sequences.gb') as f:
    print(f.read())
Explanation: Full records can be obtained in the GenBank format in plain text (and in XML).
End of explanation
import Bio.SeqIO as BSIO
Explanation: SeqIO produces SeqRecord objects
The Bio.SeqIO module contains functions to parse the retrieved sequence records into SeqRecord objects. The read function parses an input that contains a single record, whereas the parse function returns a generator that produces all records from the input sequentially. They require the source (filename or file-like object) and the format as input.
(Bio.Entrez can parse the XML output of efetch but the result is a dictionary rather than SeqRecord because Bio.Entrez is intended for all Entrez databases and is not specific to sequence data.)
End of explanation
# parse and process FASTA sequences one at the time
for record in BSIO.parse('sequences.fasta', 'fasta'):
    print(record)
    print()
Explanation: Since parse returns a generator, you can process each record at the time without storing all records in memory simultaneously.
End of explanation
# parse records in GenBank format, which contains more information than FASTA
for record in BSIO.parse('sequences.gb', 'gb'):
    print(record)
    print()
Explanation: The amount of information varies between formats. You should use the format that is the most appropriate for your application.
End of explanation
# write 'record' to file named 'record.gb' in a GenBank format
# (note that the first argument is a list of records)
BSIO.write([record], 'record.gb', 'gb')
Explanation: The SeqRecord objects can also be written to a file by Bio.SeqIO.
End of explanation
# search "nucleotide" database for the "accD" gene of "opuntia" organism using history
handle = BE.esearch(db="nucleotide",
                    term="opuntia[ORGN] accD[gene]",
                    # use history
                    usehistory="y",
                    # return ID list in XML format
                    rettype="uilist",
                    retmode="xml")
results = BE.read(handle)
Explanation: Large queries with history server and batch retrieve
NCBI expects that large queries are made by using their support for query history because it saves their computational resources. To use the history, add the usehistory="y" argument to your function call.
End of explanation
webenv = results['WebEnv']
print(webenv)
querykey = results['QueryKey']
print(querykey)
# fetch records from "nucleotide" database
handle = BE.efetch(db="nucleotide",
                   webenv=webenv,
                   query_key=querykey,
                   rettype="gb",
                   retmode="text")
# save to file
with open('sequences.gb', 'w') as f:
    f.write(handle.read())
Explanation: The result will contain WebEnv and QueryKey information, which replace the list of IDs in efetch. These pieces of information identify the result set of your search.
End of explanation
# function to query "nucleotide" database and retrieve a specific number of records in batches of a given size
def retrieve(query, retrieve_size, batch_size):
    # search using history
    handle = BE.esearch(db="nucleotide",
                        term=query,
                        usehistory="y",
                        rettype="uilist",
                        retmode="xml")
    results = BE.read(handle)
    webenv = results['WebEnv']
    querykey = results['QueryKey']
    # adjust maximum size if less than maximum number of hits available
    retrieve_size = min(int(results['Count']), retrieve_size)
    # split results into batches
    for i in range(0, retrieve_size, batch_size):
        # fetch records
        handle = BE.efetch(db="nucleotide",
                           webenv=webenv,
                           query_key=querykey,
                           # fetch N records starting from i
                           retstart=i,
                           retmax=batch_size,
                           rettype="gb",
                           retmode="text")
        for record in BSIO.parse(handle, 'gb'):
            # acts as an iterator and give records out as they come
            yield record
# example call to function
# (retrieve 10 records in batches of 5; usually these numbers are much larger)
for record in retrieve("opuntia[ORGN] accD[gene]", 10, 5):
    print(record)
    print()
Explanation: A large number of records should be retrieved in batches, which can be achieved with the retstart and retmax arguments. The former indicates how many records to skip from the beginning and the latter how many records to retrieve.
Note that the example below parses the records directly but, also in this case, the retrieved records could be stored to a file before parsing.
End of explanation
<END_TASK>