Réponse aux exercices

BASH

boucles, conditions, variables

  • Générez et affichez 100 valeurs aléatoires comprises entre 0 et 10

    for i in `seq 1 100`; do
      ((num=RANDOM%11))
      echo $num
    done
    
  • Créez un tableau de 100 valeurs aléatoires comprises entre 0 et 10.

    tab=()
    for i in `seq 1 100`; do
      ((num=RANDOM%11))
      tab[$i]=$num
    done
    
  • Créez un tableau qui contient les valeurs du tableau précédent préfixées par le texte "pair" si la valeur est paire
    tabpair=()
    for i in `seq 1 100`; do
      if ((tab[i] % 2 == 0)); then
          tabpair[$i]="pair "${tab[$i]}
      else
          tabpair[$i]=${tab[$i]}
      fi
    done
    echo "contenu du tableau : "${tabpair[@]}
    
  • Faites une boucle qui compte le nombre de fois qu'on a un nombre pair sur 100 000 nombres aléatoires. `` bash compteur=0 for i inseq 1 100000`; do if ((RANDOM % 2 == 0)); then
      ((compteur++))
    
    fi done

echo "le compteur vaut $compteur"


* Créer cent dossiers resu1 à resu100 dans le home.

``` bash
 for i in `seq 1 10`
 do
   mkdir resu$i
 done
  • Créer cent dossier resu0, resu2, resu4 et ainsi de suite avec tous les nombres pairs.
# plusieurs solutions
# en incrémentant un compteur de deux
compteur=0
for i in `seq 1 100`; do
    mkdir resu$compteur
    ((compteur+=2))
done

# ou encore en multipliant la variable de boucle par deux
for i in `seq 0 99`; do
    ((num=i*2))
    mkdir resu$num
done

# ou même en utilisant juste seq correctement
for i in `seq 0 2 198`; do
    mkdir resu$i
done
  • Ecrire un script qui prend un paramètre et crée 10 dossiers dont le nom est préfixé par ce paramètre.
#!/bin/bash

for i in `seq 1 10`; do
    mkdir ${1}${i}
done
  • exercice supplémentaire

  • Générer une séquence aléatoire de 100 bases avec une probabilité de 70% pour la base 'A' et 10% pour chacune des autres. Indication : voila un if avec un ET entre deux conditions :

if ((2 < n)) && ((n < 30)); then
s=""
for i in `seq 100`; do
    n=$[RANDOM % 100]
    echo $n
    if (( n < 70 )); then
        b=A
        echo "$n < 70"
    else if (( 80 > n)) && ((n > 70)); then
        b=T
    else if (( 90 > n )) && ((n> 80)); then
        b=C
    else if (( 100 > n)) && ((n > 90)); then
        b=G
    fi; fi;fi;fi

    s="$s$b"
done

echo $s

Le PIPE

  • Compter le nombre de lignes dans le fichier test.fastq REPONSE
wc -l test.fastq

ou bien avec un PIPE :

cat test.fastq | wc -l
  • Compter le nombre de séquences dans le fichier test.fastq
grep "^@" test.fastq | wc -l
  • Faire la même chose et stocker le résultat dans un fichier resu.txt
grep "^@" test.fastq | wc -l > resu.txt
  • Extraire uniquement les séquences (sans les identifiants) du fichier omm.fasta (utilisez grep -v)
grep -v "^>" omm.fasta
  • Quelle est la valeur max de MAPQ (5° colonne dans le fichier resu.sam) ?
cut -f 5 resu.sam | sort -n | tail -n 1

# je prend la 5e colonne | je trie les valeurs par ordre croissant | je prend la dernière donc la plus grande
  • Compter le nombre de séquences ne contenant pas le mot "rand" dans le fichier test.fastq
grep "^@" test.fastq | grep -v rand | wc -l
  • Compter le nombre de fichiers dans votre dossier personnel (utilisez ls -A pour voir tous les fichiers)
ls -A | wc -l
  • Compter le nombre de fichiers dont le nom commence par "." dans votre dossier personnel (indice : l'expression à donner à grep est : "^.")
ls -A | grep "^\." | wc -l

Fichiers

REPONSE :

#!/bin/bash

mkdir exercice1
cd exercice1
wget http://mbb.univ-montp2.fr/MBB/uploads/formationBash.tar.gz
mkdir data
mv formationBash.tar.gz data/
cd data
tar xvf formationBash.tar.gz
chmod go-rwx ./*
cd ../..
  • Comparer le taux de compression du fichier test.fastq à l'aide des commandes zip et gzip

REPONSE :

zip test.fastq.zip test.fastq
gzip test.fastq
du -s   test.fastq.zip   test.fastq.gz
  • Extraire les séquences dont l'identifiant contient "rand" du fichier Gzippé sans le décompresser (zcat)

REPONSE :

zcat test.fastq.gz | grep "rand" > rand_seq.txt

Outils système

  • Lancez un "sleep 15000" et essayez de le tuer avec un deuxième terminal en une ligne en utilisant les commandes ps, grep, kill et awk. (indice : attention au processus produit par grep dont le nom contient sleep lui aussi)

Avec grep, on filtre la sortie de ps, ensuite on enleve la ligne qui contient grep puis on prend juste la 2e colonne :

kill -9 `ps aux | grep sleep | grep -v grep | awk '{print $2}'`
  • Lancez dix fois la commande précédente (sleep), avec 300 secondes au lieu 15000, en arrière plan :
for i in `seq 1 10`; do (sleep 300 &) ; done

et vérifiez si votre kill précédent fonctionne aussi.

REPONSE Oui il marche encore parce que le ps dans grep dans awk va renvoyer tous les PIDs correspondants aux sleep. Tout ça va être donné en paramètre à kill -9 qui est capable de traiter plusieurs PID en paramètre d'après son manuel.

  • La commande "ps -ef" affiche en troisième colonne le PPID (PID du parent). A l'aide de cette information, trouvez le nom du processus parent de firefox.

REPONSE : Evidemment, il faut que firefox tourne. Avec ps -ef et grep et awk, on récupère le PID du parent de firefox en 3e colonne. Puis on cherche le nom de ce PID avec ps aux.

ps aux | grep `ps -ef | grep firefox | grep -v grep | awk '{print $3}'` | awk '{print $11}'
  • Trouvez quel est le processus qui vous prend le plus de mémoire vive avec ps, sort et head.

REPONSE Il faut trier numeriquement la 4e colonne de la sortie de ps aux

ps aux | sort -r -n -k 4,4 | head

Si vous voulez juste les noms :

ps aux | sort -r -n -k 4,4 | head | awk '{print $11}'

Cluster

SGE en détail

  • Soumettre via SGE directement un appel à calc_mean.pl en prenant comme entrée un des fichiers de nombre aléatoires précédemment généré.
qsub -N direct -cwd -b y "calc_mean.pl  numbers.txt  mean.txt"
  • Soumettre le script SGE calc_mean.sge . Il prend les mêmes paramètres que calc_mean.pl
qsub calc_mean.sge  numbers.txt  mean.txt
  • Modifiez calc_mean.sge pour utiliser la variable $TMPDIR et faire en sorte que le résultat de calc_mean.pl soit écrit dans $TMPDIR/ et ensuite copié dans votre HOME.
#!/bin/bash
#$ -cwd
#$ -j y
#$ -S /bin/bash
#$ -N moy
#$ -V

# on ecrit dans TMPDIR
./calc_mean.pl $1 $TMPDIR/$2
# on recupère le résultat dans le dossier courant (celui duquel on a lancé le job)
cp $TMPDIR/$2 ./

Exercices supplémentaires :

  • lancement d'un job multicoeur en déclarant le nombre de slots requis
#!/bin/bash
#$ -S /bin/bash
#$ -cwd
#$ -sync y
#$ -N multicore_job
#$ -q cemeb.q
#$ -pe multithread60 2

mon_programme  max_threads=$NSLOTS 2>&1
  • utiliser un environnement multicore pour soumettre un job qui demande un nombre variable de slots et faire afficher le nombre de slots effectivement alloués par SGE (variable $NSLOTS)
#!/bin/bash
#$ -S /bin/bash
#$ -cwd
#$ -sync y
#$ -N multicore_job
#$ -q cemeb.q
#$ -pe multithread60 4-15

mon_programme  max_threads=$NSLOTS 2>&1
  • extraire les iris dont le sépale est supérieur à 1,4cm
while read line; do.
    petal=`echo $line | cut -d ',' -f 3`
    superior=`echo "$petal > 5.4" | bc -l`
    if [[ $superior == 1 ]]; then
        echo $line
    fi
done <$1

Analyse de texte

grep

  • Chercher toutes les séquences contenant "Loxodonta" dans omm.fasta

REPONSE

grep "Loxodonta" omm.fasta
  • Filtrer uniquement les séquences Loxodonta du gène NUDT3

REPONSE Deux solutions : La première marche quel que soit l'ordre des mots recherchés :

grep "Loxodonta" omm.fasta | grep NUDT3

Et si on sait que Loxodonta va apparaitre avant NUDT3 :

grep "Loxodonta.*NUDT3" omm.fasta
  • Exclure toutes les séquences de l'homme (HOMO)

REPONSE : Pour être sur qu'on prend en compte toutes les façon d'écrire homo, on utilise l'option "-i" de grep. Attention on exclue les lignes qui contiennent homo donc on n'exclue pas le contenu des séquences. Il ne faut pas oublier de filtrer sur le ">" en début de ligne :

grep -v -i "homo" omm.fasta | grep "^>"
  • Afficher uniquement les identfiants et les sequences du fichier test.fastq (Dans ce fichier les identifiants commencent par @, les séquences sont formées ici par A, C, G, ou T et les autres lignes contiennent que des nombres)
grep "^@" test.fastq | grep "^[ATGC]"

outils divers

Soit le texte suivant dans un fichier :

ind1,20
ind2,20
ind3,10
ind4,20
ind2,30
ind1,10
ind4,11
ind3,20
  • Combien de valeur différentes existe-t-il dans la première colonne ? dans la deuxième ?

REPONSE :

cat fichier.txt | cut -d ',' -f 1 | sort | uniq | wc -l
cat fichier.txt | cut -d ',' -f 2 | sort | uniq | wc -l
  • Triez les valeurs de la deuxième colonne par fréquence d'apparition. (uniq -c peut être utile)

REPONSE : Cette solution va marcher mais on a tout de même des espaces en début de ligne à la sortie de uniq :

cat fichier.txt | cut -d ',' -f 2 | sort | uniq -c | sort -n

Si on connait déjà sed, on peut supprimer les suites d'espaces ou tabulations en début de ligne :

cat fichier.txt | cut -d ',' -f 2 | sort | uniq -c | sed -r 's/^\s+//' | sort -k 1,1

SED

REPONSE

sed -i 's/,/\t/g' iris.data
  • Extraire de ce fichier la colonne 3 (petal length in cm) (utiliser cut )

REPONSE : Pour CUT, le séparateur par defaut est TAB. On peut donc faire :

cut -f 3 iris.data
  • Compter le nombre d'espèces (cut et uniq)

REPONSE : L'option "-s" de cut est très importante, elle n'affiche pas les lignes où le séparateur n'a PAS été trouvé. Une solution parmis d'autres est de couper sur les tirets (-) et de prendre la deuxième colonne :

cut -s -d '-' -f 2 iris.data | sort | uniq | wc -l

Une autre solution est d'utiliser SED pour supprimer tout ce qui précède le tiret (inclus). On doit éliminer la ligne vide avec grep :

sed 's/^.*-//g' iris.data | sort | uniq | grep -v "^$" | wc -l
  • Trier par longueur de pétale croissante (regardez les options de sort)

REPONSE : Ça marche dans ce cas là parce que l'ordre alphabétique est le même que l'ordre numérique :

cat iris.data | sort -k 3,3

Si on veut préciser le type de tri à sort, "-n" pour tri numérique :

cat iris.data | sort -n -k 3,3

ATTENTION le tri numérique de sort ne comprend pas les nombres dont la virgule est un POINT. Si vous utilisez un tri numérique , assurez vous

  • Changer le nom du taxon du format genre-espèce en espèce-genre

REPONSE : Ici on doit utiliser les parenthèses des expressions régulières pour faire de la répétition. La forme de ce qu'on veut remplacer est, une suite non nulle de lettres (majuscules ou minuscules) suivie d'un tiret suivi d'une autre suite de lettres. Si on met les parenthèses autour de chaque suite de lettres, on peut les réutiliser dans la partie droite du remplacement avec \1 et \2 :

sed -r 's/([a-zA-Z]+)-([a-zA-Z]+)/\2-\1/' iris.data
  • Insérer une ligne d'entête précisant les noms des champs :
  • sepal length in cm, 2. sepal width in cm, 3. petal length in cm, 4. petal width in cm, 5. species

REPONSE : Sed possède une commande d'insertion à laquelle on peut préciser la ligne :

sed '1iTEXT' iris.data

En remplaçant TEXT par le contenu à insérer. Si vous voulez appliquer la modification au fichier, utilisez l'option "-i" de sed.

  • Changer toute suite de N en un seul gap (-) dans le fichier omm.fasta

REPONSE :

On pourrait faire comme ça mais on risque de modifier le contenu des autres lignes.

sed -r 's/N+/-/g' omm.fasta

Pour être sur qu'on ne modifie que les séquences, on dit à sed de ne faire le remplacement que sur les lignes qui commencent par une base

sed -r '/^[ATGCN]/     s/N+/-/g' omm.fasta

AWK

  • Afficher uniquement les lignes dont la longueur des pétales est supérieure à la taille moyenne

UNE REPONSE POSSIBLE :

awk -F ',' '{
                ligne[NR]=$0;
                col3[NR]=$3;
                sum+=$3
            }
         END{
                moy = sum / NR ;
                for(i=1;i<=NR;i++){
                     if (col3[i] > moy){
                         print ligne[i]
                     }
                }
             }' iris.data

SECONDE SOLUTION :

awk -F"," '{
       # NR devient > à FNR quand on commence à lire le second fichier
       if (NR<=FNR){
            tot++;
            somme = somme + $3;
       }
       # ici on est dans la seconde copie du fichier 
       else{
            # on filtre sur les valeurs > à la moyenne
            if ($3 > somme/tot) print
       }
} ' iris.data iris.data

REPONSE :

# si le numero de ligne est un multiple de 4 plus 1 ou plus 2, on affiche la ligne
cat test.fastq | awk '{if (NR%4 == 1 || NR%4 == 2) {print $0 ;} }' | sed 's/^@/>/' > test.fasta
# on peut aussi procéder comme ça (getline force la lecture d'une ligne supplémentaire) :
# on affiche deux lignes consécutives puis on saute deux lignes
cat test.fastq | awk '{print; getline; print; getline;getline }' | sed 's/^@/>/' > 2

En pur awk :

awk '{
    if (NR%4 == 1){
        gsub('@','>',$0)
        print
    }
    if (NR%4 == 2)
        print
}' test.fastq > test.fasta

Rappel sur le format sam produit par bowtie2

Field What You Find There
1 readname
2 bitflag -- for SE bt2 mapping, only 0 (pos strand), 4 (unmapped), 16 (neg strand) are relevant
3 chrName/refname
4 1-based start on fwd strand (this is the end if on neg strand)
5 MAPQ
6 CIGAR string rep of alignment (most are 50M which means 50 match+mismatch)
7 chrName/refName of mate if PE; * if SE
8 1 based start of mate; * if SE
9 Inferred fragment length if PE; * if SE
10 the read's sequence
11 the reads base call qualities (Qs)
  • Compter le nombre de reads couvrant chaque contig du fichier resu.sam

REPONSE On parcours le fichier en comptant puis dans un bloc END, on affiche le résultat du compte

awk '{nb[$1]++}
     END{
         for (contig in nb){
             print nb[contig], contig
         }
     }' resu.sam
  • filtrer uniquement les mapping dont la MAPQ est >= 20

REPONSE :

awk '{if ($5 >= 20) print}' resu.sam
  • Calculer la qualité moyenne du mapping sur chaque contig

REPONSE :

awk '{
        # on fait la somme des qualités pour chaque contig
        sumq[$1] = sumq[$1] + $5;
        # on compte le nombre d'apparition de chaque contig
        nbq[$1]++;
     }
     END{
         # pour tous les contigs, on affiche son nom et la moyenne des qualités
         for (contig in sumq){
             print contig, sumq[contig]/nbq[contig]
         }
     }' resu.sam
  • Faites vous même une simplification du fichier fastq en remplaçant les identifiants des séquences par leurs rangs dans le fichier (@blabla en @22). Quel gain de place obtenez-vous ?

REPONSE

awk 'BEGIN{num=1}
     {
         if ($1 ~ "^@"){
             $1="@"num; 
             num++
         }
         print
     }' test.fastq > compressed.fastq

# qui peut être raccourci en n'initialisant pas la variable num et en l'incrémentant au moment de son utilisation :
awk '{
         if ($1 ~ "^@"){
             $1="@"++num
         }
         print
     }' test.fastq > compressed.fastq
  • Trouvez un moyen de garder les identifiants dans un fichier à part en conservant l'association entre chaque nombre et l'identifiant de contig.

REPONSE

awk '{
         # on ne traite que les lignes qui commencent par "@"
         if ($1 ~ "^@"){
             # on affiche le numero et l identifiant correspondant
             print "@"++num, $1
         }

     }' test.fastq > identifiants.txt

Ensuite on pourra reconstruire le fastq original (avec les noms complets) à partir du fichier compressé compressed.fastq et du fichier d'identifiants identifiants.txt de cette manière :

awk '{
         # on parcours les noms pour construire le tableau de correspondance
         if (NR<=FNR){
             name[$1] = $2
         }
         # on parcours compressed
         else{
             # pour les lignes qui commencent par "@", on remplace 
             if ($1 ~ "^@"){
                 $1 = name[$1]
             }
             print
         }
     }' identifiants.txt compressed.fastq > uncompressed.fastq