☰

Como cortar sequência fasta em números específicos e gerar ORFs

tag-icon • tag-icon text-processing bioinformatics
Como cortar sequência fasta em números específicos e gerar ORFs

Tenho um arquivo conforme mostrado abaixo:

 CDS             join(36..56,37..67)
 CDS             36..183
 CDS             457..565
 CDS             join(505..519,521..596)
 CDS             join(577..591,725..770)
 CDS             join(516..591,725..899)
 CDS             508..556
 CDS             571..841
 CDS             complement(619..788)
 CDS             843..863

Quero imprimir o número específico do intervalo de nucleotídeos como no arquivo (a sequência é lida de outro arquivo "sequence.fasta"). Por exemplo, para o arquivo sequencia.fasta como:

>gi1234 HIVgenome|NC_909999.1
AACTGCGTGTGTGTCCACACAACACTGGGGGACACACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACAGGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTGTACAACACTGGGGGACAGTGACGACGACAACACTGGGGGACACGAGCGTTGTGAGCAGGTGACAACACTGGGGGACAGTGTTTTTACAACACTGGGGGACATTTTTGAGCAGCGACGCAGCGTTGTGGGGTGTGTCGGAAGGTGTGTCGTGTGTCGTGTGTC

A saídap deve ser

36  -  56   ACAACAACAACACTGGGGGAC 

37  -  67   CAACAACAACACTGGGGGACAACACTGGGAC

& breve...

até

843 - 863   GTGT....

Qual seria a maneira mais fácil de fazer isso por meio de scripts de shell?

Responder1

Esta questão requer um esforço de programação maior do que pode ser oferecido por este fórum (eu faço esse tipo de programação para viver).

OFormato de arquivo DDBJ/ENA/GenBank(o primeiro arquivo da questão) é complexo e permite que os CDSs (as partes codificantes de uma sequência genômica) sejam não apenas simples ou unidos, mas complementados e suas combinações. Além disso, as coordenadas posicionais podem termodificadoresque, para uma solução genérica, precisaria ser tratado.

Seria melhor perguntar a um bioinformático local (ou programador) ou em um fórum de bioinformática, como o StackExchangeSite de bioinformática. Eles indicarão ferramentas existentes para fazer esse tipo de coisa ou, conhecendo os bioinformáticos, fornecerão algum script BioPerl/BioPython peculiar que provavelmente funcionará com mais frequência ;-)

Umpossívelcaminho seria usar oExtrator de recursos GenBank, mas usá-lo on-line provavelmente não será a melhor opção para nada além de pequenos conjuntos de dados.

Responder2

Embora eu tenha dado minha resposta anterior, examinei o subproblema de extrair uma subsequência específica de um arquivo fasta. A solução está em duas partes:

  1. Um shscript de shell que analisa a linha de comando e chama ...
  2. Um awkscript que faz a análise do arquivo fasta.

Resolvi postar isso aqui porque mostra

  1. Como fazer a análise de opções de linha de comando em um script de shell.
  2. Que é possível escrever um awkroteiro, em vez de apenas awk“one-liners”.

Premissas:

  • O ID de sequência ocorrerá diretamente após 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, obtendo a sequência para ID de sequência gi1234entre 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 da pergunta e inseri quebras de linha a cada 80 caracteres 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 awkroteiro ( 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;
}

Responder3

Eu estava resolvendo um problema semelhante, mas costumava ter alguns problemas ao tentar usar os scripts fornecidos em uma das respostas (ambos WSL, Mac) ...

Continuei pesquisando e encontrei uma solução bem fácil: é a getfastafunção da bedtoolssuíte!

https://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html

informação relacionada