Доступ к файлу Fasta с помощью Bio::DB::Fasta

Я использовал модуль use Bio::DB::Fasta для доступа к файлам fasta (документация здесь:https://metacpan.org/pod/Bio::DB::Fasta#OBJECT-METHODS). Я считаю, что это гораздо быстрее, чем использование Samtools для извлечения позиций из файла fasta. Однако мне было интересно, знает ли кто-нибудь, что произойдет, если запрос включает позицию, превышающую максимальную длину фаста.

Сегодня в запросе я попытался получить доступ к позиции в фасте, которая превышает максимальную позицию в фасте. Однако в этом случае метод не дал ошибки. Мой файл fasta содержит базы 0/1, и возвращаемый результат был «1». Мне было интересно, является ли это ошибкой или на самом деле она дает правильный вывод, но для неправильной позиции. Я пытался просмотреть документацию, но не смог найти никакой информации о кодах ошибок.

Мой код выглядит следующим образом:

use strict;
use warnings;
use Bio::DB::Fasta;

my $maskFile = "1KG_maskfile.fa";

my $db = Bio::DB::Fasta->new($maskFile);

my $chrom = "chr1";
my $start = 300240548;
my $end = 300240548;
my $query = "$chrom:$start-$end"; 
my $seq = $db->seq($query, $start, $end); # also tried $seq = $db->seq($query); 
print $seq, "\n";

Примечание. В файле 1KG_maskfile.fa максимальная позиция — 249224750 (на основе количества символов, исключая заголовок).


fasta bioperl
person user19758    schedule 04.02.2014    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Здесь я вижу две проблемы. Первая заключается в том, что вы неправильно форматируете идентификатор запроса, если только у вас нет начальной/конечной позиции в заголовке Fasta (что было бы странно). Чтобы получить нужную последовательность по регионам, достаточно указать конкретный ID и координаты, т.е.

my $seq = $db->seq('chr1', 25000, 27000);

Другая проблема, о которой вы упомянули, выглядит как ошибка. Я не думаю, что есть какие-то явные проверки, если позиции начала/остановки выходят за пределы фактической длины последовательности. Я только что проверил это, и метод молчал. В этом коде есть много других проверок формата, об этом может быть полезно сообщить как об ошибке.

person SES    schedule 05.02.2014
comment
На самом деле, запрос должен быть очень гибким, и выполнение моей $seq = $db-›seq($query, $start, $end) должно дать тот же ответ, что и ваш запрос. Я эмпирически протестировал метод, используя данные 1 кг и сравнивая его с samtools. Однако программа выдает неверный ответ, если вы ошибаетесь за пределами фаста. Я написал автору по электронной почте и надеюсь, что исправление должно быть относительно легко реализовать. Тем временем я нашел FaSlice, который выполняет ту же функцию. search.cpan.org/~ajpage/Bio- Pipeline-Comparison-1.123050/lib/ - person user19758; 07.02.2014
comment
@ user19758, что значит дает неправильный ответ? Я не знаю, что делать с этой информацией. Пожалуйста, покажите пример. Кроме того, если вы знаете, что указываете запрос за пределами диапазона последовательности, какой смысл в этом упражнении? Наконец, вы говорите, что запрос должен быть гибким. Я не уверен, что вы имеете в виду мое это. Запрос не является гибким, это идентификатор последовательности, который вы хотите запросить. - person SES; 07.02.2014