Calculer Girl
  1. Calculer Girl
  3. Utiliser awk pour extraire deux chaînes distinctes
Intereting Posts
Quel logiciel ou service puis-je utiliser pour faire des appels téléphoniques par programmation? Ce qui constitue la sécurité asynchrone Comment puis-je transmettre des données d’un thread au processus parent? Comment puis-je désactiver les nouveaux messages thread / thread sortis dans gdb? Sécuriser les abonnements Webcal avec https Utilisation du fichier d’échange dans la programmation c ou java Impossible de démarrer le service Jenkins sur Ubuntu Programmation Unix Bash Shell si le répertoire existe Obtenir l’adresse IP et le port du client? Détails du répertoire temporaire Windows (Java) Comment envoyer une archive compressée contenant des exécutables afin que le filtre de pièces jointes de Google ne le rejette pas Est-il possible de rechercher un indicateur dans une page de manuel? Comment vérifier si sed a changé un fichier C Processus fork () Le serveur XAMPP Apache ne fonctionne pas sur le réseau interne

Utiliser awk pour extraire deux chaînes distinctes

MacOS, Unix

J’ai donc un fichier au format stockholm suivant: 

# STOCKHOLM 1.0 #=GS WP_002855993.1/5-168 DE [subseq from] MULTISPECIES: AAC(3) family N-acetyltransferase [Campylobacter] #=GS WP_002856586.1/5-166 DE [subseq from] MULTISPECIES: aminoglycoside N(3)-acetyltransferase [Campylobacter] WP_002855993.1/5-168 ------LEHNGKKYSDKDLIDAFYQLGIKRGDILCVHTELmkfgKALLT.K...NDFLKTLLECFFKVLGKEGTLLMP-TF---TYSF------CKNE------VYDKVHSKG--KVGVLNEFFRTSGgGVRRTSDPIFSFAVKGAKADIFLKEN--SSCFGKDSVYEILTREGGKFMLLGLNYG-HALTHYAEE----- #=GR WP_002855993.1/5-168 PP ......6788899999***********************9333344455.6...8999********************.33...3544......4555......799999975..68********98626999****************999865..689*********************9875.456799996..... WP_002856586.1/5-166 ------LEFENKKYSTYDFIETFYKLGLQKGDTLCVHTEL....FNFGFpLlsrNEFLQTILDCFFEVIGKEGTLIMP-TF---TYSF------CKNE------VYDKINSKT--KMGALNEYFRKQT.GVKRTNDPIFSFAIKGAKEELFLKDT--TSCFGENCVYEVLTKENGKYMTFGGQG--HTLTHYAEE----- #=GR WP_002856586.1/5-166 PP ......5566677788889999******************....**9953422246679*******************.33...3544......4455......799998876..589**********.******************99999886..689******************999765..5666***96..... #=GC PP_cons ......6677788899999999*****************9....77675.5...68889*******************.33...3544......4455......799999976..689*******998.8999**************99999876..689******************9998765.466699996..... #=GC RF xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx....xxxxx.x...xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx.xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx WP_002855993.1/5-168 ----------------------------------------------------------------------------------------------------- #=GR WP_002855993.1/5-168 PP ..................................................................................................... WP_002856586.1/5-166 ----------------------------------------------------------------------------------------------------- #=GR WP_002856586.1/5-166 PP ..................................................................................................... #=GC PP_cons ..................................................................................................... #=GC RF xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx //

Et j’ai créé un script pour extraire les ID que je veux, dans ce cas, WP_002855993.1 et WP_002856586.1, et rechercher dans un autre fichier pour extraire les séquences d’ADN avec les ID appropriés. Le script est le suivant: 

 #!/bin/bash for fileName in *.sto; do protID=$(grep -o "WP_.\{0,11\}" $fileName | sort | uniq) echo $protID file=$(echo $fileName | cut -d '_' -f 1,2,3) file=$(echo $file'_protein.faa') echo $file if [ -n "$protID" ]; then gawk "/^>/{N=0}/^.*$protID/{N=1} {if(N)print}" $file >> sequence_protein.file fi done

Et voici un exemple du type de fichier que je cherche: 

 >WP_002855993.1 MULTISPECIES: AAC(3) family N-acetyltransferase [Campylobacter] MKYFLEHNGKKYSDKDLIDAFYQLGIKRGDILCVHTELMKFGKALLTKNDFLKTLLECFFKVLGKEGTLLMPTFT >WP_002856586.1 MULTISPECIES: aminoglycoside N(3)-acetyltransferase [Campylobacter] MKYLLEFENKKYSTYDFIETFYKLGLQKGDTLCVHTELFNFGFPLLSRNEFLQTILDCFFEVIGKEGTLIMPTFT YSFCKNEVYDKINSKTKMGALNEYFRKQTGVKRTNDPIFSFAIKGAKEELFLKDTTSCFGENCVYEVLTKENGKY >WP_002856595.1 MULTISPECIES: acetyl-CoA carboxylase biotin carboxylase subunit [Campylobacter] MNQIHKILIANRAEIAVRVIRACRDLHIKSVAVFTEPDRECLHVKIADEAYRIGTDAIRGYLDVARIVEIAKACG

Ce script fonctionne si j’ai un seul identifiant, mais dans certains cas, j’obtiens deux identifiants, et je reçois une erreur, car je pense qu’il recherche un identifiant comme “WP_002855993.1 WP_002856586.1”. Est-il possible de modifier ce script afin qu’il recherche deux occurrences distinctes? Je suppose que c’est quelque chose avec la commande gawk, mais je ne sais pas exactement quoi. Merci d’avance!

une mise à jour du script original: 

 #!/usr/bin/env bash for file_sto in *.sto; do file_faa=$(echo $file_sto | cut -d '_' -f 1,2,3) file_faa=${file_faa}"_protein.faa" awk '(NR==FNR) { match($0,/WP_.\{0,11\}/); if (RSTART > 0) a[substr($0,RSTART,RLENGTH)]++ next; } ($1 in a){ print RS $0 }' $file_sto RS=">" $file_faa >> sequence_protein.file done

La partie awk peut probablement même être réduite à: 

 awk '(NR==FNR) { if ($0 ~ /^WP_/) a[$1]++; next } ($1 in a) { print RS $0 }' FS='/' $file_sto FS=" " RS=">" $file_faa

Ce script awk effectue les opérations suivantes:

  1. Définissez le séparateur de champs FS sur / et lisez le fichier $file_sto .
  2. Lors de la lecture de $file_sto le numéro d’enregistrement NR est identique au numéro d’enregistrement du fichier FNR .
  3. (NR==FNR) { if ($0 ~ /^WP_/) a[$1]++; next } (NR==FNR) { if ($0 ~ /^WP_/) a[$1]++; next } : cette ligne ne fonctionne qu’avec un seul $file_sto raison de la condition au premier plan. Il vérifie si la ligne commence par WP_ . Si c’est le cas, il stocke le premier champ $1 (séparé par FS qui est un / ) dans un tableau a ; il passe ensuite à l’enregistrement suivant dans le fichier ( next ).
  4. Si on a fini de lire le fichier $file_sto , on $file_sto le séparateur de champs en un seul espace FS=" " (voir section Expression régulière ) et le séparateur d’enregistrements RS en > et on commence à lire le fichier $file_faa Ce dernier implique que $0 contiendra toutes les lignes between > et le premier champ $1 est le protID .
  5. En lisant $file_faa , le numéro d’enregistrement du fichier FNR est redémarré à partir de 1 alors que NR n’est pas réinitialisé. Par conséquent, la première ligne awk est ignorée.
  6. ($1 in a){ print RS $0 } si le premier champ est dans le tableau a , imprimez-le avec le séparateur d’enregistrement devant lui.

fixer le script original:

Si vous souhaitez conserver votre script d’origine, vous pouvez stocker le protID dans une liste, puis en boucle: 

 #!/bin/bash for fileName in *.sto; do protID_list=( $(grep -o "WP_.\{0,11\}" $fileName | sort | uniq) ) echo ${protID_list[@]} file=$(echo $fileName | cut -d '_' -f 1,2,3) file=$(echo $file'_protein.faa') echo $file for protID in ${protID_list[@]}; do if [ -n "$protID" ]; then gawk "/^>/{N=0}/^.*$protID/{N=1} {if(N)print}" $file >> sequence_protein.file fi done done

Considérant que votre fichier de sortie est test.

Utiliser la commande suivante ne vous donne que des noms de fichiers: 

 >>cat text | awk '{print $1}' | grep -R 'WP*' | cut -d":" -f2

me donne sortie: 

 WP_002855993.1/5-168 WP_002856586.1/5-166 WP_002855993.1/5-168 WP_002856586.1/5-166

Voulez-vous une sortie comme ça?