Отфильтровать файлы FASTA по указанной длине последовательности в bash - PullRequest
Новая
       77

Отфильтровать файлы FASTA по указанной длине последовательности в bash

1 голос
/ 28 мая 2020

Там есть файл FASTA assembly.fasta, содержащий имена контигов и соответствующие последовательности: 

>contig_1
CCAATACGGGCGCGCAGGCTTTCTATCGCGCGGCCGGCTTCGTCGAGGACGGGCGGCGCA
AGGATTACTACCGCAGCGGC
>contig_2
ATATAAACCTTATTCATCGTTTTCAGCCTAATTTTCCATTTAACAGGGATGATTTTCGTC
AAAATGCTGAGGCTTTACCAAGATTTTCTACCTTGCACCTTCAGAAAAAAATCATGGCAT
TTATAGACGAAATTCTCGAGAAA
>contig_3
CGTGATCTCGCCATTCGTGCCG

Я хочу получить только контиги длиной более 30 букв и получить новый файл FASTA assembly.filtered.fasta, содержащий только эти длинные последовательности с именами контигов в этом формате: 

>contig_1
CCAATACGGGCGCGCAGGCTTTCTATCGCGCGGCCGGCTTCGTCGAGGACGGGCGGCGCA
AGGATTACTACCGCAGCGGC
>contig_2
ATATAAACCTTATTCATCGTTTTCAGCCTAATTTTCCATTTAACAGGGATGATTTTCGTC
AAAATGCTGAGGCTTTACCAAGATTTTCTACCTTGCACCTTCAGAAAAAAATCATGGCAT
TTATAGACGAAATTCTCGAGAAA

Ответы [ 3 ]

4 голосов
/ 28 мая 2020

Используя gnu-awk, вы можете использовать эту более простую версию: 

awk -v RS='>[^\n]+\n' 'length() >= 30 {printf "%s", prt $0} {prt = RT}' file


>contig_1
CCAATACGGGCGCGCAGGCTTTCTATCGCGCGGCCGGCTTCGTCGAGGACGGGCGGCGCA
AGGATTACTACCGCAGCGGC
>contig_2
ATATAAACCTTATTCATCGTTTTCAGCCTAATTTTCCATTTAACAGGGATGATTTTCGTC
AAAATGCTGAGGCTTTACCAAGATTTTCTACCTTGCACCTTCAGAAAAAAATCATGGCAT
TTATAGACGAAATTCTCGAGAAA
1 голос
/ 29 мая 2020

Очень быстрый способ достичь того, что вам нужно: 

awk -v n=30 '/^>/{ if(l>n) print b; b=$0;l=0;next }
            {l+=length;b=b ORS $0}END{if(l>n) print b }' file

Возможно, вас также заинтересует BioAwk , это адаптированная версия awk, которая настроена для обработки Файлы FASTA 

bioawk -c fastx -v '(length($seq)>30){print ">" $name ORS $seq}' file.fasta

Примечание: BioAwk основан на awk Брайана Кернигана *, который задокументирован в " Язык программирования AWK »Аль Ахо, Брайана Кернигана и Питера Вайнбергера (Аддисон-Уэсли, 1988, ISBN 0-201-07981-X) . Я не уверен, совместима ли эта версия с POSIX .

1 голос
/ 28 мая 2020

Не могли бы вы попробовать следовать, протестировать и написать с показанными образцами. 

awk '
/^>/{
  if(sign_val && strLen>=30){
    print sign_val ORS line
  }
  strLen=line=""
  sign_val=$0
  next
}
{
  strLen+=length($0)
  line=(line?line ORS:"")$0
}
END{
  if(sign_val && strLen>=30){
    print sign_val ORS line
  }
}
'  Input_file

Пояснение: Добавление подробного объяснения выше. 

awk '                                ##Starting awk program from here.
/^>/{                                ##Checking condition if line starts from > then do following.
  if(sign_val && strLen>=30){        ##Checking if sign_val is SET and steLen is SET then do following.
    print sign_val ORS line          ##Printing sign_val ORS and line here.
  }
  strLen=line=""                     ##Nullify variables steLen and line here.
  sign_val=$0                        ##Setting sign_val to current line here.
  next                               ##next will skip all further statements from here.
}
{
  strLen+=length($0)                 ##Checking length of line and keep adding it here.
  line=(line?line ORS:"")$0          ##Creating line variable and keep appending it to it with new line.
}
END{                                 ##Starting END block of this program from here.
  if(sign_val && strLen>=30){        ##Checking if sign_val is SET and steLen is SET then do following.
    print sign_val ORS line          ##Printing sign_val ORS and line here.
  }
}
' Input_file                         ##mentioning Input_file name here.
...