Embora eu tenha dado a minha resposta anterior, examinei o subproblema de extrair uma subseqüência específica de um arquivo fasta. A solução está em duas partes:
- Um script de shell
sh que faz a análise de linha de comando e chama ...
- Um script
awk que faz a análise do arquivo fasta.
Eu decidi postar isso aqui porque mostra
- Como fazer a análise de linha de comando de opções em um script de shell.
- É possível escrever um script
awk , em vez de apenas awk - "one-liners".
Suposições:
- O ID da sequência ocorrerá diretamente após o
> na linha do cabeçalho, seguido por um caractere de espaço.
- Não há espaços em nenhum lugar nos dados da sequência.
O script é chamado extract.sh .
Para executar, obtenha a sequência para a ID de sequência gi1234 entre as posições 36 e 183 (inclusive):
$ sh extract.sh -i gi1234 -a 36 -b 183 <sequence.fa
ACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACA
GGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTG
TACAACACTGGGGGACAGTGACG
A saída não está formatada. Nesse caso, peguei os dados na pergunta e inseri as quebras de linha após cada caractere 80 antes de executar o script.
O script de 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
O script 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;
}