提取行和子上的另一個文件的信息

我有一個文件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 
TCGACTAGCTACGACTCGGACTGACGAGCTACGACTACGG 
>2 
GCATCTGGGCTACGGGATCAGCTAGGCGATGCGAC 
... 
>100000 
TTTGCGAGCGCGAAGCGACGACGAGCAGCAGCGACTCTAGCTACTG

我序列信息的文件1.fasta信息我現在正在搜索一個腳本，該腳本從第一列中取出1.blast並提取這些序列ID（=第一列$1）加上序列，然後從序列本身除了t從1.fasta文件，從頭兩場比賽意味着$7和$8之間的軟管位置的輸出將

>1 
ACTAGCTACGACTCGGACTGACGAGCTACGACTACGG 
>27620 
GTAGATAGAGATAGAGAGAGAGAGGGGGGAGA 
...

（請注意，從>1前三項都沒有按照這個順序）

的標識是連續的，這意味着我可以提取這樣的所需的信息：

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

這使我然後包含在第一列的矩陣的右序列標識符RO w，在第二列中右邊的序列行（=在ID行後面），然後是應該排除的兩個座標。所以基本上是一個矩陣，其中包含1.fasta的所有必需信息將被提取

不幸的是，我沒有太多的腳本經驗，因此我現在有點迷路了，我如何輸入數據在適當的sed命令？我能得到特定的行這樣的：

sed -n 3,4p 1.fasta

而且我想刪除例如串通過

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

但我的問題是現在，我怎麼能管從第一awk呼叫到其他命令，讓他們提取右行，從序列行，那麼定座標刪除信息。因此，substr不是正確的命令，我需要一個命令remstr(string,start,stop)，從給定的字符串中移除這兩個位置之間的所有內容，但我認爲我可以在自己的腳本中完成。特別是正確的管道對我來說是個問題。

來源

2013-05-30 Daniel Fischer

由於無論是咚和垃圾人士指出，更適合的工具可用於這樣的任務，但在這裏你有一個腳本，可以教你一些關於如何與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序列不再是臨時的ñ一行。

來源

2013-05-30 13:34:03 Birei

你的最後一個假設是危險的，因爲[fasta-file格式]（https://en.wikipedia.org/wiki/FASTA_format）（它的用途往往遠不是標準化的）通常有一個固定的換行符，字符數量。 Bioperl或其他類似的原因之一是各自的方法會考慮（大部分）文件格式的變化。 – thunk

非常感謝！我仍然嘗試通過代碼並理解那裏發生的事情並對其進行調整（例如，腳本的輸出不應包含序列'2'，因爲它在'blast'文件中未提及）。但我希望自己能夠從這個角度來解決問題。 –

@thunk，有危險的假設點，但有一點我沒有在這裏提到 - 所有這些練習只是爲了運行某些爆炸搜索，我不使用任何其他工具，我只是在這裏叫fasta，但我做不要在任何其他程序中提供文件，我只需要序列進行第二次高速搜索，以前未找到的部分。 –

如果您使用生物信息學並使用DNA序列（或者更復雜的事物，如序列註釋），我會建議看看Bioperl。這顯然需要Perl的知識，但具有相當多的功能。

在你的情況，你會想要從你的fasta文件使用Bio::SeqIO module產生Bio::Seq objects。

然後，你需要閱讀fasta-entry-numbers和希望進入哈希的位置。以fasta-name爲關鍵字，值爲要提取的每個子序列的兩個值的數組。如果每個fasta條目可能有多個這樣的子序列，則必須創建一個數組數組作爲每個鍵的值條目。

有了這個數據結構，你可以繼續使用subseq method from Bio::Seq來提取序列。

我希望這是一種適合你的方式，但我確信這對純粹的bash來說也是可行的。

來源

2013-05-30 13:11:24 thunk

+1：不要重新發明輪子，特別是如果有人免費贈送輪子 – msw

感謝您的提示，我會看看它。我的背景是R編程，所以這使得事情變得目前有點難以改變思維方式......但是我認爲也值得熟悉Bioperl。 –

這不是一個答案，它是一個嘗試澄清你的問題;請讓我知道我的任務是否正確。

foreach row in blast: 
    get the proper (blast[$1]) sequence from fasta 
    drop bases (blast[$7..$8]) from sequence 
    print blast[$1], shortened_sequence

如果我有你的任務是正確的，你正在你的編程語言（慶典）和數據的特殊格式（跨行記錄拆分）步履蹣跚。 Perl或Python會更適合這項任務;的確是Perl的部分原因，因爲如果不是不可能的話，awk中的多個文件訪問真的很難。

你已經與你所知道的工具相差甚遠，但看起來你正在達到他們方便表達的極限。

來源

2013-05-30 13:11:31 msw

是的，您的澄清幾乎在頭上打了，這就是我打算做的。事實上，我將其作爲bash腳本的學習過程，但我想我會按照你的提示轉換到Perl（或者Bioperl，感謝@thunk），並在那裏盡我所能。 –

更新了答案：

awk ' 
NR==FNR && NF { 
    id=substr($1,2) 
    getline seq 
    a[id]=seq 
    next 
} 
($1 in a) && NF { 
    x=substr(a[$1],$7,$8) 
    sub(x, "", a[$1]) 
    print ">"$1"\n"a[$1] 
} ' 1.fasta 1.blast

來源

2013-05-30 13:31:44

@DanielFischer再次更新了答案。想獲得您的反饋。 –

感謝您的幫助，我仍然忙於理解@Birei的答案，因此我還沒有回答，目前我沒有我的數據，但星期一我將運行數據上的代碼。 –

謝謝，這個工程以及@Birei的解決方案，所以我想接受這兩個，但這不會工作不幸。但是這兩個答案都有助於更好地理解awk。 –

提取行和子上的另一個文件的信息

回答

相關問題