Использование awk для извлечения двух отдельных строк - PullRequest
Новая
       35

Использование awk для извлечения двух отдельных строк

0 голосов
/ 26 апреля 2018

MacOS, Unix

Итак, у меня есть файл в следующем стокгольмском формате: 

# 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
//

И я создал скрипт для извлечения необходимых мне идентификаторов, в данном случае,WP_002855993.1 и WP_002856586.1 и выполните поиск в другом файле, чтобы извлечь последовательности ДНК с соответствующими идентификаторами.Сценарий выглядит следующим образом: 

#!/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

И вот пример типа файла, который я просматриваю: 

>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

Этот сценарий работает, если у меня есть один идентификатор, но вВ некоторых случаях я получаю два идентификатора и получаю сообщение об ошибке, потому что я думаю, что он ищет такой идентификатор, как "WP_002855993.1 WP_002856586.1".Есть ли способ изменить этот скрипт так, чтобы он смотрел на два отдельных случая?Я думаю, что-то с командой gawk, но я не уверен, что именно.Заранее спасибо!

Ответы [ 2 ]

0 голосов
/ 26 апреля 2018

обновление оригинального сценария: 

#!/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

Часть awk может быть даже уменьшена до: 

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

This awkСценарий выполняет следующее:

  1. Установите разделитель полей FS на / и прочитайте файл $file_sto.
  2. При чтении $file_sto номер записи NR равентакой же, как номер записи файла FNR.
  3. (NR==FNR) { if ($0 ~ /^WP_/) a[$1]++; next }: эта строка работает только на одну $file_sto из-за условий на передней панели.Он проверяет, начинается ли строка с WP_.Если это так, оно сохраняет первое поле $1 (разделенное FS, которое является /) в массиве a;затем он переходит к следующей записи в файле (next).
  4. Если мы закончили чтение файла $file_sto, мы устанавливаем разделитель полей обратно на один пробел FS=" " (см. раздел ) Регулярное выражение ) и разделитель записей RS до > и начало чтения файла $file_faa Последнее подразумевает, что $0 будет содержать все строки между > и первым полем $1 - это protID.
  5. Чтение $file_faa, номер записи файла FNR перезапускается с 1, а NR не сбрасывается.Следовательно, первая строка awk пропускается.
  6. ($1 in a){ print RS $0 }, если первое поле находится в массиве a, распечатайте запись с разделителем записей перед ней.

исправление исходного сценария:

Если вы хотите сохранить исходный сценарий, вы можете сохранить protID в списке и затем зациклить список: 

#!/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
0 голосов
/ 26 апреля 2018

Учитывая, что ваш выходной файл тестовый.

Использование следующей команды дает вам только имена файлов: 

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

дает мне вывод: 

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

Хотите ли вывывод такой?

...