Molecular Modelling and GNU/Linux

Python : bonnes pratiques pour mettre en place un programme

S. Téletchéa — 2024-04-19T09:39:45Z

Sommaire

La langage python est récemment devenu très populaire de par sa facilité apparente d'utilisation. Il reste néanmoins un langage informatique qui comporte des notions complexes, et une fois les premiers programmes écrits (scripts), il faut assez rapidement réfléchir à structurer le programme pour qu'il devienne exploitable sur le moyen terme, au risque sinon de devenir du code mort, jetable, et donc perdu.

Pour éviter cette dérive, voici quelques conseils avec des exemples rapides d'utilisation.

Version 1 du programme, un "script"

Prenons comme exemple un petit programme en python suivant, qui permet d'afficher un simple "coucou" dans un terminal, que l'on enregistre dans le fichier programme_v1.py :

print("Coucou")

Pour voir le résultat ce petit programme, on peut l'exécuter ainsi :

python programme_v1.py

... qui affiche donc dans le terminal

Coucou

Version 2 : programme, le rendre exécutable

Cela peut devenir vite fastidieux de mettre à chaque fois le mot "python" devant chacun des programmes que l'on utilise, alors qu'au final c'est le programme dont on a besoin. Il existe une commande pour changer cela, mais qui nécessite de revoir un peu la structure de notre petit programme, il faut lui ajouter un shebang à la première ligne.

#!/usr/bin/env python3 print("Coucou")

Le nouvelle programme s'appelle maintenant programme_v1.1.py.

Remarquez que comme il s'agit d'une version mineure du programme, j'ai incrémenté son numéro simplement (passage de 1 à 1.1). Nous reviendrons sur cette notion plus loin.

Pour rendre ce programme exécutable et le lancer, il faut taper :

 chmod +x programme_v1.1.py ./programme_v1.1.py Coucou

Remarquez qu'il faut préciser l'emplacement du fichier pour que le système lance le programme, en spécifiant "./" au début du nombre du programme. Le point suivi d'une barre oblique "./" indique que l'on utilise le répertoire courant, et qu'il faut exécuter le programme "programme_v1.1.py".

Le résultat ne change pas.

Version 2 du programme, du script au programme

Le "script" par définition est un petit programme qui a vocation à être court, pour répondre à un besoin précis, on peut souvent le réécrire sans trop reprendre un ancien script. C'est l'une des raisons de l'utilisation du shell bash, c'est qu'il permet d'enchaîner plusieurs traitements rapidement, par exemple faire le déplacement et le renommage de milliers de fichiers.

Si l'on utilise python c'est souvent parce que le traitement que l'on envisage va faire beaucoup plus d'étapes qu'envisagé initialement. Il est donc utile de mieux organiser son programme dès le départ, en suivant quelques bonnes pratiques.

La version 2 de notre programme ne va pas faire évoluer grandement encore notre programme, mais elle va permettre de le structurer.

Le fichier programme_v2.py devient ainsi à partir du fichier programme_v1.1.py :

#!/usr/bin/env python3 def main() : print("Coucou") if __name__ == "__main__" : main()

Quand on rend exécutable le programme et qu'on le lance, le résultat est toujours le même, pour le moment.

Version 3 : découper les fonctionnalités du programme

Un programme est un ensemble d'étapes que l'on veut faire faire à un ordinateur. La tâche la plus compliquée est de découper les actions pour qu'elles soient réutilisables dès le début, c'est à la fois utile pour clarifier les étapes à exécuter, mais aussi pour mieux isoler les endroits du programme où il faudra réfléchir sur le meilleur moyen d'arriver au résultat, qui le bien souvent doit faire appel à un algorithme simple ou complexe.

Pour réaliser le programme, on va donc le découper en fonction, chaque fonction va contenir des commentaires au format docstrings.

# !/usr/bin/env python3 def affiche_coucou() : """ Fonction qui permet d'afficher le texte Coucou """ print("Coucou") def main() : affiche_coucou() if __name__ == "__main__" : main()

Après avoir donné les permissions d'exécution au programme, le texte "Coucou" s'affiche à nouveau dans le terminal quand lance le programme.

Documenter son programme

L'intérêt d'avoir mieux organisé son programme est qu'il est maintenant plus modulaire, ce qui permet de le rendre plus visible, donc de produire une documentation !

Il existe plusieurs moyens de documenter son programme, mais une manière simple est de faire appel à pydoc3.

Pour l'installer sur son système il faut utiliser la commande pip :

pip install pdoc3

Cela va créer dans votre espace de travail ($HOME) une commande dans le répertoire $HOME/.local/bin.

Il suffit maintenant de générer un fichier HTML du code en version 3 :

$HOME/.local/bin/pdoc3 programme_v3.py —html —force

Cette commande va générer un fichier au format HTML qu'il est possible de regarder en ligne directement avec un navigateur.

Version 4 : continuer à développer son programme

Une fois la première ébauche du programme réalisée, il faut commencer à réfléchir à toutes les étapes que l'on voudrait réaliser, les documenter avec des docstrings, puis ensuite les écrire une par une. Une fois cette implémentation préliminaire réalisée, il sera possible d'avoir des programmes python très complexes, bien documentés et réutilisables.

# !/usr/bin/env python3 def affiche_texte(letexte) : """ Fonction qui permet d'afficher le texte fourni en argument """ print(letexte) def demande_texte() : """ Fonction qui demande un mot à afficher à l'utilisateur, pour l'instant elle ne fait rien """ pass def main() : demande_texte() affiche_texte("Coucou") affiche_texte("Bonjour") if __name__ == "__main__" : main()

Dans cette version v4 du programme la fonction affiche_coucou() a été renommée et modifiée pour devenir affiche_texte(letexte), elle devient ainsi plus générique car l'argument passé à la fonction (enregistré dans la variable letexte) est le texte à afficher, et non une valeur statique.

Une seconde fonction a été ajoutée (demande_texte()), elle ne fait pour l'instant rien, mais aura comme objectif de demander à l'utilisateur le mot à afficher. Comme il y a plusieurs moyens pour interroger un utilisateur, en ligne de commande, ou via une interface graphique par exemple, ce point ne sera pas détaillé ici.

En conclusion

Le petit exemple présenté dans un article fait partie des bonnes pratiques qu'il est important de suivre à chaque fois que l'on écrit un programme en python. Avec un peu de rigueur cela permet d'obtenir rapidement un code propre, clair et fonctionnel. Cela permet aussi d'obtenir très rapidement un résultat qui correspond à ce dont on a besoin, c'est donc en plus d'une bonne pratique, une méthode efficace de développement.

Comme indiqué précédemment, un des éléments importants pour éviter d'avoir x version du même programme est de le versionner de manière incrémentale. Cependant cela devient vite compliqué de fonctionner ainsi, c'est pour cela qu'il vaut mieux tout de suite utiliser un gestionnaire de version comme git avec un hébergement sur github, gitlab ou encore un dépôt local. Cela vous permettra de mieux suivre vos changements, même en travaillant seul sur un projet, en faisant des changements importants entre chaque commit.

Le code détaillé rapidement ici est ainsi présent sur mon dépôt github public : https://github.com/TeletcheaLab/pythondemo/.

XmGrace day-to-day usage

S. Téletchéa — 2023-08-17T20:04:32Z

Before finding the wonderful XmGrace software, i'd loose time to do calculations, regressions, etc.

With this software, i've been able to do in a more robust manner and more rapidly everything i could do before and much more.

Here are some examples of what i'm using.

Some basics commands i use with XmGrace are detailed below.

Renormalizing a plot (mostly angles between -180 and +180, you want it between 0 and 360 )

Load your file into xmgrace
Go to Data->Transformations->Evaluate expressions
Select G0.S0
Fill in the box with :

y = y<=0 ? y + 360 : y

It is a basic if then else : y equals the value resulting of the test y<=0. If the test is true then y = y+ 360 otherwise take the y value

Subpopulation (for instance for an angle, where are the points above a certain value)

Load your file into xmgrace
Go to Data->Transformations->Sample points
Select the desired set (presumably S0)
Change sample type for 'Expression'
Fill in the form 'Logical expression' with your needs, for instance y=>50 will select the population where all points are equal or above 50.

Correlation between two states (for instance a sugar conformation and a dihedral angle)

Load file1 into xmgrace (it will appear as set 0)
Load file2 into xmgrace (it will appear as set 1)
Go to Data->Transformations->Evaluate expressions
Select g0.s1
Fill in the box with :
- copy s1 to s2
- s2.x=s0.y
Accept

Go back to the graph, hide s0 and s1, click on As (autoscale), you have your correlation :-) you can now fill in legends, etc.

Getting informations about points (mean, std, max, min, ...)

Load your file into xmgrace
Go to Edit->Data sets
Pick your set, the values will be displayed in the box below (rapid, isn't it ?)

Sending data through a pipe

echo -e "1 2\n2 4\n3 6" |xmgrace -pipe

Sending data through a fifo

Create the fifo

mkfifo fifo.pipe
xmgrace -npipe fifo.pipe

Send the commands through the fifo

echo "read \"data.out\"" > fifo.pipe
echo autoscale > fifo.pipe
echo redraw > fifo.pipe

Other sites with tips

http://ringo.ams.sunysb.edu/index.php/Xmgrace

How popular is a programming language (in bioinformatics) ?

S. Téletchéa — 2023-06-29T14:29:38Z

Sommaire

If in doubt, go for Python :-)

There are many ways to compare the "popularity" of programming languages used nowadays. It can be monitored from the demands of companies nearby you, from the relationships you may have access to (aka "the developpers"), or even by following the guidance of university teachers.

Since this is a complex problem, there is no easy answer to this question. For a rapid comparison, you can for instance use google trends to compare PHP and Python queries over time.

This google search bounds your analysis to google queries. Albeit being popular, they only capture part of the problem... How to get better data then ?

The TIOBE index

To compare programming data, why not asking a company dedicated to providing development solutions for hackers ?

There comes the TIOBE index : https://www.tiobe.com/tiobe-index.

From this index, Java, C, and python are the TOP 3 programming languages in september 2018. Python usage is rapidly and steadily increasing, as was also remarked on Stack Overflow last year (https://stackoverflow.blog/2017/09/06/incredible-growth-python/).

Two different sources, same conclusion, Python usage is growing. What if there was some bias from those sources ?

PYPL - PopularitY of Programming Language

If you need another way of getting the trends, PYPL analyses differently the data found. I'm not to comment on this approach, but you can find the details in their FAQ and explanations why their results are slightly different from others.

In their rankings, Python comes first, followed by Java (decreasing) and JavaScript (increasing). In short, the ranking is different, but Python is again important.

In conclusion

This rapid post is there to keep a track of reference sources where you can find the information on programming languages trends.

As of September 2018, high-level languages are still getting more and more popular, as Java+Python+PHP represent more than 50 % of the total languages referenced (PYPL) or 25% of engineer languages (TIOBE). Surprisingly, they are just followed by very-old-but-thought-to-be-obsolete C/C++ (6 % for PYPL, 22% for TIOBE).

If you understand the programming paradigms pertaining to this small list of languages + some practical expertise on some of them, it seems you are ready for the future, i.e. the next 5 years in the computational world (https://prezi.com/_mu9uml7kaap/computational-world/).

Parallélisation rapide de code (bio)python

S. Téletchéa — 2023-06-29T09:17:05Z

Sommaire

Comment exécuter facilement un programme plus rapidement sur un ordinateur moderne ?

Voici les commandes rapides pour transformer un programme monotâche en python en programme capable d'utiliser tous les processeurs de votre machine.

Programme exemple monotâche en biopython

Le code simplifié ci-dessous permet de lire l'ensemble des transcrits du transcriptome de Physarum polycephalum autrement connu sous le nom de BLOB, vous pouvez trouver une présentation sur YouTube de ce petit organisme par Audrey DUSSOTOUR : https://www.youtube.com/watch?v=wjzEMIox_a8.

#!/usr/bin/env python3 from Bio import SeqIO def processTranscrit(identifiant): try: r = SeqIO.read(f"{identifiant}.fasta", "fasta") print(r.id) except Exception as e: print(f" **** Erreur pour traiter {identifiant} : {e}") def main(): proteome = "proteome_physarum.fasta" for r in SeqIO.parse(proteome, "fasta"): processTranscrit(r.id) if __name__ == "__main__": main()

Ce programme va lire le fichier avec biopython (voir une excellente introduction par Olivier Dameron et Grégory Farrant : https://inria.hal.science/hal-01095475 ), extraire les séquences d'intérêt une par une, puis les traiter avec la fonction processTranscrit qui pour l'exemple ne fait que l'affichage de l'identifiant de la séquence.

Programme parallélisé

Pour paralléliser notre code où pour le moment une boucle monotâche est présente dans la fonction main(), il suffit de faire appel à la bibliothèque multiprocessing de python en début de code :

# Multiprocessing du calcul # https://docs.python.org/fr/3/library/multiprocessing.html import multiprocessing as mp from multiprocessing import Process, Pool, TimeoutError

La boucle de la fonction main() faisant appel à la fonction processTranscrit(identifiant) devient ainsi :

 pool = mp.Pool(mp.cpu_count()/2) pool.map(processTranscrit, [r.id for r in SeqIO.parse(proteome, "fasta")]) pool.close()

Il faut tout d'abord créer un Pool qui va indiquer le nombre de tâches à exécuter, soit le nombre de cœur divisé par deux dans mon cas, pour n'utiliser que les processeurs physiques, ajouter toutes les tâches à réaliser avec la fonction map qui fait appel aux compréhension de listes en python, puis de lancer le calcul avec la méthode close().

Le code s'exécute ainsi en parallèle sur 4 tâches, threads en anglais, sur ma machine, le temps de calcul initial qui dépassait une journée est ainsi descendu à moins de 6 heures...

Pour aller plus loin

Une explication plus complète et en anglais est disponible sur le site https://www.analyticsvidhya.com/blog/2021/04/a-beginners-guide-to-multi-processing-in-python/

Making use of BLAST, example on the COVID-19 genome

S. Téletchéa — 2023-06-28T15:50:20Z

Sommaire

BLAST is a very powerful and widely used method for querying large sequence databases for sequence proximity with a sequence of interest. This difficult problem requires a fast algorithm designed for nucleic and protein sequences only since their alphabet is limited (5 bases for nucleotides, 21 amino acids for protein), but sequence length can be very long (more than 40000 amino). An example of BLAST installation and usage is provided below.

BLAST Installation

Basic Local Alignment Search Tool executables can be found online from the NCBI server : https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ but if you are using a Linux-based system, it is often as easy as

 sudo apt install ncbi-blast+

Once downloaded and installed, you'll get a list of software for performing your sequence queries.

blastn : searching in a nucleotide database from a nucleotide query
blastp : searching in a protein database from an amino acids query
blastx : searching in a protein database from a translated nucleic sequence
tblastn : searing in a nucleotide database from a translated amino acids query
When possible, I use protein queries on protein databases, since this is more efficient :

three nucleotides define a codon, and thus one amino acid
some mutation are synonymous, i.e. they do not change the amino acid

Data retrieval : Sars-Cov-2 genome, Sars-Cov-1 Hemaglutinin-Esterase

All SARS-Cov-2 sequences are regrouped on the COVID-19 NCBI dedicated portal, and more specifically, the initial SARS-Cov-2 genome deposited on January, 13th, 2020 is accessible under the reference NC_045512.

This sequence is coming from the discovery of SARS-Cov-2 published in Nature in February : Wu, F., Zhao, S., Yu, B. et al. A new coronavirus associated with human respiratory disease in China. Nature 579, 265–269 (2020). https://doi.org/10.1038/s41586-020-2008-3

As of writing the genome has been revised on July, 18th, 2020. The reference SARS-Cov-2 reference Genome is thus NC_045512.2. The core sequence has not been modified, only the virus name.
The complete nucleotide sequence of the virus is provided below in zip format for convenience :

https://web.archive.org/web/20230325065144/http://www.steletch.org/IMG/zip/sars-cov-2-nc_045512.2.fasta.zip

SARS-Cov-2 complete genome

Since we want to search for the presence (or absence) of and Hemalutinin-Esterase gene, as found in other SARS-Cov genomes (or in other viral genomes), we take the reference protein sequence of this protein from the Human coronavirus OC43 (HCoV-OC43) strain, entry A0A482K8Y4.
The sequence (in fasta format) of the protein is easily accessible under the "Sequence" tab in Uniprot, but we also provide the link here for convenience (again as a zipped file).

https://web.archive.org/web/20230325065144/http://www.steletch.org/IMG/zip/a0a482k8y4.fasta.zip

Hemagglutinin-esterase of Human coronavirus OC43

Preparing the genome for BLAST

To execute the search in a reasonable time, BLAST needs to have the database sequence mostly in fasta format (text) pre-processed into 3 binary files :

nsq : sequence data
pin : sequence indices
phr : complete identifiers

The command to execute is :

makeblastdb -in sars-cov-2-NC_045512.2.fasta -out sars2Genome -dbtype nucl

The result of this command is :

 Building a new DB, current time: 11/20/2020 13:41:28 New DB name: /home/stephane/Nextcloud/Enseignement/Bioinformatique/BLAST/sars2Genmoe New DB title: sars-cov-2-NC_045512.2.fasta Sequence type: Nucleotide Keep MBits: T Maximum file size: 1000000000B Adding sequences from FASTA; added 1 sequences in 0.000279903 seconds.

Since we took the raw nucleotide genome, there is only one (big !) sequence found and indexed.
Please note that we have very restricted use of the makeblastdb command, more options are available in the online documentation.
To search for protein results from a protein query into a nucleotide query we use tblastn (as indicated above) :

 blast -query A0A482K8Y4.fasta -db sars2Genome

Analyzing the Results

The BLAST output by default is quite verbose. The complete log is available on the following file :

https://web.archive.org/web/20230325065144/http://www.steletch.org/IMG/zip/hequery.blast.zip

BLAST result on SARS-Cov-2, query from Hemaglutinin-Esterase of HCoV-OH43

An extract of the results are presented below :

 > NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome Length=29903 Score = 21.6 bits (44), Expect = 1.3, Method: Compositional matrix adjust. Identities = 9/24 (38%), Positives = 14/24 (58%), Gaps = 6/24 (25%) Frame = -2 Query 27 VSHVNGDWFL------FGDSRSDC 44 +SH++ W+L G +RSDC Sbjct 24274 ISHLHSKWYL*CSTCTKGPTRSDC 24203

The first line indicates the identifier found for the identified positions. Since we did not process the genome, we only have the entire genome and not the annotated one (containing ORFs, transcripts, and proteins).

The line starting with "Score" contains BLAST information concerning the reliability of the search, If you are not familiar with these scores, you can have a look at the NCBI YouTube channel.

Important information to consider is the Expect value (1.3, which is not good), but also the Frame, -2, indicative that an offset is required from the initial genome translation site to find the open reading frame.
The lines starting with Query and Sbject (Subject) indicate where a match between the query sequence (the Hemaglutinin-Esterase from HCov-OC43) and the reference genome is found.

This shows that in the HE sequence, from position 27 to position 44, a match is found in the SARS-CoV-2 genome at positions between 24274 and 24203. Please remark that the numbering is reversed in the output, this is indicative of the putative presence os this sequence on the reverse-complement strand of the RNA.

The rest of the BLAST output is shown below.

 Score = 21.6 bits (44), Expect = 1.4, Method: Compositional matrix adjust. Identities = 7/18 (39%), Positives = 11/18 (61%), Gaps = 0/18 (0%) Frame = -2 Query 271 TKAICFNKRKDFTPVQVV 288 +K IC N +DF V ++ Sbjct 23452 SKLICMNSNRDFCAVNIL 23399 Score = 21.2 bits (43), Expect = 1.9, Method: Compositional matrix adjust. Identities = 12/46 (26%), Positives = 22/46 (48%), Gaps = 3/46 (7%) Frame = -3 Query 140 YRSLTFVNVPYVYNGSAQ---ATAFCKSGSLVLNNPVYIAPQANSG 182 Y +L + V +YN + CKS V++ +++ PQA + Sbjct 19686 YFTLLSIKVSLIYNIKSHIF*TLKSCKSVPEVISIKLFVNPQAKAS 19549 Score = 21.2 bits (43), Expect = 1.9, Method: Compositional matrix adjust. Identities = 10/27 (37%), Positives = 14/27 (52%), Gaps = 6/27 (22%) Frame = +1 Query 13 IIGSLGFHNPPTNVVSHVNGDWFLFGD 39 I+G L N +NG+W+ FGD Sbjct 14041 IVGVLTLDN------QDLNGNWYDFGD 14103 Score = 19.2 bits (38), Expect = 7.1, Method: Compositional matrix adjust. Identities = 8/27 (30%), Positives = 13/27 (48%), Gaps = 0/27 (0%) Frame = +2 Query 288 VDSRWNNARQSDNMTAVACQPPYCYFR 314 +D RWN R + Q +C+F+ Sbjct 26744 LDHRWNCYRNGLSCRLDVAQLLHCFFQ 26824

In short, it is possible to find very tiny portions of the HE protein at the genomic level of SARS-Cov-2 (although barely significantly). If you still want to go deeper into BLAST interpretation, here is a nice video on it !

What happens ?

If we glue up the results, this means that we can find some matches in the SARS-Cov-2 genome, respectively on 24 amino acids + 18 + 46 +27 + 27, for a total of 142, or about one-third of the total protein length. But this difficult to imagine. This would require the virus polymerase to do multiple back-and-forth movements between frames (-2, then -3, then +1 and +2).

This seems very unlikely...

On a side note, there have been examples of Quality Control RNA Replication or Transcription "slippery" of coronaviruses, but if you want to know more about the story, please see the wonderful video of Pr. Britt Glaunsinger from Berkeley about that and coronaviruses, in general.

Going Further

To sum up, we discovered how to install and use BLAST, then applied it to an existing question and answered it. There is no Hemaglutinin-Esterase in SARS-CoV-2, contrary to some other beta coronaviruses.
If you still want to know more about BLAST but also about SARS-CoV-2 ongoing research efforts, you can visit the following links :

23 new ORFs detected by ribosome profiling. Finkel, Y., Mizrahi, O., Nachshon, A. et al. The coding capacity of SARS-CoV-2. Nature (2020).
Mutations : D614G (and more ?)
Other mutations and large virus follow-up
Deep docking protocol for large screening of potential inhibitors (1.3 billion compounds !)
Become a JEDI and participate in the 1 billion COVID challenge
Join the ongoing fight against COVID-19 research efforts in Europe https://www.exscalate4cov.eu/ or worldwide https://covid.postera.ai/
Structure information listed by the NCBI : https://www.ncbi.nlm.nih.gov/Structure/SARS-CoV-2.html
All SARS-CoV-2 proteins with an atomic structure (and more) : https://www.rcsb.org/news?year=2020&article=5e74d55d2d410731e9944f52&feature=true
The dedicated Uniprot Portal for COVID-19 : https://covid-19.uniprot.org/

Molecular Modelling and GNU/Linux

Python : bonnes pratiques pour mettre en place un programme

Sommaire

Version 1 du programme, un "script"

Version 2 : programme, le rendre exécutable

Version 2 du programme, du script au programme

Version 3 : découper les fonctionnalités du programme

Documenter son programme

Version 4 : continuer à développer son programme

En conclusion

XmGrace day-to-day usage

Table of contents

Renormalizing a plot (mostly angles between -180 and +180, you want it between 0 and 360 )

Subpopulation (for instance for an angle, where are the points above a certain value)

Correlation between two states (for instance a sugar conformation and a dihedral angle)

Getting informations about points (mean, std, max, min, ...)

Sending data through a pipe

Sending data through a fifo

Other sites with tips

How popular is a programming language (in bioinformatics) ?

Sommaire

The TIOBE index

PYPL - PopularitY of Programming Language

In conclusion

Parallélisation rapide de code (bio)python

Sommaire

Comment exécuter facilement un programme plus rapidement sur un ordinateur moderne ?

Programme exemple monotâche en biopython

Programme parallélisé

Pour aller plus loin

Making use of BLAST, example on the COVID-19 genome

Sommaire

BLAST Installation

Data retrieval : Sars-Cov-2 genome, Sars-Cov-1 Hemaglutinin-Esterase

SARS-Cov-2 complete genome

Preparing the genome for BLAST

Analyzing the Results

What happens ?

Going Further