Разбор задач №1. GFF из аннотаций Prodigal

Сложно писать что-нибудь полезное в отрыве от реальности, поэтому при подходящем случае буду делиться вариантами решения тех или иных задач.

Недавно необходимо было мне предоставить GFF файл с указанием генов, найденных в метагеноме. Всё бы хорошо, fasta-файлы с аминокислотными и нуклеотидными последовательностями сохранились, а вот ни GFF, ни GBK не сохранилось. К счастью, Prodigal в процессе поиска ORF сохраняет много данных в заголовках мультифасты, выглядит это примерно так:

>Contig_name_1_1 # 3 # 215 # -1 # ID=1_1;partial=10;start_type=ATG;rbs_motif=AGGAG;rbs_spacer=5-10bp;gc_cont=0.629

Легко понять, что разделителем полей является группа символов “ # “, при этом сами поля включают:

  1. Название контига, в котором был найден ген и его порядковый номер, разделенные последним знаком подчёркивания, наличие нестандартных символов зависит от вашего исходного файла для предсказания генов
  2. Начало ORF, целое число.
  3. Конец ORF, целое число.
  4. Направление считывания: “1” — “+ цепь”, “-1” — “- цепь”.
  5. Комментарии, включающие полезную и не очень информацию.

Надо заметить, что нумерация нуклеотидов идёт с 1, что сооветствует формату GFF и не соответствует формату BED (нумерация с 0).

Теперь рассмотрим структуру GFF3, в который мы будем конвертировать эту информацию. Формат – текстовый файл, каждому объекту соответствует строка, разбитая знаками табуляции на поля:

  1. Название последовательности (контиг)
  2. Источник аннотации
  3. Тип объекта: gene, CDS, mRNA etc.
  4. Начало объекта
  5. Конец объекта
  6. Score, численная оценка достоверности
  7. Цепь (+/-)
  8. Сдвиг (1, 2, 3)
  9. Комментарий

grep ">" proteins.faa | awk -F " # " 'BEGIN{OFS="\t"}{S="-";if($4=="1")S="+";print($1, ".", "gene", $2, $3, ".", S, ".", $5)}' | sed 's/>//' > output.gff

Здесь мы отбираем все заголовки, передаём их в awk, устанавливаем разделителем “ # “ ключом -F, ставим разделителем вывода табуляцию в блоке BEGIN. В теле скрипта, выполняемом для каждой строки, объявляем переменную S, хранящую минус по умолчанию, проверяем значение 4-го поля, меняем переменную на плюс, если оно равно “1”. Далее выводим поля, полученные разбиением исходной строки: они хранятся в переменных от $1 до $N, где N -количество полей, расставляя точки там, где у нас нет информации (в списке выделено курсивом). Вывод отправляем в sed, убирающий первое упоминание “>” в строке (если буквально, то мы заменяем на ничто).

Почти! Только вот незадача, в первом поле у нас остаётся порядковый индекс гена, его нужно убрать, например так:

В нашем примере мы знаем, что заменяемая позиция будет первой соответствующей шаблону “_[0-9]+\t” — “подчёркивание, одна или несколько цифр, знак табуляции”

sed 's/_[0-9]\+\t/\t/' output.gff > final.gff

Однако, если бы нам нужно было уточнить поле, в котором это случаетя, то на помощь придёт команда paste, которая объединяет все входящие потоки построчно в один файл, по умолчанию разделитель — табуляция.

paste <(cut -f 1 output.gff | sed 's/_[0-9]\+$//') <(cut -f 2- output.gff) > final.gff

Итого мы получаем:

grep ">" proteins.faa | awk -F " # " 'BEGIN{OFS="\t"}{S="-";if($4=="1")S="+";print($1, ".", "gene", $2, $3, ".", S, ".", $5)}' | sed 's/>//;s/_[0-9]\+\t/\t/' > final.gff

Всё это реализуемо и виде однострочника на Perl:

perl -lne 'next unless /^>/;($i,$s,$e,$d,$c)=split/ # /;@i=split/_/,$i;pop @i;$i=join("_",@i);$i=~s/>//;if($d=="-1"){$d="-"}else{$d="+"};print(join("\t",$i,".","gene",$s,$e,".",$d,".",$c))' proteins.faa > final.gff

Здесь мы отбираем строки, начинающиеся с “>”, разделяем в массив, разделяем в массив первое поле по знаку подчёркивания, отбрасываем последний элемент, соединяем в строку, определяем цепь, выводим строку, разделенную табуляцией. В Perl есть возможность сразу работать с массивом @f, содержащим элементы, полученные разбиением входной строки по разделителю, указанному в ключе -F, однако в данной ситуации у меня эта возможность не работала, возможно доступны только односимвольные разделители.

Надеюсь было познавательно, если вы заметили какие-то ошибки в коде, хотите предложить другую реализацию или узнали что-то новое и полезное – комментируйте запись 🙂

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *