3

1.blastこのような座標情報を含むファイルがあります

1       gnl|BL_ORD_ID|0 100.00  33      0       0       1        3
27620   gnl|BL_ORD_ID|0 95.65   46      2       0       1       46
35296   gnl|BL_ORD_ID|0 90.91   44      4       0       3       46
35973   gnl|BL_ORD_ID|0 100.00  45      0       0       1       45
41219   gnl|BL_ORD_ID|0 100.00  27      0       0       1       27
46914   gnl|BL_ORD_ID|0 100.00  45      0       0       1       45

1.fastaこのようなシーケンス情報を含むファイル

>1
TCGACTAGCTACGACTCGGACTGACGAGCTACGACTACGG
>2
GCATCTGGGCTACGGGATCAGCTAGGCGATGCGAC
...
>100000
TTTGCGAGCGCGAAGCGACGACGAGCAGCAGCGACTCTAGCTACTG

私は今、最初の列から取得し1.blast、それらのシーケンス ID (=最初の列) とシーケンスを抽出し、次にシーケンス自体からファイル間の位置とファイルからの$1位置を除くすべてを抽出するスクリプトを検索しています。つまり、最初の 2 つの一致から出力が$7$81.fasta

>1
ACTAGCTACGACTCGGACTGACGAGCTACGACTACGG
>27620
GTAGATAGAGATAGAGAGAGAGAGGGGGGAGA
...

>1( の最初の 3 つのエントリはこの順序ではないことに注意してください)

ID は連続しています。つまり、必要な情報を次のように抽出できます。

awk '{print 2*$1-1, 2*$1, $7, $8}' 1.blast

これにより、最初の列に正しいシーケンス識別子行、2 番目の列に正しいシーケンス行 (= ID 行の 1 つ後)、および除外する必要がある 2 つの座標を含むマトリックスが得られます。1.fastaしたがって、基本的には、要素が抽出されるすべての必要な情報を含むマトリックス

残念ながら、私はスクリプト作成の経験があまりないため、適切なsedコマンドなどで値を入力する方法がわかりません。次のような特定の行を取得できます。

sed -n 3,4p 1.fasta

そして、例えば経由で削除したい文字列

sed -n 5p 1.fasta | awk '{print substr($0,2,5)}'

しかし、私の問題は、最初のawk呼び出しからの情報を他のコマンドにパイプして、正しい行を抽出し、シーケンス行から指定された座標を削除する方法です。したがって、substr正しいコマンドではありません。特定の文字列からこれら 2 つの位置の間のすべてを削除するコマンドが必要remstr(string,start,stop)ですが、独自のスクリプトで実行できると思います。特に正しい配管は、ここで私にとって問題です。

4

4 に答える 4

1

thunkとmswのいずれかが指摘しているように、この種のタスクにはより適切なツールを利用できますが、ここでは、それを処理する方法について何かを教えてくれるスクリプトがありますawk。

script.awkの内容:

## Process first file from arguments.
FNR == NR {
        ## Save ID and the range of characters to remove from sequence.
        blast[ $1 ] = $(NF-1) " " $NF
        next
}

## Process second file. For each FASTA id...
$1 ~ /^>/ {
        ## Get number.
        id = substr( $1, 2 )

        ## Read next line (the sequence).
        getline sequence

        ## if the ID is one found in the other file, get ranges and
        ## extract those characters from sequence.
        if ( id in blast ) {
                split( blast[id], ranges )
                sequence = substr( sequence, 1, ranges[1] - 1 ) substr( sequence, ranges[2] + 1 )
                ## Print both lines with the shortened sequence.
                printf "%s\n%s\n", $0, sequence
        }

}

あなたの質問とそれをテストするために1.blastaカスタマイズされたと仮定します：1.fasta

>1
TCGACTAGCTACGACTCGGACTGACGAGCTACGACTACGG
>2
GCATCTGGGCTACGGGATCAGCTAGGCGATGCGAC
>27620
TTTGCGAGCGCGAAGCGACGACGAGCAGCAGCGACTCTAGCTACTGTTTGCGA

次のようにスクリプトを実行します。

awk -f script.awk 1.blast 1.fasta

それは以下をもたらします：

>1
ACTAGCTACGACTCGGACTGACGAGCTACGACTACGG
>27620
TTTGCGA

もちろん、私はいくつかのことを仮定しています.fastaシーケンスが1行より長くないことが最も重要です。

于 2013-05-30T13:34:03.553 に答える