Script bash pour réaliser un calcul toute les 51 lignes

[Oceane666]

Bonjour à tous !

Je possède fichier nommé : methylation_extracted_downstre am.bed
Avec ce type de données :
contig_1 21914 21915 0 0 0
contig_1 21915 21916 0 0 0
contig_1 21916 21917 12 7 132
contig_1 21917 21918 0 0 0
contig_1 21918 21919 0 0 0
contig_1 21919 21920 0 0 0
contig_1 21920 21921 0 0 0
contig_1 21921 21922 0 0 0
contig_1 21922 21923 0 0 0
contig_1 21923 21924 0 0 0
contig_1 21924 21925 0 0 0
contig_1 21925 21926 0 0 0
contig_1 21926 21927 0 0 0
contig_1 21927 21928 0 0 0
contig_1 21928 21929 0 0 0
contig_1 21929 21930 0 0 0
contig_1 21930 21931 0 0 0
contig_1 21931 21932 0 0 0
contig_1 21932 21933 0 0 0
contig_1 21933 21934 0 0 0
contig_1 21934 21935 0 0 0
contig_1 21935 21936 0 0 0
contig_1 21936 21937 0 0 0
contig_1 21937 21938 0 0 0
contig_1 21938 21939 0 0 0
contig_1 21939 21940 0 0 0
contig_1 21940 21941 0 0 0
contig_1 21941 21942 0 0 0
contig_1 21942 21943 0 0 0
contig_1 21943 21944 13 22 105
contig_1 21944 21945 0 0 0
contig_1 21945 21946 0 0 0
contig_1 21946 21947 0 0 0
contig_1 21947 21948 0 0 0
contig_1 21948 21949 0 0 0
contig_1 21949 21950 0 0 0
contig_1 21950 21951 0 0 0
contig_1 21951 21952 0 0 0
contig_1 21952 21953 0 0 0
contig_1 21953 21954 0 0 0
contig_1 21954 21955 0 0 0
contig_1 21955 21956 0 0 0
contig_1 21956 21957 0 0 0
contig_1 21957 21958 0 0 0
contig_1 21958 21959 0 0 0
contig_1 21959 21960 0 0 0
contig_1 21960 21961 0 0 0
contig_1 21961 21962 13 6 134
contig_1 21962 21963 0 0 0
contig_1 21963 21964 0 0 0
contig_1 21964 21965 0 0 0

En utilisant bash j'aimerai créer un tableau qui stocke toute les 51 lignes la valeur obtenu du calcul suivant : somme des valeurs de la colonne 4 / (somme des valeur de la colonne 4 + somme des valeurs de la colonne 5 + somme des valeurs de la colonne 6).
J'aimerai aussi mettre une autre colonne en première position qui met le span (la position de départ et de fin) donc par example dans l'example que j'ai donnée ce serait de 21914-21964 et mettre en deuxième colonne la valeur du calcul.

Merci d'avance pour votre aide, je suis un peu nouvelle en bioinformatique et je ne m'y connais pas trop toute aide sera bienvenue
-----

[micapivi]

Bonjour

J'ai essayé de n'utiliser que les fonctionnalités de bash, mais j'ai quand même utilisé bc pour faire la division.

#!/bin/bash

[[ $# -eq 0 ]] && { echo "Donnez le nom du fichier à traiter"; exit; }

nblignes=51
nbrDecimales=16  # pour la division
cmptLignes=0
depart=0
fin=0
total4=0
total5=0
total6=0

while read champ1 champ2 champ3 champ4 champ5 champ6; do 
  (( cmptLignes += 1 ))
  (( depart == 0 )) && depart=$champ2
  (( total4 += champ4 ))
  (( total5 += champ5 ))
  (( total6 += champ6 ))
  fin=$champ2
  (( cmptLignes >= nblignes )) && {
    echo $depart $fin $( bc <<< "scale=$nbrDecimales; $total4 / ( $total4 + $total5 + $total6 )" )
    depart=0
    fin=0
    cmptLignes=0
    total4=0
    total5=0
    total6=0
  }
done  < "$1"

Ce qui donne, avec un fichier de 51 lignes minimum :

mic@deb12:~$ ./testScript.sh fichier
21914 21964 .0855855855855855
mic@deb12:~$

[polo974]

Bonjour,
En fait, bash ne sait pas faire d'opérations flottantes.

Pourquoi cette contrainte?
N'y a-t-il pas un langage plus adapté disponible sur la machine?

(on peut remplacer bc par plein d'autre programmes donnant accès aux flottants, mais autant en prendre un, genre python et tout faire avec, ce sera bien plus lisible et formateur)

[Henrix]

Quel est le nombre de lignes maximum de tes fichiers ?

[A voir en vidéo sur Futura]

[pm42]

Une ligne :

awk '{s4+=$4; s5+=$5; s6+=$6; if(NR%51==1) {start = $2 }; if(NR%51==0) { end=$3; printf("%f %s-%s\n", s4/(s4+s5+s6), start, end) } }' methylation_extracted_downstre am.bed

[micapivi]

avec le script bash que je t'ai proposé, la sortie se présentera comme ça :

21914 21964 .0855855855855855

mais si tu veux qu'un tiret sépare les deux champs comme ça :

21914-21964 .0855855855855855

il te faudra remplacer dans le script :

  echo $depart $fin

par :

  echo $depart-$fin

=======================
D'autre part, s'il y a une ou plusieurs espace dans le nom de ton fichier, comme par exemple :

methylation_extracted_downstre am.bed

dans la ligne de commande, il te faudra mettre ce nom de fichier entre des doubles guillements :

"methylation_extracted_downstre am.bed"

[pm42]

Une ligne :

awk '{s4+=$4; s5+=$5; s6+=$6; if(NR%51==1) {start = $2 }; if(NR%51==0) { end=$3; printf("%f %s-%s\n", s4/(s4+s5+s6), start, end) } }' methylation_extracted_downstre am.bed

Petite correction pour remettre les variables de somme à 0 toutes les 51 lignes.

awk '{s4+=$4; s5+=$5; s6+=$6; if(NR%51==1) {start = $2 }; if(NR%51==0) { end=$3; printf("%f %s-%s\n", s4/(s4+s5+s6), start, end); s4=0;s5=0;s6=0 } }' methylation_extracted_downstre am.bed

[polo974]

awk != bash

Ah, mais...

Allez, il y a quelqu'un pour le faire en cobol?

(C'est vachement recherché comme compétence...)

[pm42]

awk != bash

Ca se discute et dans ce cas, bc != bash

L'intérêt du shell, c'est quand même justement toutes les commandes faites pour être utilisées dedans.
Et là, c'est du simple : pas de | , de ``, xargs et autres.

Qui plus est, awk est justement conçu pour ça : séparer des colonnes, faire des calculs dessus, etc.

[Henrix]

Mes deux cent's (pour Windows et pour le fun)
J'ai ajouté des données pour tester
Remplacer les données dans le fichier: "données.txt"
Lien 14 jours: https://www.mediafire.com/file/jnwea...1+102.zip/file

[micapivi]

Une erreur dans mon script,
il faudrait remplacer :

  fin=$champ2

par :

  fin=$champ3

[Oceane666]

Merci beaucoup pour votre aide, j'ai utilisé la commande de pm42 qui marche parfaitement.

[micapivi]

Avec awk au lieu de bc on y perd un peu en précision,
mais bon, c'est peut-être pas si important que ça pour ce que tu comptes faire :

mic@deb12:~$ awk 'BEGIN{ printf "%.20f\n",1/3}'
0.33333333333333331483
mic@deb12:~$ 
mic@deb12:~$ echo 'scale=20; 1/3' | bc
.33333333333333333333
mic@deb12:~$

[polo974]

Avec awk au lieu de bc on y perd un peu en précision,
mais bon, c'est peut-être pas si important que ça pour ce que tu comptes faire :

Ben oui, bc est une calculette base 10 de précision arbitraire...

Mais awk (version gnu) fait (4 fois) mieux que du double, il faut juste bien lui parler (positionner la variable PREC).

Sinon en python, il y la lib decimal qui permet aussi de la précision arbitraire en base 10, mais il faut vraiment en avoir besoin vu le coût en terme de perf et que (sauf erreur) ça ferme la porte à numpy et du coup, sans doute aussi à pandas qui l'utilise.

[pm42]

On peut aussi préciser qu'il n'y a pas des masses de domaine en science où on a besoin de plus de précision ce qu'un float ou un double donne.

Un float c'est déjà la mesure d'une tonne au gramme prêt.
Les choses sont plus pointues quand on fait du calcul numérique plus massif, des matrices avec des valeurs extrêmes mais dans l'exemple donné vu qu'on calcule des moyennes de 51 nombres entiers pas très grands, la question ne sd pose même pas.