Practical - Similary searches

Search for homologous proteins using BLASTP (NCBI BLAST)

We want to find homologs of a bat protein sequence available on moodle.

Composition

Exercise 1 What can you say about the composition of this protein sequence? What problem can be encountered in a similarity search with such a composition?

Databases

Exercise 2 Which databases are available for the search?

Select the SwissProt database.

Exercise 3 Will the search be exhaustive?

BLASTP parameters

See Algorithm parameters, at the bottom of the page.

Exercise 4 How many sequences can you get?

Exercise 5 What is the Expect threshold to display sequences? Modify the Expect threshold to obtain results with E≤10.

Exercise 6 What are the gap penalties? The word size?

Exercise 7 How are regions of low complexity processed?

Perform the similarity search (SwissProt, E≤10, default for other parameters)

Graphic summary

Exercise 8 How many sequences do you obtain with E≤10?

Exercise 9 What can be deduced from the graphical representation of BLAST results? There is no alignment on the C-terminal region. Why?

Description

Exercise 10 Which families are detected?

Taxonomy

Exercise 11 In which organisms do you find hits?

Exercise 12 In which species is the closest sequence found?

Alignments

Exercise 13 Identify a conserved motif in the superfamily.

Exercise 14 According to the alignments, are the detected sequences all homologous to our query?

Decrease in the word size

Exercise 15 What might be the consequences of a decrease in word size?

Repeat the search with the parameters used in question 4 but with a word size of 3 aa.

Exercise 16 How many “hits” are you getting now?

Exercise 17 In which taxon (taxa) are similar proteins found?

Localisation of a transcript on a genome and search for a CDS

A gene has been detected as being over-expressed in several patients with esophageal cancer. The sequence of the transcript is provided below.

>human mRNA
GTGTGGACACTCCTAGGTTAGAAAGTTTGGTATGTTGCTATACCTTTGCTTCTCCCACCT
TCCCCAATATCTAATATGTATTTCTCATTCTTAGAATAATCCAGAATGGCTACTCTGATC
TATGTTGATAAGGAAAATGGAGAACCAGGCACCCGTGTGGTTGCTAAGGATGGGCTGAAG
CTGGGGTCTGGACCTTCAATCAAAGCCTTAGATGGGAGATCTCAAGTTTCAACACCACGT
TTTGGCAAAACGTTCGATGCCCCACCAGCCTTACCTAAAGCTACTAGAAAGGCTTTGGGA
ACTGTCAACAGAGCTACAGAAAAGTCTGTAAAGACCAAGGGACCCCTCAAACAAAAACAG
CCAAGCTTTTCTGCCAAAAAGATGACTGAGAAGACTGTTAAAGCAAAAAGCTCTGTTCCT
GCCTCAGATGATGCCTATCCAGAAATAGAAAAATTCTTTCCCTTCAATCCTCTAGACTTT
GAGAGTTTTGACCTGCCTGAAGAGCACCAGATTGCGCACCTCCCCTTGAGTGGAGTGCCT
CTCATGATCCTTGACGAGGAGAGAGAGCTTGAAAAGCTGTTTCAGCTGGGCCCCCCTTCA
CCTGTGAAGATGCCCTCTCCACCATGGGAATCCAATCTGTTGCAGTCTCCTTCAAGCATT
CTGTCGACCCTGGATGTTGAATTGCCACCTGTTTGCTGTGACATAGATATTTAAATTTCT
TAGTGCTTCAGAGTTTGTGTGTATTTGTATTAATAAAGCATTCTTTATCAGAAAAAAAAA
AAAAAAA

Localization of a mRNA sequence on a genome

Localize the corresponding gene on the human genome with BLAST (BLAST Genomes section, at the bottom of the page). Use the GRCh38.p14 reference assembly.

Program

MEGABLAST is proposed by default.

Exercise 18 What is the word size in MEGABLAST? Can we use MegaBLAST in our case or should we use BLASTN?

Analysis

Exercise 19 Comment obtained results (number of hits, conservation…)

Exercise 20 Give the localization of the transcribed region (chromosome, strand, start, end).

Exercise 21 How many exons are present in this gene? Can you give the precise boundaries of exons?

Exercise 22 How can you explain the ranking provided by BLAST?

Exercise 23 Which program can be used to improve your results? Compare results.

Localization of the coding sequence (CDS)

Exercise 24 How could we localize the coding region of this transcript?

Exercise 25 Perform the similarity search in SwissProt.

Exercise 26 What are the boundaries and reading frame of the CDS?

Gene identification in metagenomic data

This exercise is based on content from A Quick Guide to Getting Started with Shell and BLAST.

The objective is to identify the existence of nitrogen fixation genes (nifH, nifD, nifK) in soil metagenome. These genes encode enzymes that convert atmospheric nitrogen (N₂) into ammonia (NH₃). This process is essential for maintaining soil fertility and supporting plant growth. In this exercise, we will focus on identifying nifH genes.

Connection to the Server

To learn how to connect to the server for practical work, refer to the Server Connection Guide.

Once you’re connected, start by creating a working directory for this session and navigating into it.

mkdir

The mkdir command stands for “make directory” and is used to create a new folder in the file system. You can create your directory like this:

mkdir <directory_name>

cd

The cd command stands for “change directory”, allows you to move into a different folder within the file system. Use it as follows:

cd <directory>
mkdir tp3
cd tp3
Help

If you need help with a command, display its help using: - <command-h or <command--help for short help messages. - man <command> to display the full manual of a command.

Exploring data

The materials used in this are located in the directory /data/blast_tp/.

The objective is to identify the existence of nitrogen fixation genes* in soil metagenome.

ls

The ls command is used to list the contents of a directory in Linux. It provides a way to see files and subdirectories within the specified directory.

ls <directory>

tree

The treecommand is used to display the contents of a directory in a hierarchical structure. It provides a visual representation of files and subdirectories within the specified directory.

tree <directory>

Exercise 27 What are in /data/blast_tp/ ? Use lscommand. Explore the subdirectory with the tree command and map out the directory structure.

ls /data/blast_tp/
tree /data/blast_tp/

Navigate to blast_tp directory thanks to cd command.

cd /data/blast_tp/

This directory contains:

  • nifh_db/: a folder containing a nucleotide sequence file of 50 nifh genes (nifh-ref.fa), it’s the database.
  • query_gene/: a folder contains metagenome files from three different Iowa crop soils (corn.fa, soybean.fa, and prairie.fa), there are the queries.

head

The head command is used to display the first few lines of a file, which is helpful when you want to quickly preview its content without opening the entire file.

head <file>

By default, head shows the first 10 lines of the file. You can customize the number of lines displayed using the -n option. This option allows you to specify the exact number of lines to display.

head -n 20 <file>

This command will display the first 20 lines of the file. You can adjust the number of lines according to your needs.

Exercise 28 How does sequence fasta file look like? Use the command head to display the begining of each file.

You can display the contents of each file individually:

head metags/corn.fa 
head metags/prairie.fa 
head metags/soybean.fa
head nifh_db/nifh-ref.fa

Or use * to apply the head command to all files at once:

head */*.fa

The * is a wildcard that matches multiple files. In this case: */*.fa expands to all .fa files located inside any subdirectory.

This approach is useful when working with multiple files in different directories, avoiding the need to list them manually.

grep

grep, short for “global regular expression print”, is a command used for searching and matching text patterns in files contained in the regular expressions.

Display the lines with a pattern:

grep "<pattern>" <file_path>

This command searches through the specified file and outputs the lines that contain the given pattern.

Count the lines with a pattern:

grep -c "<pattern>" <file_path>

The -c option tells grep to count the number of lines that match the specified pattern rather than displaying them. This is useful for quickly determining how many occurrences of the pattern exist in the file.

Exercise 29 How many sequences are there in theses fasta files?

You can count for each file individually:

grep -c ">" metags/corn.fa 
grep -c ">" metags/prairie.fa 
grep -c ">" metags/soybean.fa
grep -c ">" nifh_db/nifh-ref.fa

Or use * to apply the grep command to all files at once:

grep -c ">" */*.fa

The * is a wildcard that matches multiple files. In this case: */*.fa expands to all .fa files located inside any subdirectory.

This approach is useful when working with multiple files in different directories, avoiding the need to list them manually.

Make database

We want to find the existence of nifh genes in soil metagenome samples. Making a database for nifh genes is a smart way to index the gene sequences so that it can search and compare metagenome sequences (query) with nifh gene sequences much faster.

makeblastdb

The makeblastdb command is used to create a BLAST database from your sequence data, which serves as the reference for comparing queries. Creating a database for genes is an efficient way to index gene sequences, allowing for much faster searches and comparisons against the database.

makeblastdb -in <sequence_file> -dbtype <type> -out <name>

Explanation of the command options:

  • -in <sequence_file>: Specifies the input file containing the sequences you want to include in the database. This file should be in FASTA or other supported formats.
  • -dbtype <type>: Indicates the type of sequences in the input file. Use nucl for nucleotide sequences and prot for protein sequences.
  • -out <name>: Defines the name for the output database files. The program will create several files with this base name, which will be used for subsequent BLAST searches.
makeblastdb -in <sequence_file> -dbtype <type> -out <name>

Navigate back to the directory of the current practical and construct the BLAST database.

cd ~/tp3/
makeblastdb -in /data/blast_tp/nifh_db/nifh-ref.fa -dbtype nucl -out nifh-db.fa

Note: The tilde ~ is a shortcut for your home directory /home/<username>.

Exercise 30 Look the files generared. Use ls command.

ls

BLAST

Now, we can perform a BLAST search of the gene sequence against this database. BLAST offers different algorithms depending on the type of query and the database:

  • blastn: Used for nucleotide sequence queries against a nucleotide database.
  • blastp: Used for protein sequence queries against a protein database.
  • blastx: Used for translating a nucleotide query to protein and comparing it against a protein database.
  • tblastn: Used for protein queries against a translated nucleotide database.
  • tblastx: Used for translating both the nucleotide query and the database sequences, comparing them in protein space.

To see the list of parameters of this command use:

blastn -help

For this practical, we need of theses options:

blastn -query <query_file> -db <db_file> -out <name_of_output> -outfmt <format>

Explanation of options:

  • query <query_file>: Specifies the input file with the sequence to search.
  • db <db_file>: Specifies the BLAST database to search against.
  • out <name_of_output>: Defines the output file for the results.
  • outfmt <format>: Sets the output format to tabular for easier parsing. Format 6 provides tabular output, which includes information like query sequence, subject sequence, percent identity, alignment length, and e-value. It is commonly used for further parsing and analysis of the results. We can have a good understanding of -outfmt 6 format from this webpage.

Launch BLAST for each query file.

blastn -query /data/blast_tp/metags/corn.fa -db nifh-db.fa -out corn.out -outfmt 6

blastn -query /data/blast_tp/metags/prairie.fa -db nifh-db.fa -out prairie.out -outfmt 6

blastn -query /data/blast_tp/metags/soybean.fa -db nifh-db.fa -out soybean.out -outfmt 6

Playing around BLAST output

Exercise 31 How many results do you get for each metagenome? Use the wc -l.

wc

The wc (word count) command in Unix/Linux is a utility used to count lines, words, and bytes in a file or input stream.

Here are its main options:

  • -l: Counts the number of lines.
  • -w: Counts the number of words.
  • -c: Counts the number of bytes.
  • -m: Counts the number of characters.
  • -L: Reports the length of the longest line.
wc [options] <file>
wc -l *.out

The wildcard * matches all files ending with .out.

To understand the meaning of each column in the BLAST output (format 6), refer to this webpage.

more

You can display the file with the command more.

The more command is used to display the contents of a file one screen at a time. This is useful for viewing large files like FASTQ files, where the data is too big to fit on a single screen.

more <path_file_to_display>

Exercise 32 What is the range of percentage identity for the results of BLAST of soybean metagenome?? Similarly, what is the range of e-values? Use morecommand.

more soybean.out

sort

The sort command is used to sort the lines of a file or input stream based on specified criteria. It can sort in ascending or descending order and handle various types of data.

sort [options] [file]

Common Options:

  • -n: Sort numerically.
  • -r: Sort in reverse order.
  • -k <key>: Specifies the key (column) to sort by.
  • -u: Remove duplicate lines.

cut

The cut command is used to extract specific sections or columns from each line of a file or input stream. It is particularly useful for parsing delimited files, such as CSV or tab-separated files.

cut [options] [file]

Common Options:

  • -d <delimiter>: Specifies the delimiter character (default is tab).
  • -f <fields>: Specifies the fields to extract (e.g., -f 1,3 to extract the first and third fields).

uniq

The uniq command is used to remove (or keep) consecutive duplicate lines from a sorted file or input stream. It is often used in conjunction with sort to ensure that the input is sorted.

uniq [options] [file]

Common Options:

  • -c: Prefix lines with the number of occurrences.
  • -d: Only print duplicate lines.
  • -u: Only print unique lines.

Exercise 33 How many uniq nifh genes were hit by sequences from soybean metagenome?

There are several ways to do this:

cut -f 2 soybean.out | sort | uniq | wc -l

cut -f 2 soybean.out | sort -u | wc -l

Exercise 34 How many uniq sequences from prairie metagenome can hit nifh genes?

There are several ways to do this:

cut -f 1 prairie.out | sort -u |wc -l

cut -f 1 prairie.out | sort | uniq | wc -l

Exercise 35 How to get the sseqid, sstart, and send column (column 2, 9, 10) from the soybean blast output file?

cut -f 2,9,10 soybean.out

grep

The grep command is used to search for patterns within files or input streams. It stands for “global regular expression print” and is highly versatile for text processing.

grep [options] [pattern] [file]

Common Options:

  • -i: Ignore case.
  • -v: Invert match (show lines that do not match the pattern).
  • -c: Count the number of matching lines.
  • -n: Show line numbers.
  • -r: Recursively search directories.

Exercise 36 How many times was gene gi|909637271|emb|LN713523.1| hit by sequences from corn metagenome sample:

grep "gi|909637271|emb|LN713523.1|" corn.out | wc -l