[Resolu] gestion format fasta

valnovice · Le 16/01/2017, à 17:25

OK, j'ai bien tout suivi mais ca n'a pas marché pour les premiers. Ce dernier a l'air de marcher (on ne crie pas encore victoire!!) ca tourne encore (ca fait presque 2h). Je reviens après pour te donner des nouvelles

kholo · Le 16/01/2017, à 18:07

ah non, là t'as rien raté du tout...
c'est même plutôt l'inverse !
il faudra quand même que j'apprenne à optimiser mon code... ou à coder déjà !
merci pour la démonstration Watael

PS : attention valnovice, Watael utilise "gene.fasta" et moi "genes.fasta "

valnovice · Le 16/01/2017, à 18:20

OK merci ca tourne encore ...
Watanael : j'ai essayé mais je n'ai pas encore réussi à avoir le résultat escompté avec ce que t'as donné mais je relis tout et refais tout . j'ai perdu tous les entêtes ... en tout cas, merci beaucoup

Watael · Le 16/01/2017, à 18:35

peut-on trouver tes fichiers quelque part (ou des similaires, pour qu'on fasse des tests grandeur nature) ?
correspondent-ils exactement à l'exemple que tu as fourni ?

comment exécutes-tu les instructions que j'ai données ?

Watael · Le 16/01/2017, à 21:15

je ne sais pas.
2500 séquences, ça représente combien de lignes ? combien de Ko, Mo... ?

quelqu'un connaît-il un site d'hébergement pour gros fichiers, dans le style de pastebin (dont la limite est de « 512 kilobytes (0.5 megabytes) »)?

Dernière modification par Watael (Le 16/01/2017, à 21:21)

Watael · Le 16/01/2017, à 23:55

ah, oui, le format ne change pas, mais le fait qu'il y ait des doublons fait échouer le code que j'ai proposé.

en voilà un autre qui semble faire ce qu'il faut :

fichier=~/Téléchargements/exemple_prot.fasta
gawk 'BEGIN{PROCINFO["sorted_in"]="@ind_str_asc"}/^>/{Gene=$2; if ( $2 in compSeq); else compSeq[$2]=$0} !/^>/{Seq[Gene][cptSeq[Gene]++]=$0}END{ for (i in Seq){ print compSeq[i]; for (y in Seq[i])print Seq[i][y]}}' "$fichier"

on oublie le shell, et on passe à GNU/awk pour pouvoir utiliser des tableaux de tableaux, et des fonctions de tri que ne propose pas (en tout cas, je ne m'en souviens pas) KSH.

mais, comment tu vas faire pour vérifier que la sortie est correcte ?

valnovice · Le 17/01/2017, à 00:04

j'ai utilisé ca pour enlever les doublons : sed -e '/^>/s/$/@/' -e 's/^>/#/' file.fasta | tr -d '\n' | tr "#" "\n" | tr "@" "\t" | sort -u -t $'\t' -f -k 2,2 | sed -e 's/^/>/' -e 's/\t/\n/' > file.fasta
je travaille sur plusieurs fichiers et celui du lien est un exemple concret et si marche pour cet exemple, ça devrait marcher pour les miens. Pour la sortie, je ne sais absolument rien.

Watael · Le 17/01/2017, à 00:13

sed -e '/^>/s/$/@/' -e 's/^>/#/' file.fasta | tr -d '\n' | tr "#" "\n" | tr "@" "\t" | sort -u -t $'\t' -f -k 2,2  | sed -e 's/^/>/' -e 's/\t/\n/' > file.fasta

il doit y avoir une erreur, car là, tu as écrasé le contenu de file.fasta.
file.fasta devrait être vide finalement, parce que la redirection est effectuée en premier.

quand je parle de doublons, je parle de genes qui portent le même nom sur des lignes différentes : par exemple, il y a 32 lignes pour [gene=atp6].

valnovice · Le 17/01/2017, à 00:27

oui, c'est vrai, là j'ai écrasé.
je re-explique :
il y a à peu près 184 individus (identifiants) dans ce fichier, les identifiants commencent par : >lcl|JN555591.1 ou >lcl|NC_031377.1 ou >lcl|NC_004454.2 ou >lcl|HM189212.1 ...
et pour cet exemple chaque individu (identifiant) possède 13 gènes (gene=ATP6, gene=ATP8, gene=COX1...). Et là tu me fais peur car si tu me dis qu'il n'y a que 32 lignes qui ont [gene=atp6] alors les 184-36=148 n'ont pas le gène atp6??? chacun des identifiants devrait avoir un gene ATP6

Dernière modification par valnovice (Le 17/01/2017, à 00:29)

Watael · Le 17/01/2017, à 01:12

 $ grep '^>lcl' $fichier | wc -l
2212

 $ awk 'BEGIN{PROCINFO["sorted_in"]="@ind_str_asc"}/^>lcl/{ar[$2]++}END{for (i in ar){print i": "ar[i];n++}print "\n"n": individus"}' $fichier 
[gene=ATP6]: 139
[gene=ATP8]: 139
[gene=CO1]: 7
[gene=CO2]: 7
[gene=CO3]: 7
[gene=COIII]: 10
[gene=COII]: 10
[gene=COI]: 19
[gene=COX1]: 116
[gene=COX2]: 116
[gene=COX3]: 116
[gene=COXIII]: 5
[gene=COXII]: 6
[gene=COXI]: 6
[gene=CYB]: 3
[gene=CYTB]: 131
[gene=CytB]: 2
[gene=NAD1]: 19
[gene=NAD2]: 20
[gene=NAD3]: 20
[gene=NAD4L]: 19
[gene=NAD4]: 19
[gene=NAD5]: 19
[gene=NAD6]: 19
[gene=ND1]: 119
[gene=ND2]: 117
[gene=ND3]: 115
[gene=ND4L]: 119
[gene=ND4]: 118
[gene=ND5]: 119
[gene=ND6]: 114
[gene=atp6]: 32
[gene=atp8]: 32
[gene=cob]: 27
[gene=cox1]: 32
[gene=cox2]: 32
[gene=cox3]: 32
[gene=cytB]: 1
[gene=cytb]: 5
[gene=nad1]: 32
[gene=nad2]: 32
[gene=nad3]: 32
[gene=nad4L]: 30
[gene=nad4]: 32
[gene=nad4l]: 2
[gene=nad5]: 32
[gene=nad6]: 32

47: individus

Dernière modification par Watael (Le 17/01/2017, à 01:19)

valnovice · Le 17/01/2017, à 08:41

Watael : Ok j'ai compris. [gene=ATP6] n= [gene=atp6]. ce n'est pas uniforme la nomenclature. Faudra que j'uniformise. merci
Kholo : ça a tourné pendant 11h et le résultat : il n'y a plus de chevron donc pas d'identification d'aucune séquence, et je partais d'un fichier de 900 Ko pour arriver à 18 Mo (c'est probablement normal)

kholo · Le 17/01/2017, à 11:20

je reprend la démo de Watael en plus lisible
(enfin surtout pour moi, parce que j'ai pris une bonne claque là !) :

fichier="genes.fasta"
read header < "$fichier"; 
echo "$header"; 
while read info gene reste; 
do 
	[[ $info == '>'* ]] && { Gene=$gene; n=0;} || echo "$Gene$((n++))=$info"; 
done < "$fichier"| sort -k1,1 | cut -d '=' -f3

ça affiche tout directement dans le terminal donc il reste à le diriger dans un fichier

adapté à mon script :

#!/bin/bash
# ----------------------------------------------
# nomlogiciel="$(basename "$0")" # on peut récupérer le seul nom du script par Développement des paramètres, 
# comme tu le fais plus loin pour l'extension
nomlogiciel="${0##*/}"

FONCTION="Concatener des séquences de gènes d'une espèces et ne conserve que le premier commentaire"
VERSION="alpha"
# NOTES DE VERSIONS
# https://forum.ubuntu-fr.org/viewtopic.php?id=2003118
# merci Watael pour les conseils :
# https://forum.ubuntu-fr.org/viewtopic.php?pid=21662278#p21662278
# ----------------------------------------------
# à mettre au début d'un fichier bash
#PID=$$
#FIFO=/tmp/FIFO${PID}
#mkfifo ${FIFO}
# ----------------------------------------------
echo "lancement $nomlogiciel..." ; 

#-----------
# Variables
#-----------

# $1 : argument d'entrée du script
# $ft : fichier traité
# $ext : variable correspondant à l'extension du fichier entré en argument

# Pour débugage
if [ -n "$1" ]; then
	le_fichier="$1"
else # dans le cas de test simple on peut sans souci utiliser des ET et OU logique plutôt qu'un lourd if
	# pointer vers le fichier de test
	le_fichier="genes.fasta"
fi
# echo "le_fichier=$le_fichier"

#-------------------------
# Traitement des fichiers
#-------------------------

# Récupération de l'extension
ext="${le_fichier##*.}"
# echo "ext=$ext"
[[ $ext == fas?(ta) ]] || { echo -e "Le fichier n'est pas un ficher fasta. Sortie."; exit;}

# Création de la variable correspondant au fichier traité
ft="${le_fichier%.${ext}}_traité.${ext}"

# lecture du fichier à traité
# listSequences=$(grep '>' "$le_fichier")
# nbreSequences=$(echo "$listSequences" | grep -c ">")
# echo "nbre de sequences : $nbreSequences
# Séquences :
# $listSequences"

# clear

# reprend l'entête du fichier d'origine
# et initialise le fichier de sortie avec  
read header < "$le_fichier"; # header="$(head -1 "$le_fichier")" # celui de Watael
# on peut traiter l'entête ici
# ...
echo "$header" | tee "$ft"

# traitement et envoie dans le fichier
echo "$(
while read info gene reste; 
do 
	[[ $info == '>'* ]] && { Gene=$gene; n=0;} || echo "$Gene$((n++))=$info"; 
done < "$le_fichier"| sort -k1,1 | cut -d '=' -f3
)" | tee -a "$ft"

# Pour voir le résultat dans gedit après traitement
# décommenter la prochaine ligne
# nano "$ft"
# gedit "$ft"

exit 0

Watael · Le 17/01/2017, à 11:53

pourquoi

echo "$(while read...)"

?

kholo · Le 17/01/2017, à 12:46

plus propre comme ça ?

while read info gene reste; 
do 
	[[ $info == '>'* ]] && { Gene=$gene; n=0;} || echo "$Gene$((n++))=$info"; 
done < "$le_fichier"| sort -k1,1 | cut -d '=' -f3 | tee -a "$ft"

Watael · Le 17/01/2017, à 12:56

oui, pas de sous-process induit par la substitution de commande.

mais je ne suis pas convaincu que ça fonctionne correctement : 11 heures c'est énorme.
la ligne awk donnée en #33 traite le fichier de 900K en quelques secondes.

valnovice · Le 17/01/2017, à 16:20

kholo : j'ai testé, ça donne quelques chose mais il n'y a plus de ">" nulle part et difficile de distinguer qui est qui
watael : avec #33, effectivement, ca marche en quelques secondes et il y a des chevrons ">" (ils sont 47, j'utilise le fichier que j'ai partagé pour m'entrainer) mais tous les mêmes gènes sont "groupés dans un >" (mais ce n'est plus l'identification de l'individu) on a perdu les identifiants. et je ne sais pas à quoi correspond le n=47
en tout cas, merci beaucoup à tous les 2

Dernière modification par valnovice (Le 17/01/2017, à 16:23)

Watael · Le 17/01/2017, à 17:03

tu veux dire que la sortie, pour le gène nad6, par exemple devrait être :

>lcl|KJ957822.1_prot_AIH15205.1_8 [gene=nad6] [protein=NADH dehydrogenase subunit 6] [protein_id=AIH15205.1] [location=complement(7625..8071)]
MFLSVTLFVLFLMLESSSSPFKISLMLAFVSTMCSAFLYFKSSSIFLCCSLVISFSSGMMILFSYCSMFS
SCEYKNKAKLNLSPLILSLLMLSMLMVPQTTLFSSMEKLCSMVSVCLYMVFLMFIVIISMSAMNISFFNP
NESMNQKVNKGTMVMVITMITISMKMKPQDQMNMINNTLFMKINNSIMMIMMIIIMTMTFASKMSMNPFS
PSKSIK
lcl|JXLN01000001.1_prot_KPL93350.1_4 [gene=nad6] [protein=NADH dehydrogenase subunit 6] [protein_id=KPL93350.1] [location=3368..3682]
MIMMFFFIFNFSPLKFSFFLVFFSVLFSFFVYMFSLSLFLMSVVIISYSTGMMIIYSYCSLFSVHEEKKS
MVKYNMYMYXMYMYYYFLYFFLLFGFYKLFYTSN
lcl|KF197136.1_prot_AGW06994.1_11 [gene=nad6] [protein=NADH dehydrogenase subunit 6] [protein_id=AGW06994.1] [location=9402..9830]
MKLMTFFTMIFFIFNHPMYMLMSIIFITLMMSILIYKSMKMMWISLILILLILGGMLILFLYIISLIPNK
KLSLNKKIMMLMIFLLMSFPPFKMLETSFLLMNYNFFTSSMNMMLFMMIYLICTLIVVMNIMTSSNAPLS
TL
...

?
il faut regrouper les gènes, dans l'ordre d'apparition dans le fichier original (edit: ah,non, c'est vrai, il faut les trier par ordre alphabétique), conserver les identifiants en préservant intact le premier et en supprimant le premier caractère ( '>') des suivants ?

mais, une petite question : comment tu vas faire croire que c'est toi qui a pondu nos codes ?

Dernière modification par Watael (Le 17/01/2017, à 17:12)

kholo · Le 17/01/2017, à 17:05

bon, il faut reprendre pour être à la même vitesse...
avec

awk 'BEGIN{PROCINFO["sorted_in"]="@ind_str_asc"}/^>lcl/{ar[$2]++}END\
{for (i in ar){print i}}' \
$le_fichier

tu récupères tous les gènes ensuite tu pourras les passer un par un en maj ou min
pour le moment :

[gene=ATP6]
[gene=ATP8]
[gene=CO1]
[gene=CO2]
[gene=CO3]
[gene=COIII]
[gene=COII]
[gene=COI]
[gene=COX1]
[gene=COX2]
[gene=COX3]
[gene=COXIII]
[gene=COXII]
[gene=COXI]
[gene=CYB]
[gene=CYTB]
[gene=CytB]
[gene=NAD1]
[gene=NAD2]
[gene=NAD3]
[gene=NAD4L]
[gene=NAD4]
[gene=NAD5]
[gene=NAD6]
[gene=ND1]
[gene=ND2]
[gene=ND3]
[gene=ND4L]
[gene=ND4]
[gene=ND5]
[gene=ND6]
[gene=atp6]
[gene=atp8]
[gene=cob]
[gene=cox1]
[gene=cox2]
[gene=cox3]
[gene=cytB]
[gene=cytb]
[gene=nad1]
[gene=nad2]
[gene=nad3]
[gene=nad4L]
[gene=nad4]
[gene=nad4l]
[gene=nad5]
[gene=nad6]

ensuite,

awk 'BEGIN{PROCINFO["sorted_in"]="@ind_str_asc"}/^>lcl/{ar[$1]++}END\
{for (i in ar){print i}}' \
$le_fichier

on obtient la partie 1 des com :

...
>lcl|NC_032369.1_prot_YP_009332018.1_9 
>lcl|NC_032369.1_prot_YP_009332019.1_10 
>lcl|NC_032369.1_prot_YP_009332020.1_11 
>lcl|NC_032369.1_prot_YP_009332021.1_12 
>lcl|NC_032369.1_prot_YP_009332022.1_13 
# je sépare pour y voir plus clair
>lcl|NC_032381.1_prot_YP_009332308.1_1 
>lcl|NC_032381.1_prot_YP_009332309.1_2 
>lcl|NC_032381.1_prot_YP_009332310.1_3 
>lcl|NC_032381.1_prot_YP_009332311.1_4 
>lcl|NC_032381.1_prot_YP_009332312.1_5 
>lcl|NC_032381.1_prot_YP_009332313.1_6 
>lcl|NC_032381.1_prot_YP_009332314.1_7 
>lcl|NC_032381.1_prot_YP_009332315.1_8 
>lcl|NC_032381.1_prot_YP_009332316.1_9 
>lcl|NC_032381.1_prot_YP_009332317.1_10 
>lcl|NC_032381.1_prot_YP_009332318.1_11 
>lcl|NC_032381.1_prot_YP_009332319.1_12 
>lcl|NC_032381.1_prot_YP_009332320.1_13

à la fin "032381.1" est un individu et "0093323xx.x" doivent se rapporter à ses gènes

donc la vraie demande serait
prendre un individu, organiser ses gènes par ordre alphabétique, sortir une seule entête (réécrite ou adaptée), puis la concaténation de ses gènes et ainsi de suite pour les autres individus...
ça colle ?

le problème reste de ne plus avoir de séparation entre les gènes de chaque individus,
donc ne vaut il pas mieux se 'contenter' de mettre les gènes de chaque individus dans l'ordre alphabétique du nom des gènes.... sauf que là, on va avoir les réf des gènes dans le désordre...

Watael · Le 17/01/2017, à 17:14

est-ce qu'on pourrait avoir un exemple de ce à quoi doit ressembler la sortie pour deux ou trois gènes donnés ?

pingouinux · Le 17/01/2017, à 17:48

Watael #46 a écrit :

est-ce qu'on pourrait avoir un exemple de ce à quoi doit ressembler la sortie pour deux ou trois gènes donnés ?

C'est en gros ce que je demandais en #2…

valnovice · Le 17/01/2017, à 17:55

Watael : merci pour ta remarque mais je n'ai pas besoin de faire croire à quiconque que c'est moi qui ai pondu VOS codes, je ne m'en vanterai pas, je ne suis plus étudiante, c'est un travail personnel et cette étape de nettoyage de fichier n'est que le début d'une très très longue histoire. Je ne serais pas ici si je savais comment faire.
Kholo : merci !!
Le résultat attendu ressemble à ce que kholo a fait dans #11 mais là, les gènes étaient encore en désordre.
>lcl|HM189212.1_prot_ADJ68021.1_9 --> avec les gènes concatenés en Ordre (par ordre alphabétique puis numérique) derrière cet identifiant (ATP6, ATP8, COX1, COX2 COX3....)
>lcl|NC_004454.2_prot_ND2_16830_1 --> idem
>lcl|KX452726.1_prot_AOR08478.1_8 --> idem

Il se pourrait que certains identifiants n'ont pas tous les gènes mais ça ce sera une autre histoire, et disons que ce n'est pas grave

valnovice · Le 17/01/2017, à 18:10

pingouinux : vraiment desolee c'est vrai que je ne t'ai pas répondu. et là je pense que j'ai répondu ta question?!

kholo · Le 17/01/2017, à 18:27

donc
étape 1 : normaliser le nom des gènes (tout en majuscule) pour avoir
[gene=XXXXX] où XXXXX est le nom du gène en maj

étape 2 : séparer les individus ;

>lcl|NC_031377.1_prot_YP_...

ici l'individu serait NC_031377.1
peut être mettre chaque individu dans un fichier séparé ?

étape 3 : reprendre les gènes de chaque individu indépendamment et les ordonner par ordre alphab.
ça c'est ce qu'on a fait avant...

valnovice · Le 17/01/2017, à 18:32

kholo : exactement! dans cet exemple, il faudra faire une normalisation. Oui, on pourrait mettre chaque individu dans un fichier séparé mais je pense qu'il y a 184 individus !!

Dernière modification par valnovice (Le 17/01/2017, à 18:34)

pingouinux · Le 17/01/2017, à 18:37

Donc à l'étape 3, tu perds toutes les informations qui étaient sur les lignes (sauf la première) commençant par >. Notamment, les noms des gènes ont disparu (sauf le premier). Est-ce bien ce que tu veux ?

Remarque : Merci d'utiliser les balises-code <> pour renvoyer les informations volumineuses.

Édité : Je parle de l'étape 3 de ton message #49

Dernière modification par pingouinux (Le 17/01/2017, à 18:39)

Ubuntu-fr

Navigation

Liens de recherche

Annonce

#26 Le 16/01/2017, à 17:25

Re : [Resolu] gestion format fasta

#27 Le 16/01/2017, à 18:07

Re : [Resolu] gestion format fasta

#28 Le 16/01/2017, à 18:20

Re : [Resolu] gestion format fasta

#29 Le 16/01/2017, à 18:35

Re : [Resolu] gestion format fasta

#30 Le 16/01/2017, à 21:15

Re : [Resolu] gestion format fasta

#31 Le 16/01/2017, à 23:55

Re : [Resolu] gestion format fasta

#32 Le 17/01/2017, à 00:04

Re : [Resolu] gestion format fasta

#33 Le 17/01/2017, à 00:13

Re : [Resolu] gestion format fasta

#34 Le 17/01/2017, à 00:27

Re : [Resolu] gestion format fasta

#35 Le 17/01/2017, à 01:12

Re : [Resolu] gestion format fasta

#36 Le 17/01/2017, à 08:41

Re : [Resolu] gestion format fasta

#37 Le 17/01/2017, à 11:20

Re : [Resolu] gestion format fasta

#38 Le 17/01/2017, à 11:53

Re : [Resolu] gestion format fasta

#39 Le 17/01/2017, à 12:46

Re : [Resolu] gestion format fasta

#40 Le 17/01/2017, à 12:56

Re : [Resolu] gestion format fasta

#41 Le 17/01/2017, à 16:20

Re : [Resolu] gestion format fasta

#42 Le 17/01/2017, à 17:03

Re : [Resolu] gestion format fasta

#43 Le 17/01/2017, à 17:05

Re : [Resolu] gestion format fasta

#44 Le 17/01/2017, à 17:14

Re : [Resolu] gestion format fasta

#45 Le 17/01/2017, à 17:48

Re : [Resolu] gestion format fasta

#46 Le 17/01/2017, à 17:55

Re : [Resolu] gestion format fasta

#47 Le 17/01/2017, à 18:10

Re : [Resolu] gestion format fasta

#48 Le 17/01/2017, à 18:27

Re : [Resolu] gestion format fasta

#49 Le 17/01/2017, à 18:32

Re : [Resolu] gestion format fasta

#50 Le 17/01/2017, à 18:37

Re : [Resolu] gestion format fasta

Pied de page des forums