Несмотря на то, что я дал свой предыдущий ответ, я рассмотрел подпроблему извлечения конкретной подпоследовательности из файла fasta. Решение состоит из двух частей:
sh
, который выполняет разбор командной строки и вызывает... awk
, который выполняет разбор файла fasta. Я решил разместить это здесь, потому что это показывает
awk
скрипт, в отличие от просто awk
-"однострочников". Предположения:
>
в строке заголовка, за которым следует символ пробела. Скрипт называется extract.sh
.
Запуск, получение последовательности для идентификатора последовательности gi1234
между позициями 36 и 183 (включительно):
$ sh extract.sh -i gi1234 -a 36 -b 183 <sequence.fa
ACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACA
GGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTG
TACAACACTGGGGGACAGTGACG
Вывод неформатирован. В этом случае я взял данные из вопроса и вставил символы перевода строки после каждого 80-го символа перед запуском скрипта.
Скрипт shell (extract.sh
):
#!/bin/sh
usage() {
cat <<USAGE_END
Usage:
sh $0 -i seq_id -a start -b end <sequence.fa
USAGE_END
}
while getopts "a:b:i:" opt; do
case $opt in
a) start_pos="$OPTARG" ;;
b) end_pos="$OPTARG" ;;
i) seq_id="$OPTARG" ;;
*) usage; exit 1 ;;
esac
done
if [ "x$seq_id" = "x" ]; then
echo "Missing sequence ID (-i seq_id)"
usage
exit 1
fi
if [ "x$start_pos" = "x" -o "x$end_pos" = "x" ]; then
echo "Missing either start or end coordinates (-a, -b)"
usage
exit 1
fi
if [ ! -f "extract.awk" ]; then
echo "Can not find the extract.awk AWK script!"
exit 1
fi
awk -v id="$seq_id" -v a="$start_pos" -v b="$end_pos" -f extract.awk
Скрипт awk
(extract.awk
):
# state == 0: not yet at wanted sequence entry in file
# state == 1: found sequence entry, not yet at start position
# state == 2: found start position, not yet at end position
# off: offset in sequence read so far
BEGIN { state = 0 }
state == 0 && $0 ~ "^>" id " " {
# We found the sequence entry we're looking for.
state = 1;
off = 0;
next;
}
state > 0 && /^>/ {
# New sequence header found while parsing sequence, terminate.
exit;
}
state == 1 {
len = length($0);
if (len + off < a) {
# Start position is not on this line.
off += len;
next;
}
state = 2;
if (len + off >= b) {
# Whole sequence is found on this line.
print(substr($0, a - off, b - a + 1))
exit;
}
# Output partial start of sequence.
print(substr($0, a - off , len - (a - off) + 1))
off += len;
next;
}
state == 2 {
len = length($0);
if (off > b) {
# (we should not arrive here)
exit;
}
if (len + off < b) {
# End position is not on this line.
print($0);
off += len;
next;
}
# We found the end position, output line until end position
# and terminate.
print(substr($0, 1, b - off));
exit;
}