特定の数のファスタ配列を切断してORFを生成する方法

特定の数のファスタ配列を切断してORFを生成する方法

次のファイルがあります。

 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

ファイル内の特定の数のヌクレオチド範囲を印刷したいです(他のファイル「sequence.fasta」からシーケンスを読みます)。たとえば、Sequence.fasta ファイルは次のようになります。

>gi1234 HIVgenome|NC_909999.1
AACTGCGTGTGTGTCCACACAACACTGGGGGACACACAACAACAACACTGGGGGACACACTGGGACAACACTGGGGGACAGGACACTGTACAACACTGGGTGTGTCGGGACAGTACACATGTTGGGGGGGTGTGTCGGACAACACTGGGGGACATGTGTGTACAACACTGGGGGACAGTGACGACGACAACACTGGGGGACACGAGCGTTGTGAGCAGGTGACAACACTGGGGGACAGTGTTTTTACAACACTGGGGGACATTTTTGAGCAGCGACGCAGCGTTGTGGGGTGTGTCGGAAGGTGTGTCGTGTGTCGTGTGTC

出力pは次のようになります。

36  -  56   ACAACAACAACACTGGGGGAC 

37  -  67   CAACAACAACACTGGGGGACAACACTGGGAC

&まもなく...

~まで

843 - 863   GTGT....

シェルスクリプトでこれを行う最も簡単な方法は何ですか?

答え1

この問題を解決するには、このフォーラムで提供するよりも多くのプログラミング努力が必要です(私はこの種のプログラミングを専門としています)。

これDDBJ/ENA/GenBank ファイル形式(質問の最初のファイル)は複雑で、CDS(ゲノム配列のコード部分)が単純または連結されているだけでなく、相補的であり、それらの組み合わせでもあります。また、位置座標には次のものがあります。修飾子一般的な解決策では、この問題を処理する必要があります。

地域の生物情報学者(またはプログラマー)に問い合わせるか、StackExchangeなどの生物情報学フォーラムに連絡する方が良いでしょう。生物情報学ウェブサイト。この種のタスクを実行するための既存のツールを教えたり、生物情報学者を知っていれば、より効率的である可能性のある奇妙なBioPerl / BioPythonスクリプトを提供します。

一つ可能パスを使用しています。GenBank特徴抽出器しかし、オンラインで使用することは、小規模なデータセット以外の場合は最善の選択ではないかもしれません。

答え2

以前の答えにもかかわらず、fastaファイルから特定のサブシーケンスを抽出するサブ問題を調査しました。解決策は2つの部分に分けられます。

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

見せてくれるのでここに投稿することにしました。

  1. シェルスクリプトでオプションのコマンドライン解析を実行する方法。
  2. 単純な「一行」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

スクリプトawkextract.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;
}

答え3

同様の問題を解決していますが、過去に答え(WSL、Mac)の1つに提供されているスクリプトを使用しようとしたときにいくつかの問題がありました...

だから私は探し続け、とても簡単な解決策を見つけました。それがキットgetfastaに含まれていましたbedtools

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

関連情報