Bio informatics: Evolution
Session 1: Database querying with BLAST, multiple sequence alignment
and phylogenetic reconstruction
Jupyter notebook: exerciseYFV.blast.search.AliView
- AliView: not on the exam
Goal of this session:
- Query GenBank database via web interface and through command line
- Look up sequences similar to a query in a database with BLAST via the
command line
- create and edit a multiple sequence alignment
- reconstruct a phylogeny with the NJ algorithm
Workflow:
- Look up patients virus sequence => YFV sequence
- Fetch a set of highly similar sequences
- Create and edit multiple sequence alignment
- Reconstruct NJ tree
The goal is to investigate the genetic relationship of a patient’s genome with publicly
available sequences (focusing on samples from 1 outbreak), to determine whether
the patient’s virus strain is more closely related to strains from outbreak or if it
represents a district lineage from a different region.
1) Load Python modules to be used
= from Bio.Blast import NCBIWWW, NCBIXML
from Bio import Entrez, SeqIO, Phylo, AlignIO
import pandas as pd
import time, os, alv, re
from IPython.display import Image
*Change nothing*
2) Load interface to run R embedded in a Python process:
= %load_ext rpy2.ipython
3) Load R packages to be used:
= %%R
,library("ape")
4) Define variables for input names that will only will have to be changed once
= baseName = "namefile(without .fasta)"
Base name is main part of the filename without details
= folderName = "exercise_xxx"
Name of directory where these files are stored
5) Push that variable to R using Rpush
= %Rpush baseName folderName
6) Check if this works:
= %%R
print(paste0("baseName: ", baseName, " -- folderName: ", folderName))
7) Let’s check this for Bash:
= %%bash -s "$baseName" "$folderName"
echo baseName: "$1" -- folderName: "$2"
8) Create a folder where analysis files can be saved:
= %%bash -s "$folderName"
= mkdir "$1"
If folder already exists, try $mkdir -p "$1
Database querying
First step = look up patient’s sequence in GenBank with accession number MF347613
=> download it => browse to https://www.ncbi.nlm.nih.gov/nuccore and save file as
fasta-format => upload in Jupyter
(1) Go to GenBank: paste accession number in search bar
(2) Download as fasta & change name
(3) Upload in Jupyter
,1) Move uploaded file into the folder for this exercise or upload sequence file to the
exercise folder (skip this command then)
= %%bash -s "$baseName" "$folderName"
mv ${1}.fasta ${2}/${1}.fasta
2) Change directory to exercise’s folder
= os.chdir(folderName)
3) Use query sequence to fetch the 100 most similar sequence available in NCBI
nucleotide database:
= inputFile = baseName+'.fasta' #your sequence file
query = SeqIO.read(inputFile, format="fasta") #Read sequence
search_type = "blastn"
database = "nt"
max_nr_hits = 100 #Adjust if you want more/less results
before = round(time.time(), ndigits=0)
blast_stream = NCBIWWW.qblast(program = search_type,
database = database,
sequence = query.seq,
alignments = 1,
hitlist_size = max_nr_hits)
after = round(time.time(), ndigits=0)
text = f"The BLAST search took {after - before} seconds."
print(text)
This script performs BLAST search using NCBIWWW.qblast function from
Biopython, which queries the NCBI BLAST database
Defines input file
Reads query sequence: uses SeqIO.read() from Biopython to read a FASTA file
& stores the sequence object in query
Uses BLASTN for nucleotide sequence comparison
Searches against NCBI’s nucleotide database
Limits results to 100 matches
Stores current time before BLAST search starts
NCBIWWW.qblast() submits the query sequence (query.seq) to NCBI’s online
BLAST service
, alignments = 1 → Retrieves only the top alignment per hit
hitlist_size = max_nr_hits → Limits the result list to 100 hits
captures the time after the BLAST search finishes
computes time taken for the BLAST search
In Python, a stream refers to a sequence of data elements, often used for input and
output operations. In this case, blast_stream holds the result of our search
4) Save the result to file:
= blastn_result_file = baseName+".blastn.result.xml"
with open(blastn_result_file, 'w') as file:
file.write(blast_stream.read())
Saves BLAST search results in a XML file
5) Extract the sequences in fasta format from the BLAST search results
=all_data_fasta = baseName+".blastn.result.fasta"
with open(blastn_result_file) as result_handle:
blast_records = NCBIXML.parse(result_handle)
with open(all_data_fasta, 'w') as fasta_handle:
for blast_record in blast_records:
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
db_seq = hsp.sbjct
db_id = alignment.accession
fasta_handle.write(f">{db_id}\n{db_seq}\n")
Extracts the BLAST search results from previously saved XML file and saves
the matched sequences in FASTA format
Multiple sequence alignment, visualization and trimming
1) [Bash]: Check whether required files are present
= %%bash -s "$all_data_fasta"
ls | grep ${1}
Should give you the .blastn.result.fasta file
2) [Bash]: perform alignment using mafft: mafft –auto inputFile.fas > outputFile.fas
Session 1: Database querying with BLAST, multiple sequence alignment
and phylogenetic reconstruction
Jupyter notebook: exerciseYFV.blast.search.AliView
- AliView: not on the exam
Goal of this session:
- Query GenBank database via web interface and through command line
- Look up sequences similar to a query in a database with BLAST via the
command line
- create and edit a multiple sequence alignment
- reconstruct a phylogeny with the NJ algorithm
Workflow:
- Look up patients virus sequence => YFV sequence
- Fetch a set of highly similar sequences
- Create and edit multiple sequence alignment
- Reconstruct NJ tree
The goal is to investigate the genetic relationship of a patient’s genome with publicly
available sequences (focusing on samples from 1 outbreak), to determine whether
the patient’s virus strain is more closely related to strains from outbreak or if it
represents a district lineage from a different region.
1) Load Python modules to be used
= from Bio.Blast import NCBIWWW, NCBIXML
from Bio import Entrez, SeqIO, Phylo, AlignIO
import pandas as pd
import time, os, alv, re
from IPython.display import Image
*Change nothing*
2) Load interface to run R embedded in a Python process:
= %load_ext rpy2.ipython
3) Load R packages to be used:
= %%R
,library("ape")
4) Define variables for input names that will only will have to be changed once
= baseName = "namefile(without .fasta)"
Base name is main part of the filename without details
= folderName = "exercise_xxx"
Name of directory where these files are stored
5) Push that variable to R using Rpush
= %Rpush baseName folderName
6) Check if this works:
= %%R
print(paste0("baseName: ", baseName, " -- folderName: ", folderName))
7) Let’s check this for Bash:
= %%bash -s "$baseName" "$folderName"
echo baseName: "$1" -- folderName: "$2"
8) Create a folder where analysis files can be saved:
= %%bash -s "$folderName"
= mkdir "$1"
If folder already exists, try $mkdir -p "$1
Database querying
First step = look up patient’s sequence in GenBank with accession number MF347613
=> download it => browse to https://www.ncbi.nlm.nih.gov/nuccore and save file as
fasta-format => upload in Jupyter
(1) Go to GenBank: paste accession number in search bar
(2) Download as fasta & change name
(3) Upload in Jupyter
,1) Move uploaded file into the folder for this exercise or upload sequence file to the
exercise folder (skip this command then)
= %%bash -s "$baseName" "$folderName"
mv ${1}.fasta ${2}/${1}.fasta
2) Change directory to exercise’s folder
= os.chdir(folderName)
3) Use query sequence to fetch the 100 most similar sequence available in NCBI
nucleotide database:
= inputFile = baseName+'.fasta' #your sequence file
query = SeqIO.read(inputFile, format="fasta") #Read sequence
search_type = "blastn"
database = "nt"
max_nr_hits = 100 #Adjust if you want more/less results
before = round(time.time(), ndigits=0)
blast_stream = NCBIWWW.qblast(program = search_type,
database = database,
sequence = query.seq,
alignments = 1,
hitlist_size = max_nr_hits)
after = round(time.time(), ndigits=0)
text = f"The BLAST search took {after - before} seconds."
print(text)
This script performs BLAST search using NCBIWWW.qblast function from
Biopython, which queries the NCBI BLAST database
Defines input file
Reads query sequence: uses SeqIO.read() from Biopython to read a FASTA file
& stores the sequence object in query
Uses BLASTN for nucleotide sequence comparison
Searches against NCBI’s nucleotide database
Limits results to 100 matches
Stores current time before BLAST search starts
NCBIWWW.qblast() submits the query sequence (query.seq) to NCBI’s online
BLAST service
, alignments = 1 → Retrieves only the top alignment per hit
hitlist_size = max_nr_hits → Limits the result list to 100 hits
captures the time after the BLAST search finishes
computes time taken for the BLAST search
In Python, a stream refers to a sequence of data elements, often used for input and
output operations. In this case, blast_stream holds the result of our search
4) Save the result to file:
= blastn_result_file = baseName+".blastn.result.xml"
with open(blastn_result_file, 'w') as file:
file.write(blast_stream.read())
Saves BLAST search results in a XML file
5) Extract the sequences in fasta format from the BLAST search results
=all_data_fasta = baseName+".blastn.result.fasta"
with open(blastn_result_file) as result_handle:
blast_records = NCBIXML.parse(result_handle)
with open(all_data_fasta, 'w') as fasta_handle:
for blast_record in blast_records:
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
db_seq = hsp.sbjct
db_id = alignment.accession
fasta_handle.write(f">{db_id}\n{db_seq}\n")
Extracts the BLAST search results from previously saved XML file and saves
the matched sequences in FASTA format
Multiple sequence alignment, visualization and trimming
1) [Bash]: Check whether required files are present
= %%bash -s "$all_data_fasta"
ls | grep ${1}
Should give you the .blastn.result.fasta file
2) [Bash]: perform alignment using mafft: mafft –auto inputFile.fas > outputFile.fas
