Python, Biopython and Data Libraries in Bioinformatics

Python is widely used in bioinformatics to analyze DNA sequences, protein sequences, gene expression data, and other biological information. Its simple syntax and extensive collection of libraries make it useful for processing biological data and automating research tasks.

Biopython provides specialized tools for handling biological sequences, sequence files, databases, alignments, and phylogenetic trees. Other Python libraries, such as NumPy and Pandas, help process numerical data and tables, while Matplotlib and Seaborn are used to visualize results.

In this article, we will learn how to use Python and its major libraries in bioinformatics, explore practical code examples, and implement the Needleman–Wunsch algorithm for global sequence alignment.

What Is Python in Bioinformatics?

Python in bioinformatics refers to using the Python programming language to process, analyze, and visualize biological data.

Researchers use Python to automate repetitive tasks, compare biological sequences, analyze gene expression, retrieve information from public databases, and study evolutionary relationships.

Common applications include:

• DNA and protein sequence analysis.
• Genome and transcriptome data processing.
• Sequence alignment and comparison.
• Gene expression analysis.
• Phylogenetic tree analysis.
• Biological data visualization.

Installing Python Libraries for Bioinformatics

Before running the examples, install the required libraries using the following command:

pip install biopython numpy pandas matplotlib seaborn

This command installs Biopython for biological data, NumPy for numerical calculations, Pandas for tabular data, and Matplotlib and Seaborn for visualization.

You can run the examples in a Python environment such as VS Code, Jupyter Notebook, or another environment that supports the required libraries.

1. Biopython: Python Library for Bioinformatics

Biopython is an open-source collection of Python tools designed for computational biology and bioinformatics. It simplifies common tasks such as manipulating sequences, reading biological file formats, retrieving database records, and analyzing alignments.

One of its basic components is the Seq object, which provides operations specifically designed for biological sequences.

Working with DNA and Protein Sequences

The following example demonstrates how to calculate the reverse complement, transcribe DNA into RNA, and translate DNA into a protein sequence.

from Bio.Seq import Seq

dna = Seq("ATGGCCATTGTAATGGGCCGCTGA")

print(dna.reverse_complement())
print(dna.transcribe())
print(dna.translate())
 

Output:

TCAGCGGCCCATTACAATGGCCAT AUGGCCAUUGUAAUGGGCCGCUGA MAIVMGR*

Explanation:

• reverse_complement() returns the reverse complement of the DNA sequence.
• transcribe() converts DNA into RNA by replacing thymine (T) with uracil (U).
• translate() converts codons into amino acids using the genetic code.

The * symbol represents a stop codon in the translated protein sequence.

2. SeqIO: Reading and Writing Biological Sequence Files

Biological sequences are commonly stored in formats such as FASTA and GenBank. Biopython's SeqIO module provides functions for reading, writing, and converting these files.

Reading a FASTA File

Suppose genes.fasta contains multiple DNA sequences. The following code reads each sequence and displays its identifier and length.

from Bio import SeqIO

for record in SeqIO.parse("genes.fasta", "fasta"):
   print(record.id, len(record.seq))
 

Explanation:

• SeqIO.parse() reads multiple sequence records.
• record.id provides the sequence identifier.
• len(record.seq) returns the sequence length.

Important SeqIO Functions

Converting Sequence Formats

The following example converts a GenBank file into FASTA format.

from Bio import SeqIO

SeqIO.convert(
   "plasmid.gb",
   "genbank",
   "plasmid.fasta",
   "fasta",
)

print("File conversion completed.")
 

The input file must exist and contain valid GenBank data.

Calculating GC Content

GC content represents the percentage of guanine (G) and cytosine (C) bases in a DNA sequence. It is useful for studying genome composition and comparing sequences.

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

for record in SeqIO.parse("genes.fasta", "fasta"):
   gc_percent = gc_fraction(record.seq) * 100
   print(record.id, round(gc_percent, 2), "%")
 

The gc_fraction() function calculates the fraction of G and C bases. Multiplying the result by 100 converts it into a percentage.

3. Entrez: Accessing NCBI Databases

Entrez is a system for searching and retrieving biological information from databases maintained by the National Center for Biotechnology Information (NCBI).

Biopython's Bio.Entrez module allows Python programs to access these databases without manually searching through web pages.

For example, researchers can search nucleotide records, retrieve DNA sequences, and use the results in their analysis.

Searching and Downloading DNA Sequences

The following example searches the NCBI nucleotide database for human BRCA1 records and retrieves one matching sequence.

from Bio import Entrez, SeqIO

Entrez.email = "your.name@example.com"

handle = Entrez.esearch(
   db="nucleotide",
   term="BRCA1[Gene] AND Homo sapiens[Organism]",
   retmax=3,
)

result = Entrez.read(handle)
handle.close()

ids = result["IdList"]

if ids:
   handle = Entrez.efetch(
       db="nucleotide",
       id=ids[0],
       rettype="fasta",
       retmode="text",
   )

   record = SeqIO.read(handle, "fasta")
   handle.close()

   print("Sequence ID:", record.id)
   print("Sequence length:", len(record.seq))
else:
   print("No matching records found.")
 

Explanation:

1. Entrez.email identifies the contact email associated with the requests.
2. Entrez.esearch() searches the selected NCBI database.
3. Entrez.read() reads the search results.
4. Entrez.efetch() retrieves a sequence record.
5. SeqIO.read() parses the downloaded FASTA record.

4. AlignIO: Working with Multiple Sequence Alignments

Multiple sequence alignment arranges three or more biological sequences so that corresponding positions appear in the same columns. It helps researchers identify conserved regions and sequence differences.

Biopython's AlignIO module reads and writes alignment files. Alignment programs such as MUSCLE and Clustal Omega can generate the alignments that are then processed using AlignIO.

Reading an Alignment File

Suppose aligned.fasta contains the following aligned sequences:

>Human
ATGGTC-CTA
>Chimp
ATGGTC-CTA
>Mouse
ATGATCACTA
 

The sequences must have equal lengths for this alignment format.

from Bio import AlignIO

alignment = AlignIO.read("aligned.fasta", "fasta")

print("Sequences:", len(alignment))
print("Alignment length:", alignment.get_alignment_length())
print("Column 4:", alignment[:, 3])
 

Output:

Sequences: 3

Alignment length: 10

Column 4: GGA

The expression alignment[:, 3] extracts the fourth column because Python indexing starts at zero.

Counting Conserved Columns

A conserved column contains the same character in every sequence.

from Bio import AlignIO

alignment = AlignIO.read("aligned.fasta", "fasta")

conserved = 0

for i in range(alignment.get_alignment_length()):
   column = alignment[:, i]

   if len(set(column)) == 1:
       conserved += 1

print("Conserved columns:", conserved)
 

Output:

Conserved columns: 8

5. Phylo: Analyzing Phylogenetic Trees

A phylogenetic tree represents evolutionary relationships among organisms, genes, or other biological sequences. It shows how the analyzed groups are connected according to the tree and its underlying data. Biopython's Bio.Phylo module supports reading, displaying, and analyzing phylogenetic trees in formats such as Newick.

Reading and Displaying a Phylogenetic Tree

from io import StringIO

from Bio import Phylo

newick = "((Human:0.1,Chimp:0.1):0.2,(Mouse:0.3,Rat:0.3):0.1);"

tree = Phylo.read(StringIO(newick), "newick")

Phylo.draw_ascii(tree)

print("Number of terminal taxa:", tree.count_terminals())
print("Terminal names:", [leaf.name for leaf in tree.get_terminals()])
print("Human to Mouse distance:", tree.distance("Human", "Mouse"))
 

Explanation:

• Phylo.read() loads the tree from Newick data.
• Phylo.draw_ascii() displays a text-based tree.
• count_terminals() counts the terminal taxa.
• get_terminals() returns the terminal nodes.
• distance() calculates the path length between two taxa using the supplied branch lengths.

Building a Tree from a Sequence Alignment

Biopython can also construct a basic distance-based tree from a multiple sequence alignment.

from Bio import AlignIO, Phylo
from Bio.Phylo.TreeConstruction import (
   DistanceCalculator,
   DistanceTreeConstructor,
)

alignment = AlignIO.read("aligned.fasta", "fasta")

calculator = DistanceCalculator("identity")
distances = calculator.get_distance(alignment)

constructor = DistanceTreeConstructor()
tree = constructor.nj(distances)

Phylo.draw_ascii(tree)
 

6. Pairwise Alignment in Biopython

Pairwise sequence alignment compares two DNA, RNA, or protein sequences to identify matching positions and differences.

Two common approaches are:

• Global alignment: Aligns sequences across their full lengths.
• Local alignment: Identifies the highest-scoring matching regions within sequences.

Alignment algorithms assign scores to matches, mismatches, and gaps. The algorithm then selects an alignment with the best score under the chosen scoring rules.

Using PairwiseAligner

Bio.Align.PairwiseAligner is the recommended modern interface for pairwise alignment in Biopython.
 

from Bio import Align
aligner = Align.PairwiseAligner()

aligner.mode = "global"
aligner.match_score = 1
aligner.mismatch_score = -1
aligner.open_gap_score = -2
aligner.extend_gap_score = -2

seq1 = "ACCGGT"
seq2 = "ACGT"

alignments = aligner.align(seq1, seq2)

print("Best score:", alignments.score)
print(alignments[0])
 

Explanation:

• mode selects global or local alignment.
• match_score specifies the score for matching characters.
• mismatch_score specifies the score for different characters.
• open_gap_score and extend_gap_score specify the penalties for opening and extending gaps.

With these settings, the best global alignment score is 0. Multiple alignments may share the same score.

Aligning Protein Sequences with BLOSUM62

Protein alignments often use substitution matrices because different amino acid substitutions have different biological significance.

BLOSUM62 is a commonly used protein substitution matrix.

from Bio import Align
from Bio.Align import substitution_matrices

aligner = Align.PairwiseAligner()

aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
aligner.open_gap_score = -10
aligner.extend_gap_score = -0.5

score = aligner.score("HEAGAWGHEE", "PAWHEAE")

print("Alignment score:", score)
 

This example calculates an alignment score using BLOSUM62 and the specified gap penalties.

What Is pairwise2?

Bio.pairwise2 is an older Biopython module that appears in many existing tutorials and scripts. New code should generally use Bio.Align.PairwiseAligner, while older code may still need to be understood or maintained.

7. NumPy: Numerical Computing in Bioinformatics

NumPy is a Python library for numerical calculations. It provides arrays and mathematical functions that help researchers process large sets of biological measurements efficiently.

In bioinformatics, NumPy is useful for working with gene expression values, read counts, statistical measurements, and numerical matrices.

Calculating Gene Expression Changes

Consider the following gene expression measurements for control and treated samples.

import numpy as np

control = np.array([10, 25, 8, 40])
treated = np.array([30, 24, 2, 41])

fold_change = np.log2(treated / control)

print(fold_change.round(2))
 

Output:

[ 1.58 -0.06 -2.    0.04]

Explanation:

• np.array() creates numerical arrays.
• np.log2() calculates the base-2 logarithm.
• round(2) displays the results to two decimal places.

A log2 fold change of approximately 1.58 represents a threefold ratio between the treated and control values. A value of -2 represents a treated-to-control ratio of 0.25.

Creating a Numerical Matrix

NumPy also supports multidimensional arrays, which are useful in scoring systems and matrix-based algorithms.

import numpy as np
matrix = np.zeros((4, 5), dtype=int)
matrix[0, :] = np.arange(5) * -2

print(matrix.shape)
print(matrix[0])
 

Output:

(4, 5) [ 0 -2 -4 -6 -8]

The example creates a matrix with four rows and five columns, then fills the first row with values from 0 to -8.

8. Pandas: Analyzing Biological Data Tables

Pandas is a Python library used to organize, filter, summarize, and analyze tabular data. Its main data structure, the DataFrame, stores information in rows and named columns. Researchers use Pandas to analyze gene lists, expression tables, sequence statistics, and annotation results.

Creating a Sequence Information Table

The following example extracts sequence identifiers, lengths, and GC percentages from a FASTA file.

import pandas as pd

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

rows = []

for record in SeqIO.parse("genes.fasta", "fasta"):
   rows.append({
       "id": record.id,
       "length": len(record.seq),
       "gc_percent": round(
           100 * gc_fraction(record.seq),
           2,
       ),
   })

df = pd.DataFrame(rows)

print(df.head())
print(df.describe())

df.to_csv("genes.csv", index=False)
 

Explanation:

• pd.DataFrame() converts the collected records into a table.
• df.head() displays the first five rows.
• df.describe() summarizes numerical columns.
• df.to_csv() saves the results as a CSV file.

You can also filter sequences longer than 1,000 bases:

print(df[df["length"] > 1000])

Analyzing Gene Expression Data

Pandas and NumPy can be combined to calculate expression changes and categorize genes.

import numpy as np
import pandas as pd

df = pd.DataFrame({
   "gene": ["geneA", "geneB", "geneC", "geneD"],
   "control": [10, 25, 8, 40],
   "treated": [30, 24, 2, 41],
})

df["log2_fc"] = np.log2(df["treated"] / df["control"])

df["change"] = np.where(
   df["log2_fc"] > 1,
   "up",
   np.where(
       df["log2_fc"] < -1,
       "down",
       "same",
   ),
)

print(df)
 

9. Matplotlib: Visualizing Biological Data

Matplotlib is a Python library for creating graphs and charts. It helps researchers present numerical results in a visual format.

Common applications include plotting gene expression values, sequence lengths, GC content, and experimental measurements.

Creating a Gene Expression Bar Chart

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

df = pd.DataFrame({
   "gene": ["geneA", "geneB", "geneC", "geneD"],
   "control": [10, 25, 8, 40],
   "treated": [30, 24, 2, 41],
})

df["log2_fc"] = np.log2(df["treated"] / df["control"])

fig, ax = plt.subplots(figsize=(6, 4))

ax.bar(df["gene"], df["log2_fc"])
ax.axhline(0, linewidth=0.8)

ax.set_xlabel("Gene")
ax.set_ylabel("Log2 fold change")
ax.set_title("Gene Expression Changes")

fig.tight_layout()
fig.savefig("fold_change.png", dpi=300)

plt.show()
 

The chart displays the log2 fold change for each gene. Positive values indicate higher treated-to-control ratios, while negative values indicate lower ratios.

Creating a Histogram of Sequence Lengths

import matplotlib.pyplot as plt

lengths = [
   320, 450, 510, 780, 910,
   1200, 1500, 640, 720, 880,
]

plt.hist(
   lengths,
   bins=5,
   edgecolor="black",
)

plt.xlabel("Sequence length (bp)")
plt.ylabel("Number of sequences")
plt.title("Distribution of Sequence Lengths")

plt.show()
 

A histogram shows how sequence lengths are distributed across different ranges.

For clear biological graphs, include descriptive titles, labeled axes, and appropriate units. High-resolution images can be saved using dpi=300 for reports and presentations.

10. Seaborn: Statistical Data Visualization

Seaborn is a statistical visualization library built on Matplotlib. It provides convenient functions for creating plots from Pandas DataFrames.

It is commonly used to visualize distributions, compare experimental groups, and display gene expression patterns.

Comparing Gene Expression Between Conditions

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

df = pd.DataFrame({
   "gene": ["geneA", "geneB", "geneC", "geneD"],
   "control": [10, 25, 8, 40],
   "treated": [30, 24, 2, 41],
})

long_df = df.melt(
   id_vars="gene",
   value_vars=["control", "treated"],
   var_name="condition",
   value_name="expression",
)

sns.barplot(
   data=long_df,
   x="gene",
   y="expression",
   hue="condition",
)

plt.title("Gene Expression by Condition")
plt.show()
 

The melt() function converts the expression table into a long format, making it easier to compare control and treated values.

Common Seaborn Plots in Bioinformatics

Creating a Gene Expression Heatmap

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

df = pd.DataFrame({
   "gene": ["geneA", "geneB", "geneC", "geneD"],
   "control": [10, 25, 8, 40],
   "treated": [30, 24, 2, 41],
})

expression = df.set_index("gene")

sns.heatmap(
   expression,
   annot=True,
   cmap="viridis",
)

plt.title("Gene Expression Heatmap")
plt.show()
 

The heatmap uses color intensity to represent expression values. For comparisons across genes, researchers may need to normalize or scale the values depending on the analysis objective.

Matplotlib vs. Seaborn: Matplotlib offers detailed control over graph elements, while Seaborn simplifies the creation of statistical visualizations.

11. Implementing the Needleman–Wunsch Algorithm in Python

The Needleman–Wunsch algorithm is a dynamic programming method used for global sequence alignment. It compares two sequences across their full lengths by finding an alignment with the highest score under a specified scoring system.

Implementing this algorithm helps explain how alignment programs evaluate matches, mismatches, and gaps.

How Does the Needleman–Wunsch Algorithm Work?

The algorithm follows three main steps:

1. Initialize the score matrix: Create a matrix and fill the first row and column with cumulative gap penalties.
2. Fill the matrix: Calculate each cell using the best score from a diagonal match or mismatch, a gap in one sequence, or a gap in the other sequence.
3. Trace back: Start at the bottom-right cell and follow the scoring decisions to construct the alignment.

Python Implementation

import numpy as np

def needleman_wunsch(seq1, seq2, match=1, mismatch=-1, gap=-2):
   n, m = len(seq1), len(seq2)

   score = np.zeros((n + 1, m + 1), dtype=int)

   score[:, 0] = np.arange(n + 1) * gap
   score[0, :] = np.arange(m + 1) * gap

   for i in range(1, n + 1):
       for j in range(1, m + 1):
           s = match if seq1[i - 1] == seq2[j - 1] else mismatch

           score[i, j] = max(
               score[i - 1, j - 1] + s,
               score[i - 1, j] + gap,
               score[i, j - 1] + gap,
           )

   a1, a2 = [], []
   i, j = n, m

   while i > 0 or j > 0:
       if i > 0 and j > 0:
           s = match if seq1[i - 1] == seq2[j - 1] else mismatch

           if score[i, j] == score[i - 1, j - 1] + s:
               a1.append(seq1[i - 1])
               a2.append(seq2[j - 1])
               i -= 1
               j -= 1
               continue

       if i > 0 and score[i, j] == score[i - 1, j] + gap:
           a1.append(seq1[i - 1])
           a2.append("-")
           i -= 1
       else:
           a1.append("-")
           a2.append(seq2[j - 1])
           j -= 1

   return (
       "".join(reversed(a1)),
       "".join(reversed(a2)),
       int(score[n, m]),
   )

first, second, best = needleman_wunsch("ACCGGT", "ACGT")

print(first)
print(second)
print("Best score:", best)
 

Output:

ACCGGT A-C-GT Best score: 0

Understanding the Result

The example uses the following scoring rules:

• Match: +1
• Mismatch: -1
• Gap: -2

The displayed alignment contains four matches and two gaps. Its score is therefore:

[
4(1) + 2(-2) = 0
]

Other alignments may produce the same optimal score. This implementation returns one optimal alignment and uses a simple linear gap penalty.

Conclusion

Python provides a flexible environment for analyzing biological data. Biopython simplifies sequence processing, database access, alignment analysis, and phylogenetic tree handling. NumPy and Pandas support numerical calculations and data organization, while Matplotlib and Seaborn help visualize biological results. The Needleman–Wunsch algorithm demonstrates how dynamic programming can identify an optimal global sequence alignment using a scoring matrix and traceback.