☰

So schneiden Sie eine Fasta-Sequenz an bestimmten Stellen ab und generieren ORFs

tag-icon • tag-icon text-processing bioinformatics
So schneiden Sie eine Fasta-Sequenz an bestimmten Stellen ab und generieren ORFs

Ich habe eine Datei wie unten gezeigt:

 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

Ich möchte die spezifische Nummer des Nukleotidbereichs wie in der Datei drucken (die Sequenz wird aus einer anderen Datei „sequence.fasta“ gelesen). Zum Beispiel für die Datei sequence.fasta wie folgt:

>gi1234 HIVgenome|NC_909999.1
AACTGCGTGTGTGTCCACACAACACTGGGGGACACACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACAGGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTGTACAACACTGGGGGACAGTGACGACGACAACACTGGGGGACACGAGCGTTGTGAGCAGGTGACAACACTGGGGGACAGTGTTTTTACAACACTGGGGGACATTTTTGAGCAGCGACGCAGCGTTGTGGGGTGTGTCGGAAGGTGTGTCGTGTGTCGTGTGTC

Outputp sollte sein

36  -  56   ACAACAACAACACTGGGGGAC 

37  -  67   CAACAACAACACTGGGGGACAACACTGGGAC

& bald...

bis

843 - 863   GTGT....

Was wäre der einfachste Weg, dies über Shell-Skripte zu tun?

Antwort1

Diese Frage erfordert einen größeren Programmieraufwand, als dieser in diesem Forum möglich wäre (ich verdiene meinen Lebensunterhalt mit dieser Art der Programmierung).

DerDDBJ/ENA/GenBank-Dateiformat(die erste Datei in der Frage) ist komplex und ermöglicht CDSs (die codierenden Teile einer Genomsequenz) nicht nur einfach oder verbunden, sondern auch ergänzt und Kombinationen davon. Darüber hinaus können die PositionskoordinatenModifikatorendie für eine generische Lösung behandelt werden müssten.

Am besten fragen Sie einen lokalen Bioinformatiker (oder Programmierer) oder in einem Bioinformatik-Forum wie StackExchange.Bioinformatik-Site. Sie werden Sie auf vorhandene Tools für derartige Aufgaben hinweisen oder Ihnen, wenn Sie sich mit Bioinformatik auskennen, ein ausgefallenes BioPerl/BioPython-Skript geben, das wahrscheinlich in den meisten Fällen funktioniert ;-)

EinsmöglichDer Weg wäre, dieGenBank-Feature-Extraktor, aber die Online-Verwendung ist wahrscheinlich nur für kleine Datensätze die beste Option.

Antwort2

Obwohl ich meine vorherige Antwort gegeben habe, habe ich mir das Unterproblem des Extrahierens einer bestimmten Untersequenz aus einer Fasta-Datei angesehen. Die Lösung besteht aus zwei Teilen:

  1. Ein shShell-Skript, das die Befehlszeilenanalyse durchführt und aufruft …
  2. Ein awkSkript, das die Analyse der Fasta-Datei durchführt.

Ich habe mich entschlossen, dies hier zu veröffentlichen, weil es zeigt

  1. So führen Sie eine Befehlszeilenanalyse von Optionen in einem Shell-Skript durch.
  2. Dass es möglich ist, ein Skript zu schreiben awk, und nicht nur awkein paar „Einzeiler“.

Annahmen:

  • Direkt danach steht >in der Kopfzeile die Sequenz-ID, gefolgt von einem Leerzeichen.
  • In den Sequenzdaten sind nirgends Leerzeichen vorhanden.

Das Skript heißt extract.sh.

So führen Sie es aus: Abrufen der Sequenz für die Sequenz-ID gi1234zwischen den Positionen 36 und 183 (einschließlich):

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

Die Ausgabe ist unformatiert. In diesem Fall habe ich die Daten aus der Frage genommen und vor dem Ausführen des Skripts nach jedem 80. Zeichen einen Zeilenumbruch eingefügt.

Das Shell-Skript ( 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

Das awkSkript ( 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;
}

Antwort3

Ich habe ein ähnliches Problem gelöst, aber beim Versuch, die in einer der Antworten bereitgestellten Skripts zu verwenden (sowohl WSL als auch Mac), traten einige Schwierigkeiten auf ...

Also suchte ich weiter und fand eine ganz einfache Lösung: Es ist die getfastaFunktion aus der bedtoolsSuite!

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

verwandte Informationen