Сложно писать что-нибудь полезное в отрыве от реальности, поэтому при подходящем случае буду делиться вариантами решения тех или иных задач.
Недавно необходимо было мне предоставить 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
Легко понять, что разделителем полей является группа символов “ # “, при этом сами поля включают:
- Название контига, в котором был найден ген и его порядковый номер, разделенные последним знаком подчёркивания, наличие нестандартных символов зависит от вашего исходного файла для предсказания генов
- Начало ORF, целое число.
- Конец ORF, целое число.
- Направление считывания: “1” — “+ цепь”, “-1” — “- цепь”.
- Комментарии, включающие полезную и не очень информацию.
Надо заметить, что нумерация нуклеотидов идёт с 1, что сооветствует формату GFF и не соответствует формату BED (нумерация с 0).
Теперь рассмотрим структуру GFF3, в который мы будем конвертировать эту информацию. Формат – текстовый файл, каждому объекту соответствует строка, разбитая знаками табуляции на поля:
- Название последовательности (контиг)
- Источник аннотации
- Тип объекта: gene, CDS, mRNA etc.
- Начало объекта
- Конец объекта
- Score, численная оценка достоверности
- Цепь (+/-)
- Сдвиг (1, 2, 3)
- Комментарий
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, однако в данной ситуации у меня эта возможность не работала, возможно доступны только односимвольные разделители.
Надеюсь было познавательно, если вы заметили какие-то ошибки в коде, хотите предложить другую реализацию или узнали что-то новое и полезное – комментируйте запись 🙂