bpy27. Counting Letters of a biological sequence

We use the createRecords function, that we had used in the last example, by using an appropriate import.

We also define a new function, count, which will count unique letters and return their number of occurrences. This function, in turn, can be used by other modules if they import it.

In the function count, we try to add 1 to number already there, for the appropriate dictionary entry. Should that entry not exist, it is created and set to 1.

# bpy27.py
from __future__ import print_function, division
from bpy26 import createRecords

def count(seq):
    count = {}
    for i in seq:
        try: count[i] = count[i] + 1
        except: count[i] = 1
    return count
    
if __name__ == '__main__':
    records = createRecords('data')
    for record in records:
        print('Entry Name:',record.entry_name)
        print('Sequence counts:')
        print(count(record.sequence))

#Entry Name: IGF1R_HUMAN
#Sequence counts:
#{'A': 72, 'C': 44, 'E': 113, 'D': 62, 'G': 91, 'F': 47,
# 'I': 73, 'H': 23, 'K': 69, 'M': 39, 'L': 116, 'N': 85,
# 'Q': 35, 'P': 81, 'S': 97, 'R': 80, 'T': 69, 'W': 21,
# 'V': 87, 'Y': 63}
#Entry Name: IGF1R_MOUSE
#Sequence counts:
#{'A': 68, 'C': 44, 'E': 114, 'D': 63, 'G': 91, 'F': 48,
# 'I': 73, 'H': 24, 'K': 67, 'M': 40, 'L': 115, 'N': 87,
# 'Q': 36, 'P': 83, 'S': 92, 'R': 82, 'T': 73, 'W': 22,
# 'V': 90, 'Y': 61}

bpy26. Parsing SwissProt files using Biopython

In the last example, we saved the insulin proteins to 'data' subfolder.

We have a function that takes one parameter, the subfolder where *.txt files are stored, each corresponding to a UniProtKB text record.

We could have given our files a different extension during saving, for example .dat or .swiss which would require that term in the list comprehension. We do not have to have a filter term if only UniProtKB records are in the subfolder; in this case we have fils = os.listdir(fol)

Several attributes are printed of the records.

# bpy26.py
from __future__ import print_function, division
import os 
from Bio import SwissProt

def createRecords(fol):
    records = []
    fils = [fil for fil in os.listdir(fol) if fil.endswith('.txt')]
    for fil in fils:
        handle = open(fol + '/' + fil)
        record = SwissProt.read(handle)
        records.append(record)
        handle.close()
    return records

if __name__ == '__main__':
    records = createRecords('data')
    for record in records:
        print('Entry Name:',record.entry_name)
        print('Organism:',record.organism)
        print('Length:',record.sequence_length)
        first_crossref = record.cross_references[0]
        print('First cross ref:')
        for i in first_crossref:
            print('\t',i)
        print()
    
#Entry Name: IGF1R_HUMAN
#Organism: Homo sapiens (Human).
#Length: 1367
#First cross ref:
#         EMBL
#         X04434
#         CAA28030.1
#         -
#         mRNA
#
#Entry Name: IGF1R_MOUSE
#Organism: Mus musculus (Mouse).
#Length: 1373
#First cross ref:
#         EMBL
#         AF056187
#         AAC12782.1
#         -
#         mRNA

bpy25. Saving Proteins from UniProtKB using Biopython

We save insulin proteins, from the protein database UniProtKB, for both human and mouse.

The accession strings are stored in the accessions list.

Later, we will parse these files. However, it is better to save the files, rather than always loading them. We save them to two separate files in the subfolder 'data'.

# bpy25.py
import os 
from Bio import ExPASy
accessions = ['P08069','Q60751']
records = []
for accession in accessions:
    handle = ExPASy.get_sprot_raw(accession)
    fpath = os.path.join('data',accession+'.txt')
    with open(fpath,'w') as fout:
        fout.write(handle.read())

bpy24. The Smith-Waterman algorithm in Python

The Smith–Waterman algorithm is used for local alignment of two sequences. This algorithm is described at Wikipedia. We use the two sequences from that link.

The last example, bpy23.py, is modified. There is no initialization as the fill for first row and first column is zero, and np.zeros is used to construct the matrices. For the first step, we find the matrices H and T. In the Wikipedia page, you can see the arrows, for the pointer matrix T, which makes the traceback path clear. In the program the various directions are labelled as integers, such as 3 for diagonal.

Here the high score, 12, happens to be the last column of last row. In general, this will not be true and the highest scoring entry could be located anywhere in the matrix. From the highest entry, we proceed to 0, using the pointer information in T. At each step to 0, we record the letter of sequence 1 and sequence2. There might be cases it is null, where a gap is recorded for the particular entry. This is shown in the printouts, where sometimes the string '-' is added.

# bpy24.py

from __future__ import division, print_function
import numpy as np

pt ={'match': 2, 'mismatch': -1, 'gap': -1}

def mch(alpha, beta):
    if alpha == beta:
        return pt['match']
    elif alpha == '-' or beta == '-':
        return pt['gap']
    else:
        return pt['mismatch']

def water(s1, s2):
    m, n = len(s1), len(s2)
    H = np.zeros((m+1, n+1))    
    T = np.zeros((m+1, n+1))
    max_score = 0
    # Score, Pointer Matrix
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            sc_diag = H[i-1][j-1] + mch(s1[i-1], s2[j-1])
            sc_up = H[i][j-1] + pt['gap']
            sc_left = H[i-1][j] + pt['gap']
            H[i][j] = max(0,sc_left, sc_up, sc_diag)
            if H[i][j] == 0: T[i][j] = 0
            if H[i][j] == sc_left: T[i][j] = 1
            if H[i][j] == sc_up: T[i][j] = 2
            if H[i][j] == sc_diag: T[i][j] = 3
            if H[i][j] >= max_score:
                max_i = i
                max_j = j
                max_score = H[i][j];
    
    print('H=\n',H,'\n')
    print('T=\n',T,'\n')
    align1, align2 = '', ''
    i,j = max_i,max_j
    
    #Traceback
    while T[i][j] != 0:
        if T[i][j] == 3:
            a1 = s1[i-1]
            a2 = s2[j-1]
            i -= 1
            j -= 1
        elif T[i][j] == 2:
            a1 = '-'
            a2 = s2[j-1]
            j -= 1
        elif T[i][j] == 1:
            a1 = s1[i-1]
            a2 = '-'
            i -= 1
        print('%s ---> a1 = %s\t a2 = %s\n' % ('Add',a1,a2))
        align1 += a1
        align2 += a2

    align1 = align1[::-1]
    align2 = align2[::-1]
    sym = ''
    iden = 0
    for i in range(len(align1)):
        a1 = align1[i]
        a2 = align2[i]
        if a1 == a2:                
            sym += a1
            iden += 1
        elif a1 != a2 and a1 != '-' and a2 != '-': 
            sym += ' '
        elif a1 == '-' or a2 == '-':          
            sym += ' '
    
    identity = iden / len(align1) * 100
    print('Identity = %f percent' % identity)
    print('Score =', max_score)
    print(align1)
    print(sym)
    print(align2)

if __name__ == '__main__':
    water('AGCACACA','ACACACTA')

#H=
# [[  0.   0.   0.   0.   0.   0.   0.   0.   0.]
# [  0.   2.   1.   2.   1.   2.   1.   0.   2.]
# [  0.   1.   1.   1.   1.   1.   1.   0.   1.]
# [  0.   0.   3.   2.   3.   2.   3.   2.   1.]
# [  0.   2.   2.   5.   4.   5.   4.   3.   4.]
# [  0.   1.   4.   4.   7.   6.   7.   6.   5.]
# [  0.   2.   3.   6.   6.   9.   8.   7.   8.]
# [  0.   1.   4.   5.   8.   8.  11.  10.   9.]
# [  0.   2.   3.   6.   7.  10.  10.  10.  12.]] 
#
#T=
# [[ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
# [ 0.  3.  2.  3.  2.  3.  2.  2.  3.]
# [ 0.  1.  3.  1.  3.  1.  3.  3.  1.]
# [ 0.  1.  3.  2.  3.  2.  3.  2.  2.]
# [ 0.  3.  1.  3.  2.  3.  2.  2.  3.]
# [ 0.  1.  3.  1.  3.  2.  3.  2.  2.]
# [ 0.  3.  1.  3.  1.  3.  2.  2.  3.]
# [ 0.  1.  3.  1.  3.  1.  3.  2.  2.]
# [ 0.  3.  1.  3.  1.  3.  1.  3.  3.]] 
#
#Add ---> a1 = A  a2 = A
#
#Add ---> a1 = -  a2 = T
#
#Add ---> a1 = C  a2 = C
#
#Add ---> a1 = A  a2 = A
#
#Add ---> a1 = C  a2 = C
#
#Add ---> a1 = A  a2 = A
#
#Add ---> a1 = C  a2 = C
#
#Add ---> a1 = G  a2 = -
#
#Add ---> a1 = A  a2 = A
#
#Identity = 77.777778 percent
#Score = 12.0
#AGCACAC-A
#A CACAC A
#A-CACACTA

by23. Needleman–Wunsch algorithm in Python

The Needleman–Wunsch algorithm is used for global alignment of two sequences. The best global alignment can be scored, and we may find the number of identities..

In the pt dictionary, the three scores for match, mismatch and gap are given.

The first two steps in the algorithm are initialization and fill, which will create the score matrix.

The last column on the last row shows the final score, and we start traceback to first column in first row, and in so doing find the best scoring alignment.

# bpy23.py
from __future__ import print_function, division
import numpy as np

pt ={'match': 1, 'mismatch': -1, 'gap': -1}

def mch(alpha, beta):
    if alpha == beta:
        return pt['match']
    elif alpha == '-' or beta == '-':
        return pt['gap']
    else:
        return pt['mismatch']

def needle(s1, s2):
    m, n = len(s1), len(s2)
    score = np.zeros((m+1, n+1))
    
    #Initialization
    for i in range(m+1):
        score[i][0] = pt['gap'] * i
    for j in range(n+1):
        score[0][j] = pt['gap'] * j
    
    #Fill
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            diag = score[i-1][j-1] + mch(s1[i-1], s2[j-1])
            delete = score[i-1][j] + pt['gap']
            insert = score[i][j-1] + pt['gap']
            score[i][j] = max(diag, delete, insert)

    print('score matrix = \n%s\n' % score)
    align1, align2 = '', ''
    i,j = m,n
    
    #Traceback
    while i > 0 and j > 0:
        score_current = score[i][j]
        score_diag = score[i-1][j-1]
        score_left = score[i][j-1]
        score_up = score[i-1][j]
        
        print('score_current: ',score_current)
        print('score_diag: ',score_diag)
        print('score_left: ',score_left)
        print('score_up: ',score_up)

        if score_current == score_diag + mch(s1[i-1], s2[j-1]):
            print('diag')
            a1,a2 = s1[i-1],s2[j-1]
            i,j = i-1,j-1
        elif score_current == score_up + pt['gap']:
            print('up')
            a1,a2 = s1[i-1],'-'
            i -= 1
        elif score_current == score_left + pt['gap']:
            print('left')
            a1,a2 = '-',s2[j-1]
            j -= 1
        print('%s ---> a1 = %s\t a2 = %s\n' % ('Add',a1,a2))
        align1 += a1
        align2 += a2
            

    while i > 0:
        a1,a2 = s1[i-1],'-'
        print('%s ---> a1 = %s\t a2 = %s\n' % ('Add',a1,a2))
        align1 += a1
        align2 += a2
        i -= 1
        
    while j > 0:
        a1,a2 = '-',s2[j-1]
        print('%s --> a1 = %s\t a2 = %s\n' % ('Add',a1,a2))
        align1 += a1
        align2 += a2
        j -= 1
    
    align1 = align1[::-1]
    align2 = align2[::-1]
    seqN = len(align1)
    sym = ''
    seq_score = 0
    ident = 0
    for i in range(seqN):
        a1 = align1[i]
        a2 = align2[i]
        if a1 == a2:
            sym += a1
            ident += 1
            seq_score += mch(a1, a2)
    
        else: 
            seq_score += mch(a1, a2)
            sym += ' '
        
    ident = ident/seqN * 100
    
    print('Identity = %2.1f percent' % ident)
    print('Score = %d\n'% seq_score)
    print(align1)
    print(sym)
    print(align2)

if __name__ == '__main__':
    needle("TCGA","ACCGT")

#score matrix = 
#[[ 0. -1. -2. -3. -4. -5.]
# [-1. -1. -2. -3. -4. -3.]
# [-2. -2.  0. -1. -2. -3.]
# [-3. -3. -1. -1.  0. -1.]
# [-4. -2. -2. -2. -1. -1.]]
#
#score_current:  -1.0
#score_diag:  0.0
#score_left:  -1.0
#score_up:  -1.0
#diag
#Add ---> a1 = A  a2 = T
#
#score_current:  0.0
#score_diag:  -1.0
#score_left:  -1.0
#score_up:  -2.0
#diag
#Add ---> a1 = G  a2 = G
#
#score_current:  -1.0
#score_diag:  -2.0
#score_left:  0.0
#score_up:  -3.0
#diag
#Add ---> a1 = C  a2 = C
#
#score_current:  -2.0
#score_diag:  -1.0
#score_left:  -1.0
#score_up:  -2.0
#diag
#Add ---> a1 = T  a2 = C
#
#Add --> a1 = -   a2 = A
#
#Identity = 40.0 percent
#Score = -1
#
#-TCGA
#  CG 
#ACCGT

bpy22. Using ECitMatch in Python

If we have these 5 items for an article: (1) Journal Name, (2) Year (3) Volume, (4) Page, (5) Author, we can find the pubmed id for the article using ecitmatch. This will rarely be needed as esearch will usually find the information. However parsing the html string from get requests is important in general.

We do not have a built-in Bio.Entrez utility to deal with ecitmatch yet. But esearch can also return pubmed IDs, and the search term may use limits by Author, by Journal, by Date, etc.

The tuples journal, year, vol, page, author contain strings for the first and second article. We replace spaces with +, add pipes between different items and join the two terms using '%0D'

# bpy22.py
from __future__ import print_function, division
import urllib2

ids = 'article1','article2'
journal = "Front Cell Neurosci","J Biomed Sci"
year = "2014","2013"
vol = "8","20"
page = "47","92"
author = "Iyengar BR", "Skipper KA"
base = ("http://eutils.ncbi.nlm.nih.gov/entrez/"
        "eutils/ecitmatch.cgi?db=pubmed&retmode=xml&bdata=")
s = [None]*len(ids)
for i in range(len(ids)):
    fmt = (journal[i].replace(' ','+'), year[i], vol[i],
           page[i], author[i].replace(' ','+'), ids[i])
    s[i] = '%s|%s|%s|%s|%s|%s' % fmt
addr = base + '%0D'.join(s)    
response = urllib2.urlopen(addr)
html = response.read()
pIDs = []
for i in html.split('\n'):
    if i == '': continue
    pIDs.append(i.split('|')[-1])
print("pIDs=\n",pIDs)

#pIDs=
# ['24605084', '24320156'] 

bpy21. Using elink EUtil in Biopython

We can use Bio.Entrez.elink to find related IDs from one database to another, or from same database. Thus elink requires we give name of database from where the IDs are located. We also need name of database where the new IDs are located. Finally we have to provide a string, or a list of strings, for the IDs.

In the example, the database from is 'protein'. We have three IDs of tyrosine kinases from human, cow and mouse.

elink returns the IDs in the 'gene' database for these proteins.

# bpy21.py
from __future__ import print_function
from Bio import Entrez

Entrez.email = "A.N.Other@example.com"
IDs = ['15718680','157427902','119703751']
handle = Entrez.elink(dbfrom = 'protein', db='gene', id = IDs)
values = Entrez.read(handle)
link_IDs = []
for value in values:
    for linkset in value['LinkSetDb']:
        for link in linkset['Link']:
            ID = link['Id']
            print('Link:', ID)
            link_IDs.append(ID)
print("link_ids=\n",link_IDs)

#Link: 3702
#Link: 522311
#Link: 16428
#link_ids=
# ['3702', '522311', '16428']

bpy20. Using espell EUtil in Biopython

We can use Bio.Entrez.espell to find correct spelling for a term. The default database is Pubmed.

In the dictionary returned by Entrez.read(), there are 2 keys. The search term has the key of "Query" and the corrected term has the key "CorrectedQuery".

# bpy20.py
from __future__ import print_function, division
from Bio import Entrez
Entrez.email = "Your.Name.Here@example.org"
handle = Entrez.espell(term="cance")
record = Entrez.read(handle)
print("Query:",record["Query"])
print("Corrected Query:",record["CorrectedQuery"])

#Query: cance
#Corrected Query: cancer

bpy19. Using epost EUtil in Biopython

We can use Bio.Entrez.epost to post IDs to NCBI, so we may use them later. This example is similar to the last, except now we do not use the usehistory='y' keyword.

After getting a list of IDs using esearch, with the same search term as in last example, we selected a random subset of the IDs, of length 500, and posted them.

Then, they were used in a Bio.Entrez.esummary.

# bpy19.py
from __future__ import print_function
from Bio import Entrez
import numpy as np

Entrez.email = "A.N.Other@example.com"
term = ('"mitochondria"[MeSH Terms] OR "mitochondria"[All Fields])'
        ' AND (Review[ptyp] AND "loattrfree full text"[sb] AND'
        ' ("2012/01/01"[PDAT] : "2014/12/31"[PDAT]))')
handle = Entrez.esearch("pubmed", retmax = 1500, term = term)
values = Entrez.read(handle)
print("Count:",values['Count'])
print("len(IdList)",len(values['IdList']))
randIDs = np.random.choice(values['IdList'],500,replace=False)
print("len(randIDs)",len(randIDs))
handle1 = Entrez.epost("pubmed",id=','.join(randIDs))
values1 = Entrez.read(handle1)
print("QueryKey:", values1['QueryKey'])
print("len(WebEnv):", len(values1['WebEnv']))
summary = Entrez.esummary(db="pubmed", retmax=500,
                          webenv=values1["WebEnv"],
                          query_key=values1["QueryKey"])
records = Entrez.read(summary)
for i,record in enumerate(records):
    if i%100 != 0: continue
    print('%d. %s' % (i+1,record['Title']))
    print('')

#Count: 1223
#len(IdList) 1223
#len(randIDs) 500
#QueryKey: 1
#len(WebEnv): 77
#1. Mitochondrial poly(A) polymerase and polyadenylation.
#
#101. The role of mitochondrial dysfunction in sepsis-induced
# multi-organ failure.
#
#201. The involvement of the sigma-1 receptor in neurodegeneration
# and neurorestoration.
#
#301. Provitamin A metabolism and functions in mammalian biology.
#
#401. Peroxisome proliferator activated receptor α ligands as
# anticancer drugs targeting mitochondrial metabolism.

bpy18. Posting IDs in a NCBI EUtil using Biopython

We can use the keyword parameter usehistory='y' to Bio.Entrez.esearch.

We searched for mitochondria review articles, that have free full text available, from years 2012 through 2014.

Because of the additional parameter, two new keywords were returned and we can use in a subsequent eUtil. We print the title of a few records, from the 1223 that were returned.

# bpy18.py

from __future__ import print_function
from Bio import Entrez
Entrez.email = "A.N.Other@example.com"
term = ('"mitochondria"[MeSH Terms] OR "mitochondria"[All Fields])'
        ' AND (Review[ptyp] AND "loattrfree full text"[sb] AND'
        ' ("2012/01/01"[PDAT] : "2014/12/31"[PDAT]))')
handle = Entrez.esearch("pubmed", retmax = 1500, term = term,
                        usehistory = 'y' )

values = Entrez.read(handle)
print("Count:",values['Count'])
print("len(IdList)",len(values['IdList']))
print("QueryKey:", values['QueryKey'])
print("len(WebEnv):", len(values['WebEnv']))

summary = Entrez.esummary(db="pubmed", retmax=values['Count'],
                          webenv=values["WebEnv"],
                          query_key=values["QueryKey"])
records = Entrez.read(summary)

for i,record in enumerate(records):
    if i%300 != 0: continue
    print('%d. %s' % (i+1,record['Title']))
    print('')

#Count: 1223
#len(IdList) 1223
#QueryKey: 1
#len(WebEnv): 78
#1. The role of SIGMAR1 gene mutation and mitochondrial dysfunction
# in amyotrophic lateral sclerosis.
#
#301. You are what you eat: multifaceted functions of autophagy
# during elegans development.
#
#601. Mitochondrial deficiency in Cockayne syndrome.
#
#901. Transcriptional repression of mitochondrial function in
# aging: a novel role for the silencing mediator of retinoid
# and thyroid hormone receptors co-repressor.
#
#1201. Signal transduction by mitochondrial oxidants.

bpy17. Using efetch EUtil in Biopython

We can use Bio.Entrez.efetch to get records from a database, given an ID or a set of IDs. We use the set of IDs that we used earlier in bpy14.

There are two programs here. The first file saves the records to a text file. The rettype and retmode key arguments are given at the NCBI Entrez site. The variable out is the string corresponding to the text file returned. This particular combination is for getting a text file in Medline format.

The second program uses Bio.Medline.parse to read the Medline file. Here we print the Title and first 200 characters of Abstract.

# bpy17a.py
from __future__ import print_function
from Bio import Entrez
Entrez.email = "Your.Name.Here@example.org"
handle = Entrez.efetch(db="pubmed", id="25741283,25798216",
                       rettype ='medline', retmode = 'text')
out = handle.read()
handle.close()
with open("bpy17.txt","w") as fout:
    fout.write(out)

# bpy17b.py
from __future__ import print_function
from Bio import Medline
fin = open('bpy17.txt')
records = Medline.parse(fin)

for record in records:
    print('Title:',record['TI'])
    print('Abstract\n',record['AB'][:200], end = '')
    print('...\n')
    
fin.close()

#Title: Nitric oxide and mitochondria in metabolic syndrome.
#Abstract
# Metabolic syndrome (MS) is a cluster of metabolic disorders
# that collectively increase the risk of cardiovascular disease.
# Nitric oxide (NO) plays a crucial role in the pathogeneses of
# MS components a...
#
#Title: Maternal ancestry and population history from whole
# mitochondrial genomes.
#Abstract
# MtDNA has been a widely used tool in human evolutionary and
# population genetic studies over the past three decades. Its
# maternal inheritance and lack of recombination have offered
# the opportunity to e...

bpy16. Using esearch EUtil in Biopython

We can use Bio.Entrez.esearch to find IDs from a database. Besides the database, we have to give the search term. We may, additionally, give the maximum number of IDs to return, as a keyword argument.

20 IDs are returned since that is the default value for the maximum number of IDs returned, that is, RetMax is equal to 20. We can also see the total number of matches.

We may search any of the fields, however here we search in [All Fields].

# bpy16.py
from __future__ import print_function, division
from Bio import Entrez
Entrez.email = "Your.Name.Here@example.org"
handle = Entrez.esearch("omim","cancer")
record = Entrez.read(handle)
handle.close()
print("There are a total of %s results." % record['Count'])
r = record['IdList']
NumIds = len(r)
print("%d Ids were returned. They are:" % NumIds)
for i in range(0,NumIds,4):
    r_fmt = i+1,r[i],i+2,r[i+1],i+3,r[i+2],i+4,r[i+3]
    print("%d: %s\t%d: %s\t%d: %s\t%d: %s" % r_fmt)

#There are a total of 2802 results.
#20 Ids were returned. They are:
#1: 616093       2: 616092       3: 616073       4: 616043
#5: 616041       6: 616016       7: 616015       8: 616000
#9: 615993       10: 615977      11: 615943      12: 615930
#13: 615920      14: 615909      15: 615908      16: 615874
#17: 615869      18: 615854      19: 615845      20: 615825

bpy15. Writing a query function using Biopython

Using Bio.Entrez.egquery, that we used earlier in bpy12.py, we will encapsulate the table creation into a function, which will write that table to a csv file.

Additionally, we can give a boolean value to indicate, that if additionally, we want to print or not.

In later examples, we will use it as an import. However, it is good to have test code, should the module be run directly. The test code just prints the table returned for 'asthma' and saves result to count.csv. We will note that there are more Pubmed and Pubmed Central records than in our previous search.

# query.py

from __future__ import print_function, division
from Bio import Entrez

def query(term, fname, to_print):
    Entrez.email = "Your.Name.Here@example.org"
    handle = Entrez.egquery(term=term)
    record = Entrez.read(handle)
    handle.close()
    X=[]
    Y=[]
    out=[]
    for row in record["eGQueryResult"]:
        for i in row:
            if i == 'DbName':
               Y.append(row[i])
            elif i=='Count':
                X.append(row[i])
    for i in range(len(X)):
        if int(X[i])==0:
            continue
        out.append("%s,%s\n" % (Y[i],X[i]))
    if fname == None:
        pass
    else:
        with open(fname+'.csv','w') as fout:
            fout.writelines(out)
    if to_print:
        for i in out:
            print(i, end = '')
  
if __name__ == '__main__':
    query('asthma','count',True)
    
#    pubmed,149164
#    pmc,106784
#    mesh,11
#    books,7165
#    pubmedhealth,1461
#    omim,270
#    ncbisearch,75
#    nuccore,32446
#    nucgss,91
#    nucest,210
#    protein,2559
#    genome,1
#    structure,390
#    dbvar,9
#    gene,1036
#    sra,8245
#    biosystems,125
#    unigene,1
#    cdd,22
#    popset,7
#    geoprofiles,731699
#    gds,4417
#    homologene,18
#    pccompound,122
#    pcsubstance,2125
#    pcassay,3515
#    nlmcatalog,2610
#    probe,34
#    gap,2387
#    bioproject,211
#    biosample,12542  

bpy14. Using esummary EUtil in Biopython

We can use Bio.Entrez.esummary to find summary (not details) of records in a database. The parameter are the ids of records we want to retrieve as well as the database name. If we want information from more than one ID, we can pass the ids as a comma delimited string.

You should give your email for Entrez.email. Biopython returns a list of length corresponding to the number of ids that are provided in the id string.

For the 2 ids, we get title of article, the authors, and journal name. We will have to find ids by using other Entrez eUtils.

# bpy14.py
from __future__ import print_function
from Bio import Entrez
Entrez.email = "Your.Name.Here@example.org"
handle = Entrez.esummary(db="pubmed", id="25741283,25798216")
record = Entrez.read(handle)
handle.close()
print('type(record) =',type(record))
print('len(record) =',len(record))
for i in range(len(record)):
    print("record:",i)
    print("title:",record[i]["Title"])
    for author in record[i]["AuthorList"]:
        print("author:",author)
    print("journal:",record[i]["FullJournalName"])
    print("")

#type(record) = <class 'Bio.Entrez.Parser.ListElement'>
#len(record) = 2
#record: 0
#title: Nitric oxide and mitochondria in metabolic syndrome.
#author: Litvinova L
#author: Atochin DN
#author: Fattakhov N
#author: Vasilenko M
#author: Zatolokin P
#author: Kirienkova E
#journal: Frontiers in physiology
#
#record: 1
#title: Maternal ancestry and population history from whole
# mitochondrial genomes.
#author: Kivisild T
#journal: Investigative genetics

bpy13. Using einfo EUtil in Biopython

We can use Bio.Entrez.einfo to find information about databases. If we call einfo without any parameter, it will return the name of all the databases. However, if we put in a particular db, such as 'pubmed', it will return information about this database.

For Entrez.email, you should give your email. Biopython returns a dictionary of length 1, with the result in key 'DbInfo". Some information about the database is printed such as its name and count.

The dictionary within the result is iterated over (with i being the current subkey). We list all the field names and their description.

# bpy13.py
from __future__ import print_function, division
from Bio import Entrez

Entrez.email = "Your.Name.Here@example.org"
handle = Entrez.einfo(db="pubmed")
record = Entrez.read(handle)
handle.close()
db = record['DbInfo']
name = db['MenuName']
dbName = db['DbName']
desc = db['Description']
count = db['Count']
print(" - ".join([name,dbName,desc,count]))
for i in db['FieldList']:
    print(i['Name'],end='\t')
    print(i['Description'])
    
#PubMed - pubmed - PubMed bibliographic record - 24738242
#ALL     All terms from all searchable fields
#UID     Unique number assigned to publication
#FILT    Limits the records
#TITL    Words in title of publication
#WORD    Free text associated with publication
#MESH    Medical Subject Headings assigned to publication
#MAJR    MeSH terms of major importance to publication
#AUTH    Author(s) of publication
#JOUR    Journal abbreviation of publication
#AFFL    Author's institutional affiliation and address
#ECNO    EC number for enzyme or CAS registry number
#SUBS    CAS chemical name or MEDLINE Substance Name
#PDAT    Date of publication
#EDAT    Date publication first accessible through Entrez
#VOL     Volume number of publication
#PAGE    Page number(s) of publication
#PTYP    Type of publication (e.g., review)
#LANG    Language of publication
#ISS     Issue number of publication
#SUBH    Additional specificity for MeSH term
#SI      Cross-reference from publication to other databases
#MHDA    Date publication was indexed with MeSH terms
#TIAB    Free text associated with Abstract/Title
#OTRM    Other terms associated with publication
#INVR    Investigator
#COLN    Corporate Author of publication
#CNTY    Country of publication
#PAPX    MeSH pharmacological action pre-explosions
#GRNT    NIH Grant Numbers
#MDAT    Date of last modification
#CDAT    Date of completion
#PID     Publisher ID
#FAUT    First Author of publication
#FULL    Full Author Name(s) of publication
#FINV    Full name of investigator
#TT      Words in transliterated title of publication
#LAUT    Last Author of publication
#PPDT    Date of print publication
#EPDT    Date of Electronic publication
#LID     ELocation ID
#CRDT    Date publication first accessible through Entrez
#BOOK    ID of the book that contains the document
#ED      Section's Editor
#ISBN    ISBN
#PUBN    Publisher's name
#AUCL    Author Cluster ID
#EID     Extended PMID
#DSO     Additional text from the summary
#AUID    Author Identifier
#PS      Personal Name as Subject

bpy12. Using egquery EUtil in Biopython

We can use Bio.Entrez.egquery to find counts in different databases, for a query term.

For Entrez.email, you should give your email. The query term, here, is 'asthma'. Biopython returns a dictionary of length 2, with the result in key 'eGQueryResult". The result inside 'eGQueryResult' are also a dictionary.

The dictionary within the result is iterated over (with row being the current subkey). Usually, we are only interested in two subkeys; the database name 'MenuName' and the count 'Count'.

We populate lists X and Y with the strings from the 2 subkeys. Finally, we print the two lists, provided count is not zero.

# bpy12.py
from __future__ import print_function, division
from Bio import Entrez

Entrez.email = "Your.Name.Here@example.org"
handle = Entrez.egquery(term="asthma")
record = Entrez.read(handle)
handle.close()
X=[]
Y=[]
for row in record["eGQueryResult"]:
    for i in row:
        if i == 'MenuName':
           Y.append(row[i])
        elif i=='Count':
            X.append(row[i])
for i in range(len(X)):
    if int(X[i])==0:
        continue
    print("%20s\t%s" % (Y[i],X[i]))

#              PubMed    149142
#      PubMed Central    106753
#                MeSH    11
#               Books    7165
#       PubMed Health    1461
#                OMIM    270
#         Site Search    75
#          Nucleotide    32446
#                 GSS    91
#                 EST    210
#             Protein    2558
#              Genome    1
#           Structure    390
#               dbVar    9
#                Gene    1034
#                 SRA    8245
#          BioSystems    125
#             UniGene    1
#   Conserved Domains    22
#              PopSet    7
#        GEO Profiles    731699
#        GEO DataSets    4417
#          HomoloGene    18
#    PubChem Compound    122
#   PubChem Substance    2125
#    PubChem BioAssay    3515
#         NLM Catalog    2610
#               Probe    34
#               dbGaP    2387
#          BioProject    209
#           BioSample    12542

bpy11. Using Entrez EUtils in Biopython

We can use Bio.Entrez to use the different EUtils at NCBI Entrez.

There are 9 utilities, and currently 8 of them can be accessed using Bio.Entrez. The 9th is the most recent, ECitMatch. Other than using Biopython, you can also use HTML requests with appropriate query term, name of EUtil, etc. It can return data as XML, Python object, etc.

Before going into details of the different EUtils, this program just lists some documentation. You can get complete documentation by removing 2 from docL term, in which case it will print all the lines.

# bpy11.py

from __future__ import print_function, division
from Bio import Entrez

eUtils = [e for e in Entrez.__dict__ if e[0]=='e' and e != 'email']
for i,eUtil in enumerate(eUtils):
    print('%d.%s' % (i+1,eUtil))
    exec('doc=Entrez.%s.__doc__' % eUtil)
    docL = doc.split('\n')
    print(''.join(docL[:2]))
    print('')

#1.espell
#ESpell retrieves spelling suggestions, returned in a results handle.
#
#2.egquery
#EGQuery provides Entrez database counts for a global search.
#
#3.esummary
#ESummary retrieves document summaries as a results handle.
#
#4.efetch
#Fetches Entrez results which are returned as a handle.
#
#5.esearch
#ESearch runs an Entrez search and returns a handle to the results.
#
#6.elink
#ELink checks for linked external articles and returns a handle.
#
#7.einfo
#EInfo returns a summary of the Entez databases as a results handle.
#
#8.epost
#Post a file of identifiers for future use.

bpy10. Reading a 1-record Genbank file using Biopython

Search for 'GenBank: BC135714.1' at NCBI.

In the results, there is Literature, Health, Genomes, and other sections. In the Genomes section, select the Nucleotide. If there is only 1 nucleotide, it will go to the Genbank file with that Locus ID.

In the upper right, select down arrow for send. Make sure complete record is selected, and then choose destination of File. Download options will come, and download the Genbank file.

Rename the file to BC135714.1.gb and save it to the working directory or a subfolder, such as data, under the working directory.

In this program, the function Bio.SeqIO.read is used to parse the text file. We print some information along with translated protein from a feature. The selected feature is the gene for the protein (ppm1a).

# bpy10.py
from __future__ import print_function, division
from Bio import SeqIO

gb = SeqIO.read('data/BC135714.1.gb','genbank')
print('type(gb) =',type(gb), end='\n\n')
print('id =', gb.id, end='\n\n')
print('desc =',gb.description, end='\n\n')
for annot in gb.annotations:
    print("--> ",annot)
    if annot == 'source':
        print(gb.annotations[annot], end='\n\n')
for ft in gb.features:
    print('--ft loc-->',ft.location)
dna = ft.extract(gb.seq) # ft is last feature ('ppm1a protein')
print("3rd feature translation: ",dna.translate())

#type(gb) = <class 'Bio.SeqRecord.SeqRecord'>
#
#id = BC135714.1
#
#desc = Xenopus tropicalis protein phosphatase 1A (formerly 2C),
# magnesium-dependent, alpha isoform, mRNA (cDNA clone MGC:121664
# IMAGE:7625995), complete cds.
#
#-->  comment
#-->  sequence_version
#-->  source
#Xenopus (Silurana) tropicalis (western clawed frog)
#
#-->  taxonomy
#-->  keywords
#-->  references
#-->  accessions
#-->  data_file_division
#-->  date
#-->  organism
#-->  gi
#--ft loc--> [0:3385](+)
#--ft loc--> [0:3385](+)
#--ft loc--> [336:1488](+)
#3rd feature translation:  
#MGAFLDKPKMEKHNAHGQGNGLRYGLSSMQGWRVEMEDAHTAVIGLPNGLD
#AWSFFAVYDGHAGSQVAKYCCEHLLDHITSNQDFKGTDGHLSVWSVKNGIR
#TGFLQIDEHMRVISEKKHGADRSGSTAVGVMTSPNHIYFINCGDSRGLLCR
#SKKVHFFTQDHKPSNPLEKERIQNAGGSVMIQRVNGSLAVSRALGDFDYKC
#VHGKGPTEQLVSPEPEVYEIERSEEDDQFIILACDGIWDVMGNEELCDFVW
#SRLEVTDDLERVCNEIVDTCLYKGSRDNMSVILICFPSAPKVLPEAVKREA
#ELDKYLEGRVEDIIKKQGEEGVPDLVHVMRTLASENIPNLPPGGELASKRS
#VIEAVYNRLNPYRNDDTDSASTDDMW*

bpy9. Counting k-mer words in Python

The function count, below, will find the number of times a pattern appears.

For every position in text, minus length of pattern, it checks for equality with pattern, and returns the total number of equalities.

# bpy9.py

from __future__ import print_function, division
import numpy as np

def rand_dna(n):
    dna = ['A','C','T','G']
    seq = np.random.choice(dna,n)
    return "".join(seq)

def count(text, pattern):
    cnt = 0
    for i in range(len(text)-len(pattern)):
        if text[i:i+len(pattern)] == pattern: cnt+=1
    return cnt

np.random.seed(34)
dna = rand_dna(50)
dna_pattern = rand_dna(3)
c = count(dna, dna_pattern)
print("dna = ",dna)
print("dna_pattern = ", dna_pattern)
print("c =",c)

# dna =  CTTCTCAGCGTTCTCCACTCTACATCAGGCCGAAGTCAATTCATCTTTGT
# dna_pattern =  CTT
# c = 2

bpy8. Translating DNA using Biopython

There are 2 python files in this example. The first (gc.py) has only data, and the second one (bp8.py) has the actual code. We can load our data using:

from gc import *

This will import all data into the namespace. Also it will load functions, in case, there are any present.

The dictionary gc is the standard genetic code.

In the loop, we read the next 3 characters, and then use the dictionary to lookup what protein to append to the out list.

We can also use a list comprehension and make the program shorter. However, we can just use the translate method in Biopython.

# gc.py
gc = {'TTT':'F', 'TCT':'S', 'TAT':'Y', 'TGT':'C',
      'TTC':'F', 'TCC':'S', 'TAC':'Y', 'TGC':'C',
      'TTA':'L', 'TCA':'S', 'TAA':'*', 'TGA':'*',
      'TTG':'L', 'TCG':'S', 'TAG':'*', 'TGG':'W',
      'CTT':'L', 'CCT':'P', 'CAT':'H', 'CGT':'R',
      'CTC':'L', 'CCC':'P', 'CAC':'H', 'CGC':'R',
      'CTA':'L', 'CCA':'P', 'CAA':'Q', 'CGA':'R',
      'CTG':'L', 'CCG':'P', 'CAG':'Q', 'CGG':'R',
      'ATT':'I', 'ACT':'T', 'AAT':'N', 'AGT':'S',
      'ATC':'I', 'ACC':'T', 'AAC':'N', 'AGC':'S',
      'ATA':'I', 'ACA':'T', 'AAA':'K', 'AGA':'R',
      'ATG':'M', 'ACG':'T', 'AAG':'K', 'AGG':'R',
      'GTT':'V', 'GCT':'A', 'GAT':'D', 'GGT':'G',
      'GTC':'V', 'GCC':'A', 'GAC':'D', 'GGC':'G',
      'GTA':'V', 'GCA':'A', 'GAA':'E', 'GGA':'G',
      'GTG':'V', 'GCG':'A', 'GAG':'E', 'GGG':'G' }

# bpy8.py

from __future__ import print_function, division
from Bio.Seq import Seq
import numpy as np
from gc import *

def rand_dna(n):
    dna = ['A','C','T','G']
    seq = np.random.choice(dna,n)
    return "".join(seq)

data = rand_dna(45)
out = []
L = len(data)
assert(L%3==0)
for i in range(L//3):
    in_data = data[3*i:3*(i+1)]
    out.append(gc[in_data])

print("DNA: %s" % data)
print("According to gc, the protein: %s" % ("".join(out)))

dna = Seq(data)
prot = dna.translate()
print("According to Bio, the protein: %s" % prot)

# DNA: ATCTATTACTCTGCGGATATCTCAAAGGGCCGGGTGCGGCCCTGA
# According to gc, the protein: IYYSADISKGRVRP*
# According to Bio, the protein: IYYSADISKGRVRP*