☰

Cómo cortar una secuencia fasta en números específicos y generar ORF

tag-icon • tag-icon text-processing bioinformatics
Cómo cortar una secuencia fasta en números específicos y generar ORF

Tengo un archivo como se muestra a continuación:

 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

Quiero imprimir el número específico del rango de nucleótidos como en el archivo (la secuencia se lee desde otro archivo "sequence.fasta"). Por ejemplo, para el archivo secuencia.fasta como:

>gi1234 HIVgenome|NC_909999.1
AACTGCGTGTGTGTCCACACAACACTGGGGGACACACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACAGGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTGTACAACACTGGGGGACAGTGACGACGACAACACTGGGGGACACGAGCGTTGTGAGCAGGTGACAACACTGGGGGACAGTGTTTTTACAACACTGGGGGACATTTTTGAGCAGCGACGCAGCGTTGTGGGGTGTGTCGGAAGGTGTGTCGTGTGTCGTGTGTC

La salida debe ser

36  -  56   ACAACAACAACACTGGGGGAC 

37  -  67   CAACAACAACACTGGGGGACAACACTGGGAC

& pronto...

hasta

843 - 863   GTGT....

¿Cuál sería la forma más sencilla de hacerlo mediante scripts de shell?

Respuesta1

Esta pregunta requiere un esfuerzo de programación mayor que el que puede ofrecer este foro (me dedico a este tipo de programación).

ElFormato de archivo DDBJ/ENA/GenBank(el primer archivo de la pregunta) es complejo y permite que los CDS (las partes codificantes de una secuencia genómica) no sean simplemente simples o unidos, sino complementados y combinaciones de los mismos. Además, las coordenadas posicionales pueden tenermodificadoresque, para una solución genérica, tendrían que ser manejados.

Sería mejor que le preguntaras a un bioinformático (o programador) local o en un foro de bioinformática, como StackExchange.Sitio de bioinformática. Le indicarán las herramientas existentes para hacer este tipo de cosas o, conociendo a los bioinformáticos, le darán algún script peculiar de BioPerl/BioPython que probablemente funcionará la mayoría de las veces ;-)

UnoposibleLa ruta sería utilizar elExtractor de funciones de GenBank, pero usarlo en línea probablemente no sea la mejor opción para nada que no sean conjuntos de datos pequeños.

Respuesta2

Aunque di mi respuesta anterior, analicé el subproblema de extraer una subsecuencia particular de un archivo fasta. La solución consta de dos partes:

  1. Un shscript de shell que analiza la línea de comando y llama...
  2. Un awkscript que analiza el archivo fasta.

Decidí publicar esto aquí porque muestra

  1. Cómo realizar un análisis de opciones en la línea de comandos en un script de shell.
  2. Que es posible escribir un awkguión, en lugar de simplemente awk"frases ingeniosas".

Supuestos:

  • El ID de la secuencia aparecerá directamente después de >en la línea del encabezado, seguido de un carácter de espacio.
  • No hay espacios en ninguna parte de los datos de la secuencia.

El guión se llama extract.sh.

Para ejecutar, obteniendo la secuencia para el ID de secuencia gi1234entre las posiciones 36 y 183 (inclusive):

$ sh extract.sh -i gi1234 -a 36 -b 183 <sequence.fa
ACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACA
GGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTG
TACAACACTGGGGGACAGTGACG

La salida no está formateada. En este caso, tomé los datos de la pregunta e inserté saltos de línea después de cada 80 caracteres antes de ejecutar el script.

El 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

La awksecuencia de comandos ( 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;
}

Respuesta3

Estaba resolviendo un problema similar, pero solía tener algunos problemas al intentar usar los scripts proporcionados en una de las respuestas (tanto WSL como Mac)...

Así que seguí buscando y encontré una solución bastante sencilla: ¡es la getfastafunción de la bedtoolssuite!

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

información relacionada