特定の番号で fasta 配列をカットして ORF を生成する方法

Question 1

この質問には、このフォーラムで提供されるものよりも大きなプログラミング作業が必要です (私はこの種のプログラミングを生業としています)。

のDDBJ/ENA/GenBank ファイル形式（質問の最初のファイル）は複雑で、CDS（ゲノム配列のコード部分）を単純に結合したり、補完したり、組み合わせたりすることができます。さらに、位置座標には修飾子一般的なソリューションでは、これを処理する必要があります。

地元のバイオインフォマティクス専門家（またはプログラマー）に尋ねるか、StackExchangeなどのバイオインフォマティクスフォーラムで尋ねるのが良いでしょう。バイオインフォマティクスサイト彼らは、このようなことを実行するための既存のツールを紹介したり、バイオインフォマティクスに精通した専門家が、たいていはうまく機能するであろう、風変わりな BioPerl/BioPython スクリプトを紹介したりします ;-)

1つ可能ルートはGenBank 特徴抽出ツールただし、小規模なデータセット以外では、オンラインで使用するのはおそらく最適なオプションではありません。

Answer

この質問には、このフォーラムで提供されるものよりも大きなプログラミング作業が必要です (私はこの種のプログラミングを生業としています)。

のDDBJ/ENA/GenBank ファイル形式（質問の最初のファイル）は複雑で、CDS（ゲノム配列のコード部分）を単純に結合したり、補完したり、組み合わせたりすることができます。さらに、位置座標には修飾子一般的なソリューションでは、これを処理する必要があります。

地元のバイオインフォマティクス専門家（またはプログラマー）に尋ねるか、StackExchangeなどのバイオインフォマティクスフォーラムで尋ねるのが良いでしょう。バイオインフォマティクスサイト彼らは、このようなことを実行するための既存のツールを紹介したり、バイオインフォマティクスに精通した専門家が、たいていはうまく機能するであろう、風変わりな BioPerl/BioPython スクリプトを紹介したりします ;-)

1つ可能ルートはGenBank 特徴抽出ツールただし、小規模なデータセット以外では、オンラインで使用するのはおそらく最適なオプションではありません。

Question 2

前回の回答でも述べましたが、Fasta ファイルから特定のサブシーケンスを抽出するというサブ問題を検討しました。解決策は 2 つの部分に分かれています。

shコマンドライン解析を実行し、呼び出すシェルスクリプト...
awkfasta ファイルの解析を実行するスクリプト。

これをここに投稿することにしたのは、

シェルスクリプトでオプションのコマンドライン解析を行う方法。
単なる「ワンライナー」awkではなく、スクリプトを書くことも可能です。awk

前提:

シーケンス ID は>ヘッダー行のの直後に表示され、その後にスペース文字が続きます。
シーケンスデータ内にはスペースはありません。

スクリプトはと呼ばれますextract.sh。

gi1234実行するには、位置 36 から 183 (含む) までのシーケンス ID のシーケンスを取得します。

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

出力はフォーマットされていません。この場合、質問のデータを取得し、スクリプトを実行する前に 80 文字ごとに改行を挿入しました。

シェルスクリプト（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;
}

Answer