Browsing genomes with EnsemblGenomes

As reference genome, we will use the strain K12 MG1655 of Escherichia coli.

1. In a web browser, open a connexion to Ensembl Bacteria (https://bacteria.ensembl.org).

2. In the Search for a genome box, type MG1655. The interface will display a list of the genomes matching this string. Select the strain with taxonomic ID (TAXID) 511145.

EnsemblGenomes Bacteria organism choice.

3. After having selected your bacteria of interest, type aspartokinase in the search box and start the search.

EnsemblGenomes gene search.

4. Read the annotation of the three first matches. Note the names of the genes coding for these three enzymes.

EnsemblGenomes search result for aspartokinase.

Downloading genomic features and protein sequences

We will use blast on the command line to search for homologs of E.coli aspartokinases in other bacterial genomes.

Each student will chose a different bacteria to perform the search.

Downloading protein sequences and gene annotations from Ensembl genome

1. In a separate window, open a connection to EnsemblGenomes Bacteria (https://bacteria.ensembl.org/).

2. Click on the **Downloads* link (you can also access it directly here: https://bacteria.ensembl.org/info/website/ftp/).

3. In the filter box, type MG1655. This select the subset of genomes whose name contains this string.

EnsemblGenomes download filter.

4. Click on the FASTA link in the column Protein sequences/

5. and select Copy link address.

EnsemblGenomes Link to download all protein sequences for a selected organism.

5. Come back to the Unix terminal of the IFB NNCR.

6. Create a directory for the current practical.

mkdir -p ~/bacteria_blast
cd ~/bacteria_blast

7. type wget --no-clubber and paste the content of the clipboard. This should give the following command (note: at the tiem you run this tutorial, the release may be different, this one is from October 2018).

wget --no-clobber ftp://ftp.ensemblgenomes.org/pub/release-41/bacteria//fasta/bacteria_0_collection/escherichia_coli_str_k_12_substr_mg1655/pep/Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa.gz

If you want to know the meaning of the --no-clobber option you can checek wget manual with the command.

man wget

We will now check if the file has been correctly downloaded, and uncompress it.

## Check that the file is there
ls -l

## Uncompress the file
gunzip Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa.gz

# Measure the size of the uncompressed file
du -k Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa

# Check the 20 fist lines of the file
head -n 20 Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa

# Check the 20 last lines of the file
tail -n 20 Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa

Downloading genome annotations

After having downloaded the protein sequences, we will download gene annotations. In EnsemblGenomes, annotations are available in different formats, including GTF, a text file with columns separated by tabulations.

1. We In the browser, come back to the EnsemblGenomes Bacteria Download site, select MG1655 and follow the same instructions as above to download and uncompress the GTF file.

2. Uncompress the gzip archive.

## Download genome annotation file (GTF format)
wget --no-clobber ftp://ftp.ensemblgenomes.org/pub/release-41/bacteria//gtf/bacteria_0_collection/escherichia_coli_str_k_12_substr_mg1655/Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.41.gtf.gz

## Uncompress genome annotation file
gunzip Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.41.gtf.gz

## Check that the uncompressed file is there
ls Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.41.gtf

## Check the head of the annotaiton file
head -n 20 Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.41.gtf

## Check the tail
tail -n 10 Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.41.gtf

Exercise: genome annotations

As an exercise, you will now elaborate Unix commands to explore the content of E.coli genome annotations.

For this, we recommmend to use and combine the following commands. In order to fine-tune the parameters, read their manual on the Unix terminal by typing man [command_name] (for example man wc).

You can learn more Unix commands by typing “Linux cheat sheet” in any Web search engine.

We ask you to write all the answers in a text file, with comments, in order to enable other people to reproduce all your commands.

Activating the conda environment for blast on the NNCR cluster

In addition to facilitating the installation of software tools on different operating system, conda also permits to create modules that regroup a set of tools usually used together (for example NGS analysis, protein structure, proteomics, …).

The IFB NNCR cluster is already equipped with a series of conda environments with the classical bioinformatics tools. Before being able to use a tool, you must activate the corresponding environment.

## List the conda environments available on the system
$CONDA_HOME/bin/conda env list

## FInd the blast environment
$CONDA_HOME/bin/conda env list | grep -i blast

## Activate the blast environment. 
source $CONDA_HOME/bin/activate blast-2.7.1

## Check the location of the blastp command
which blastp

## Deactivate the current environment. 
source $CONDA_HOME/bin/deactivate
## Note the change in the shell prompt

## Re-activate the blast environment
source $CONDA_HOME/bin/activate blast-2.7.1
## Note the change in the shell prompt

Note: if you are working on a machine that does not belong to the federation of NNCR clusters, you can create the blast environment in your own conda folder with the following command (assuming that the conda package has been installed and that its executable files are in your path).

## ATTENTION: DO NOT RUN THIS COMMAND IF YOU ARE ON THE NNCR CLUSTER,
## because this would occupy a big space to reinstall BLAST in your 
## own account, although it is already available for all users. 

## Install blast
conda create -n blast-2.7.1 blast==2.7.1

## Install emboss
conda create -n emboss-6.6.0 emboss==6.6.0

Formatting the blast database for the reference genome

## Check the path of the program makeblastdb (it should be in your conda folder)
which makeblastdb

## Get the usage of makeblastdb (formal specification of the command syntax)
makeblastdb -h

## Get the help for makeblastdb (explanation of the options)
makeblastdb -help

## Build a BLAST database with all the protein sequences of E.coli
makeblastdb -in Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa  -dbtype prot

## Check the new files that were created with this commands. 
## For this we list the files in reverse (-r) temporal (-t) order. 
ls -1tr

The program makeblastdb created three files besides the input fasta file.

Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa.pin
Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa.phr
Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa.psq

These files contain an index of all the oligopeptites found in all the sequences of the fasta file. This “k-mer” index will be used by BLAST to rapidly find all the sequences in the database that match a query sequence, and use these short correspondences to start an alignment.

Read blastp manual

Since we want to search a protein database with protein query sequences, we will use the blastp command.

## Command  syntax
blastp -h

## Help
blastp -help

OK, the first contact is a bit frightening, because blastp has many options. This is because this program allows you to run queries in different modes, with different parameters, and to export the results in different formats.

In this tutorial we will show you the most usual ways to use the tool, and you will then be able to refine your search by combing back to the manual and understanding the use of additional options.

Downloading a query protein

We will first search similarities from a given protein of interest. We will start from the “askartokinase 3” protein from Escherichia coli, which is the product of the gene lysC.

Before starting the analysis, read the Uniprot record (https://www.uniprot.org/uniprot/P08660) in order to know the function and regulation of this protein.

We will use BLAST to identify all the proteins of Escherichia coli that have significant similarities with the aspartokinase 3.

For this we first need to download the sequence of this protein on our computer. The sequence can be downloaded directly from Uniprot.

wget --no-clobber https://www.uniprot.org/uniprot/P08660.fasta

Similiraty search for one sequence against a database

The following command search similarities for one or several query sequences provided in a fasta file (P08660.fasta) in a protein sequence database (Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa). By default, blastp prints the result on the screen, but we prefer to redirect it to a file in order to keep a trace of the result and inspect it later.

## Search similarities for aspartokinase 3 in E.coli proteome
blastp \
  -db Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa  \
  -query P08660.fasta \
  > lysC_vs_Ecoli.txt
  
## Check the result
less lysC_vs_Ecoli.txt

Getting a tabular synthesis of blast results

By default, BLAST displays the detailed results of a similaty search, starting with a summary of the matches, followed by all the alignemnts. It can be convenient to generate a synthetic result indicating only the relevant statistics for each alignment. This can be done by tuning BLAST options.

Read blastp manual and find a way to tune the formatting options in order to get a tabular output with comments.

## Search similarities for aspartokinase 3 in E.coli proteome.
## Return the result as a synthetic table with the matching statistics for each hit. 
blastp -outfmt 7 \
  -db Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa  \
  -query P08660.fasta \
  > lysC_vs_Ecoli_synthesis.tsv

## Check the result
more lysC_vs_Ecoli_synthesis.tsv

Shuffling sequences with EMBOSS shuffleseq

## Activate the conda environment with the emboss software suite
source  $CONDA_HOME/bin/activate emboss-6.6.0

## Check that shuffleseq is available
which shuffleseq

## Get the suffle manual 
shuffleseq -h

## Shuffle the apartokinase 3 sequence
shuffleseq -sequence P08660.fasta -outseq shuffled_lysC.fa

## Check the result
cat shuffled_lysC.fa

## Compare the shuffled seq with the original one
cat P08660.fasta

Here is an example of result (I just simplified the header of the shuffle output fasta sequence for the sake of readability).

Of course, your result should differ since each shuffling produces a randomized sequence.

>shuffled_lysC
VEAINNAVKRELSALGDVADLFVLEFDNDRLGWADAIIRVLLNTEMNGELKPAEVASAIV
GLIDEPELDTLVYSPLAGIFFLMCLIDFRGLGRREPTVPHSIMAGSGNALCRYFVTGLLL
LDTNAFAIRETKNLEAFVAEFQAADNFAQARNPANLPKLILCLSVNGVLQASPFFNLNEG
ADVVSESSAASLKADEDQSVITLEEDIVVTEGASSFDLALGLLMHLNVDRKVGAEWLEDE
VSPVMVEVLTLVLTQEALTVTSIALLRGDTRKRVGLVLNSFTSHDLIAPSLILIGTSLVT
RHFLVGAKVRIQERKASTSFRRLLAGRTETYRGLEIYAPCISMTTAVRRMQAGGIAHREI
QETIACRKSTELLSTIASPHAEVGFEAANGTLLGVLTADHVSMALQAVNGGESLPDEHRR
TVGDILTERVARLLLPFKEGPARTSSTLA

Exercise

1. Shuffle the aspartokinase 3 sequence with the tool shuffleseq (found in the conda environment emboss).

2. Run a blastp search with this shuffled sequence against all proteins of Escherichia coli.

## Run blast with the shuffled lysC sequence against all E.coli proteins
blastp \
  -db Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa  \
  -query shuffled_lysC.fa \
  > shuffled_lysC_vs_Ecoli.txt
## Check the result
less shuffled_lysC_vs_Ecoli.txt

## Run the same query with a tabular output
blastp -outfmt 7 \
  -db Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa  \
  -query shuffled_lysC.fa \
  > shuffled_lysC_vs_Ecoli_synthesis.tsv

## Check the result
less shuffled_lysC_vs_Ecoli_synthesis.tsv

Downloading a query genome

1. Open a connection to the EnsemblGenomes Download page https://bacteria.ensembl.org/info/website/ftp/

2. Identify the genome of your organism of interest (for this tutorial, I will choose Bacillus subtilis subsp. subtilis str. 168, but each student should take a different organism).

3. Download the genome annotation (GTF file) and protein sequences (fasta file).

## Download genome annotations for organism of interest
wget --no-clobber ftp://ftp.ensemblgenomes.org/pub/release-41/bacteria//gtf/bacteria_0_collection/bacillus_subtilis_subsp_subtilis_str_168/Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.41.gtf.gz

## Uncompress genome annotations for the organism of interest
gunzip -f Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.41.gtf.gz

## Download protein sequeces for query organism
wget --no-clobber ftp://ftp.ensemblgenomes.org/pub/release-41/bacteria//fasta/bacteria_0_collection/bacillus_subtilis_subsp_subtilis_str_168/pep/Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.pep.all.fa.gz

## Uncompress protein sequences for the organism of interest
gunzip -f Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.pep.all.fa.gz

## Check the result files
ls -l Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1*


## Format the blas database for the organism of interest
makeblastdb -in Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.pep.all.fa -dbtype  prot

Identifying aspartokinase homologs in the proteome of interest

To identify homologs for a given protein, we can run a simple blastp search against the proteome of the organism of interest (in my case, Bacillus subtilis).

blastp \
  -db Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.pep.all.fa \
  -query P08660.fasta \
  > lysC_vs_Bacillus_subtilis.txt
  
blastp -outfmt 7\
  -db Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.pep.all.fa \
  -query P08660.fasta \
  > lysC_vs_Bacillus_subtilis_outfmt7.tsv
  

However this would return many false positives, since BLAST sets its default E-value to 10 (we thus expect ~10 false positives on average for such a search). To avoid this, blastp enables to control the expected number of false positives with a threshold on the E-value, using the option -evalue.

blastp  -evalue 0.00001 \
  -db Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.pep.all.fa \
  -query P08660.fasta \
  > lysC_vs_Bacillus_subtilis_eval0.00001.txt

## Check the alignments of aspartokinase homologs
less lysC_vs_Bacillus_subtilis_eval0.00001.txt

blastp -evalue 0.00001 -outfmt 7 \
  -db Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.pep.all.fa \
  -query P08660.fasta \
  > lysC_vs_Bacillus_subtilis_eval0.00001_outfmt7.tsv

## Check the matching statistics for the aspartokinase homologs
cat lysC_vs_Bacillus_subtilis_eval0.00001_outfmt7.tsv

All proteins of B.subtilis against all proteins of E. coli

## Note : we prepend the command "time" to measure the computation time, just for the sake of curiosity
time blastp -outfmt 7 \
  -db Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa  \
  -query Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.pep.all.fa \
  > Bacillus_subtilis__vs__Escherichia_coli_blastp_outfmt7.tsv

## Create a separate file without comments
grep -v '^#' Bacillus_subtilis__vs__Escherichia_coli_blastp_outfmt7.tsv > Bacillus_subtilis__vs__Escherichia_coli_blastp_table.tsv  

## Count the number of hits
wc -l Bacillus_subtilis__vs__Escherichia_coli_blastp_table.tsv  

## Note : we prepend the command "time" to measure the computation time, just for the sake of curiosity
time blastp -outfmt 6 -evalue 0.00001 \
  -db Escherichia_coli_str_k_12_substr_mg1655.ASM584v2.pep.all.fa  \
  -query Bacillus_subtilis_subsp_subtilis_str_168.ASM904v1.pep.all.fa \
  > Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv

## Count the number of hits
wc -l Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv

## Identify  the column containing the ID of the hits in the proteome of interest
head Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv

## Cut the first column (containing the ID of the hits in the query proteome)
cut -f 1 Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv | more

## Count the number of hits
cut -f 1 Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv | wc -l

cut -f 1 Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv | sort -u | wc -l

cut -f 1 Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv | sort | uniq | wc -l

cut -f 1 Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv | sort | uniq -c | less

cut -f 1 Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv | sort | uniq -c |sort -nr -k 1 | less

cut -f 1 Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv | sort | uniq -c |sort -nr -k 1 | head -n 10

## A first attempt to sort by E-value, not very convincint !
sort -k 11,11 -n Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv | less 

**Note: after hacing seen the result, type "q"** to quit `less`.

...
CAB16047        AAC73707        56.673  517     ...   517     0.0     585
CAB16048        AAC75939        65.828  477     ...   479     0.0     676
CAB16138        AAC76764        52.564  624     ...   629     0.0     635
CAB11917        AAC76323        56.034  116     ...   116     1.00e-47        146
CAB12080        AAC76505        28.638  213     ...   229     1.00e-21        87.8
CAB12153        AAC73552        26.360  239     ...   564     1.00e-12        64.3
CAB12183        AAC73788        34.934  229     ...   223     1.00e-36        125
CAB12414        AAC74372        26.061  165     ...   195     1.00e-07        50.4
CAB12414        AAC76505        24.895  237     ...   223     1.00e-13        68.9
CAB12437        AAC74386        29.075  227     ...   220     1.00e-13        65.1
CAB12735        AAC76480        22.707  229     ...   242     1.00e-13        65.5
...

## Sort BLAST hits on a column with significance, computed with awk
awk '{sig= -log($11) / log(10); print $0"\t" sig}' \
  Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv \
  | awk '{ sub("inf", "200"); print }' \
  | sort -k 13,13 -n -r \
  | less 

- Primary key: decreasing significance (as computed above) 
- Secondary key: decreasing bit score (the $12^{th}$ column of BLAST output)

## Define an environment variable with the name of the sorted hits (this will be useful for further processing)
BLAST_HITS_SORTED=Bacillus_subtilis__vs__Escherichia_coli_blastp_sorted.tsv

## Sort BLAST hits by decreasing significance and decreasing bit score
awk '{sig= -log($11) / log(10); print $0"\t" sig}' \
  Bacillus_subtilis__vs__Escherichia_coli_blastp_eval0.00001.tsv \
  | awk '{ sub("inf", "200"); print }' \
  | sort -k 13,13 -k 12,12 -n -r \
  > ${BLAST_HITS_SORTED}
  
## Check the result
less ${BLAST_HITS_SORTED}

- a rank of the hits according to the query protein (`qrank`),
- another rank according to the "subject" protein (`srank`), i.e. the protein found in the blast database of your reference organism. 

awk '{ qrank[$1]++ ; print $0"\t"qrank[$1]}'  ${BLAST_HITS_SORTED} | less

BLAST_HITS_RANKED=Bacillus_subtilis__vs__Escherichia_coli_blastp_ranks.tsv
awk '{ qrank[$1]++ ; srank[$2]++ ; print $0"\t"qrank[$1]"\t"srank[$2]}' \
  ${BLAST_HITS_SORTED} > ${BLAST_HITS_RANKED}
  
## Check the result
less ${BLAST_HITS_RANKED}

- the best hit for each protein of the query organism (query BH);
- the best hit for each protein of the reference organism (subject BH);
- the bidirectional best hit (BBH)

## Best hits per query sequence
BLAST_QUERY_BH=Bacillus_subtilis__vs__Escherichia_coli_blastp_query_BH.tsv
awk '$14 == 1' ${BLAST_HITS_RANKED} > ${BLAST_QUERY_BH}

## Best hits per subject sequence
BLAST_SUBJECT_BH=Bacillus_subtilis__vs__Escherichia_coli_blastp_subject_BH.tsv
awk '$15 == 1' ${BLAST_HITS_RANKED} > ${BLAST_SUBJECT_BH}

## Bidirectional best hits
BLAST_BBH=Bacillus_subtilis__vs__Escherichia_coli_blastp_BBH.tsv
awk '$14 == 1 && $15 == 1' ${BLAST_HITS_RANKED} > ${BLAST_BBH}

## Count number of best hits
wc -l  ${BLAST_QUERY_BH} ${BLAST_SUBJECT_BH} ${BLAST_BBH}

All proteins from E. coli versus E. coli

Exercise: build a command that compares the sequences of all the E.coli proteins against each other, and store the results in a tab-delimited file with a single line per hit (we do not need the alignments).

Informations

• Ce devoir-maison peut être réalisé en combinant quelques commandes en ligne (cut, grep, awk, sort, uniq, wc, diff) et l’outil bioinformatique blastp.

  Nous vous recommandons de lire les manuels de ces commandes, ou de consulter quelques tutoriels. En particulier, il sera utile de connaître la syntaxe de awk pour manipuler des chaînes de caractères (substitution, impression de colonnes sélectionnées) et pour numéroter les lignes. On trouve par exemple des informations ici: http://www.catonmat.net/blog/awk-one-liners-explained-part-two/.

• Le rendu de ce TP inclut un script bash rédigé de façon à être compréhensible par un lecteur non averti (commentez le code, nommez les variables de façon indicative, …). Ce script devra également être paramétrisable pour analyser d’autres organismes que le vôtre. A cette fin, définissez dès le début des variables pour les organismes d’intérêt (QUERY_ORG) et de référence (REF_ORG), ains que pour les noms de fichiers d’input et output. En principe l’enseignant devra pouvoir refaire tourner toutes vos analyses à partir de votre script.

• La seconde pièce du rendu sera un rapport d’interprétation synthétique.

• Ci-dessous (et dans votre rapport) nous désignerons les hits avec une E-valeur inférieure à 1e-5 comme “homologues prédits” ou plus simplement “homologues”.

Modèle de rapport

Je fournis ci-dessous 2 modèles de rapport pour le devoir-maison du TP proteome BLAST.

• Format docx (peut s’ouvrir avec LibreOffice ou avec un éditeur de texte commercial).

• Pour les geeks: un fichier R markdown (extension Rmd), que j’ai utilisé pour générer le modèle de rapport docs. Ce fichier Rmd peut être édité et compilé avec RStudio de façon relativement intuitive, pour produire un rapport HTML, pdf ou docx. Comme vous n’avez pas encore eu l’occasion d’apprendre R markdown, il n’y a aucune obligation à utiliser ce format, mais je vous le fournis déjà pour ceux qui seraient tentés de se lancer. Si c’est votre cas, soumettez à la fois le fichier Rmd et le résultat de votre compilation (html, word ou docx).

• Modèle de rapport: docx Rmd

Questions

1. Faites tourner blast pour toutes les protéines de votre organisme d’intérêt contre toutes celles de E.coli, en sélectionnant un fichier de sortie tabulaire sans commentaires, et un seuil de \(10^{-5}\) sur la E-valeur. Triez les hits par signiticativité (clé primaire) et par bit score (clé secondaire), puis calculez les rangs des hits en fonction des séquences “query” et “subject” comme nous l’avons fait dans le tutoriel.

2. Comptez

   • le nombre total de hits
   • le nombre total de paires de protéines avec un hit
   • le nombre de protéines de votre ortanisme d’intérêt ayant au moins un hit chez E.coli
   • le nombre de protéines d’E.coli pour lesquelles au moins un hit a été trouvé à partir de votre organisme d’intérêt
   • le nombre de best hits pour l’organisme d’intérêt
   • le nombre de best hits pour l’organisme de référence (E.coli)
   • le nombre de best hits bidirectionnels

3. Générez un fichier qui indique le nombre d’homologues prédits par protéine de l’organisme d’intérêt. Ce fichier contiendra une ligne par protéine, et deux colonnes, indiquant respectivement le nombre de hits et l’ID de la protéine (uniq -c). Triez les résultats par nombre décroissant de hits (les protéines ayant le plus grand nombre de hits s’affichent en premier lieu). Afficher les 10 premières lignes du fichier.

Question d’interprétation : Utilisez Uniprot pour savoir quelle est la fonction biologique des protéines présentant un très grand nombre de hits. Cette fonction permet-elle de comprendre ce nombre excessif de hits ?

4. Calculez le nombre de protéines de l’organisme d’intérêt sans homologue prédit.

Question d’interprétation : Comment interprétez-vous ces protéines sans hit ?

5. Calculer la distribution du nombre d’homologues par protéine (combien de protéines ont 0, 1, 2, 3, … homologues).

Question d’interprétation : quel est le nombre moyen de hits ? Combien de protéines ont un seul hit ? Comment interpréteriez-vous ces deux nombres en termes évolutifs ?

6. Identifiez dans les deux comparaisons (votre organisme vs E.coli et E.coli vs votre organisme) toutes les homologies concernant des aspartokinase. Sur base de la règle du bidirectional best hit, tentez d’établir les relations d’orthologie et de paralogie entre chaque paire de protéines, et représentez les résultats dans un tableau avec les gènes de votre organisme sur les lignes, les gènes de E.coli sur les colonnes, et en indiquant “P”, “O”, ou “I” pour paralogie, orthologie ou identité, respectivement.

Question d’interprétation : comment interprétez-vous ces résultats en terme de voies métaboliques ? N’hésitez pas à constuler des ressources complémentaires (KEGG, MetaCyc, …).