Python for Bioinformatics (Video)

0:03 Introduction

Now, with the advent of AI, it's very common for us to think about learning languages as not that important. People may say, why learn Python? Why learn R or any other language for bioinformatics? And my answer to that is that it really is important until the models get really, really good. It's very important to learn at least fundamentals of a programming language to have a kickstart in bioinformatics. I think that is very important.

0:34 Why learn Python

First, we will look at why to learn Python. Computational biology and bioinformatics both depend on programs that researchers write. This video builds a first program in Python. Python is a high-level, general-purpose programming language. Guido van Rossum developed the language at the CWI research institute in the Netherlands. In February 1991, after just over a year of work, he posted the first version to USENET, a network of online discussion groups. Python's own FAQ recommends Python as a first language for people who have never programmed, and it gives the simple and consistent syntax as the reason. In this script, each line holds one statement, and Python runs the lines in order from top to bottom. Because the language is general-purpose, a Python program can read a sequence file, run a statistical test, draw a figure or run an analysis pipeline. A script also records the exact steps of an analysis. You or a colleague can then run the same steps again on new data and compare the results. I request that you watch the video till the end because then we can make more of these videos as the algorithm recommends it to more people because of increased retention time. And that really helps us a lot.

1:51 Python in bioinformatics goes beyond sequences

Second, we will see that Python in bioinformatics goes beyond sequences. Python does those jobs through libraries. A library is a collection of code that other people wrote and shared. You load a library into your own script with an import statement, a line at the top of the script that names the library. The map sorts these tools into groups by the kind of job they do. Each tool is either written in Python or used from Python. Biopython started as a project in 1999, and this video uses Biopython for its sequence examples. Biopython reads and writes sequence files and alignments, which line up related sequences position by position. It also works with protein structures and with online databases. Biological data also arrive as numbers and tables. NumPy is the main Python library for arrays, which hold many numbers in a list or a grid. Pandas handles tables with labeled rows and columns, such as gene expression counts. Matplotlib turns such numbers into figures, including interactive figures that you can zoom and pan. Seaborn builds statistical graphics, such as box plots, on top of Matplotlib. SciPy provides algorithms for scientific computing, and its first release came in 2001. Scikit-learn is a library for machine learning, in which a program learns patterns from data. It supports supervised learning, which learns from examples with known labels, and unsupervised learning, which finds groups in data without labels. A Jupyter notebook keeps code and its results in one file, with explanations in plain language between the pieces of code. The notebook divides its content into blocks, the cells, and you run the code cells one at a time. Several Python tools work directly on sequencing data. Sequencing produces reads, short stretches of DNA sequence that FASTQ files store. Pysam works with FASTQ files and with the SAM and BAM formats, which record where each read matches a reference genome, the standard genome sequence of a species. Pysam also handles VCF files, which list the positions where a sample differs from the reference. Cutadapt removes adapter sequences, short synthetic pieces of DNA added to each fragment before sequencing, and other unwanted sequence such as primers from the reads. MultiQC gathers the quality checks that many tools run on many samples into one report. Cutadapt and MultiQC install with pip, the package installer for Python. Single-cell RNA sequencing measures gene expression in individual cells, and Scanpy is a Python toolkit for these data. Scanpy prepares the data and groups cells with similar expression into clusters. It also tests which genes differ in expression between clusters, and its Python implementation handles data sets of more than one million cells. When an analysis has many steps, Snakemake describes the workflow, the ordered chain of steps, in a Python-based language. The same workflow runs on a laptop or on a cluster of linked computers without changes. Python tools also work outside sequencing. A molecular dynamics simulation calculates how the atoms of a molecule move over time and saves their positions as a trajectory. MDAnalysis is a Python library that analyzes these trajectories. RDKit is a toolkit for cheminformatics, the computer analysis of chemical compounds such as drug molecules. Its core is written in C++, and you use RDKit from Python. For protein structure, ColabFold predicts the three-dimensional shape of a protein from its amino acid sequence. ColabFold combines MMseqs2, a program that searches databases for related sequences, with the prediction program AlphaFold2. You can run ColabFold free of charge in Google Colaboratory, a notebook service in the web browser. A beginner does not need these tools at the start. To use one of them on your own computer, you first install Python and then install the tool with pip.

5:54 Installing Python

Third, we will install Python. Installing Python comes first, because each tool on the map needs Python to run. You download Python from python.org, the official Python website. On Windows, python.org now recommends the Python install manager, which you install from python.org/downloads or from the Microsoft Store. After you install it, the command python works in any terminal, the window where you type commands as text. The older full installer is still available for Python 3.15, but the Python team is phasing it out and will not produce it for Python 3.16 or later. On Windows the older route needs one setting. PATH is the list of folders that the terminal searches when you type the name of a command. If you use the full installer, tick the box Add python.exe to PATH on its first page before you click Install Now. Without that tick, the terminal does not find the python command. On macOS, you download the installer package from python.org and double-click the downloaded file. When the installer finishes, you open the Python 3.15 folder in Applications and double-click Install Certificates.command to complete the installation. If you already use Homebrew, a package manager that adds programs to a Mac from the terminal, the command brew install python3 installs Python without the python.org package. On Windows you open Terminal or PowerShell from the Start menu, and on a Mac you open Terminal from the Utilities folder inside Applications. You check the installation by typing python --version on Windows or python3 --version on macOS. Python prints its version number, for example Python 3.15.0. The terminal then returns to the prompt, the line where you type the next command.

7:42 A first Python script

The fourth part writes a first Python script. With Python installed, you need a place to write programs. Visual Studio Code, or VS Code, is a code editor from Microsoft for Windows and macOS. You then install the Python extension from the Extensions panel inside VS Code. An extension is an add-on, and this one adds Python support to the editor. A Python file is a plain text file that holds Python statements, and its name ends in dot p y. Any text editor can open such a file. In VS Code you create a new file and name the file hello.py. The ending dot p y tells VS Code to treat the file as a Python program. You type one line that uses the print function. A function is a named block of code that does one job, and print writes a message to the terminal. In this script the message is Hello DNA, between quotation marks, and you save the file when the line is complete. To run the file, you open a terminal in the same folder and type python3 hello.py on macOS or python hello.py on Windows. Python reads the file and executes its statements from top to bottom. VS Code also has a Run Python File button in the top right of the editor, and the button runs the same command in its terminal panel. If you misspell print as p r n t, Python stops with a NameError, the error for a name that Python does not recognize. The message repeats the misspelled name and gives the number of the line that failed, and you correct that line and run the script again.

9:15 Installing Biopython with pip

Fifth, we will install Biopython with pip. The sequence examples in this video need Biopython, which is not part of the standard library, the large collection of code that comes with Python. You add Biopython and other libraries with pip, which downloads packages from the Python Package Index, or PyPI. Python from python.org or Homebrew already includes pip. You install Biopython with the command python3 -m pip install biopython on macOS, or py -m pip install biopython on Windows. The interpreter is the program that runs Python code, and the dash m form runs pip for the exact interpreter that you name. On a computer with more than one copy of Python, that choice decides which copy receives Biopython. To confirm that Biopython works, you import the Bio package in Python and print its version number. The same command, with a different package name, installs pandas, Matplotlib and the other libraries from the map.

10:19 A sequence is a string

Sixth, we will treat a sequence as a string. The second script works on a short DNA sequence, and its first steps use plain Python without Biopython. A DNA sequence is text, so Python stores the sequence as a string, an ordered series of characters. You write the name seq, an equals sign, and the 12 bases ATGGCTTGGTAA between quotation marks. This line makes seq a variable, a name that refers to a stored value. When you apply the len function to seq, Python returns 12, the number of bases. The count method, a function that you attach to the string with a dot, returns 4 when you give it the base G. Python numbers the characters of a string from 0. So seq with 0 in square brackets gives the first base, A, and seq with 1 gives the second base, T. A codon is a group of three bases that the ribosome reads together. A slice takes the characters between two positions. The slice from 0 to 3 gives ATG, the first codon, because a slice includes its start position and leaves out its end position. A script displays nothing unless you ask for output, so you wrap each of these steps in print.

11:37 GC content

Seventh, we will calculate GC content. Counting single bases leads to the first measurement in the script. GC content is the percentage of bases in a sequence that are G or C. GC content differs widely between organisms. In bacteria, genomic GC content ranges from below 20 percent to above 70 percent. In Python you add the count of G to the count of C, divide the total by the length, and multiply by 100. For this sequence, 4 G and 1 C make 5 bases. Five divided by twelve is 0.4167, and multiplying by 100 gives 41.7 percent. Python prints the unrounded result. To keep one decimal place, you pass the result and the number 1 to the round function, which gives 41.7.

12:27 Reverse complement and translation

The eighth part covers the reverse complement and translation. Some sequence operations take several lines with plain strings, and Biopython provides them as methods. For example, the opposite strand of a DNA molecule is the reverse complement of a sequence. You swap A with T and G with C, then read the result backwards. Biopython keeps a sequence in a Seq object, a value that comes with sequence methods. You import Seq from the Bio.Seq module, one file of code inside Biopython, and you pass your string to Seq. The reverse_complement method of the new object returns TTACCAAGCCAT. The translate method reads the sequence codon by codon with the standard genetic code, the table that links each codon to an amino acid. For this sequence it returns M, A and W, the one-letter codes for methionine, alanine and tryptophan, followed by an asterisk. A stop codon ends translation, and Biopython marks a stop codon with an asterisk by default. In this sequence TAA is the stop codon.

13:36 Reading a FASTA file

The last part reads a FASTA file. Real sequences arrive in files, so the last step reads a file instead of typing a sequence into the script. You keep the file in the same folder as your script, and you open a terminal in that folder. FASTA is a plain text format for sequences. A definition line starts with a greater-than sign followed by a sequence identifier, and the lines after the definition line hold the sequence. Biopython reads a FASTA file with the parse function of its SeqIO module. You give the function the file name and the format name fasta, and it returns one record for each sequence in the file. Each record holds the identifier from the definition line in its id attribute and the sequence in its seq attribute. A for loop repeats a block of code once for each item, here once for each record, and prints the identifier and the length of each sequence. You run the script with python3 read_fasta.py, in the same way as hello.py. Because the loop runs once per record, the same three lines read a file of three sequences or a file of ten thousand. Here is what you learned in this video. You installed Python from python.org and ran a first script, hello.py, from the terminal. You then installed Biopython with pip and worked with a short DNA sequence, first as a string and then as a Biopython Seq object. Finally, a for loop with the parse function of SeqIO read the records in a FASTA file and printed the identifier and length of each sequence. Here is a question for you to test what you learned. What does the slice from 3 to 6 return for the sequence ATGGCTTGGTAA? Pause here and try it before the answer appears. Python counts positions from zero and stops before the end index, so the slice from 3 to 6 returns positions 3, 4 and 5, which hold G C T.