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）加上序列，然后从序列本身中提取除了文件之间和文件中的位置之外的所有位置，$7这意味着从前两个匹配输出将是$81.fasta

>1
ACTAGCTACGACTCGGACTGACGAGCTACGACTACGG
>27620
GTAGATAGAGATAGAGAGAGAGAGGGGGGAGA
...

（请注意，前三个条目>1不在此序列中）

ID 是连续的，这意味着我可以像这样提取所需的信息：

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

这给了我一个矩阵，它在第一列中包含正确的序列标识符行，在第二列中包含正确的序列行（= 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)从给定字符串中删除这两个位置之间的所有内容的命令，但我认为我可以在自己的脚本中执行此操作。特别是正确的管道对我来说是个问题。

4

4 回答 4

2

如果您从事生物信息学并使用 DNA 序列（甚至更复杂的东西，如序列注释），我建议您看看Bioperl。这显然需要 Perl 知识，但具有相当多的功能。

在您的情况下，您可能希望使用模块从您的 fasta 文件生成Bio::Seq对象。Bio::SeqIO

然后，您需要将所需的 fasta-entry-numbers 和位置读取到哈希中。以 fasta-name 作为键，值是要提取的每个子序列的两个值的数组。如果每个 fasta 条目可以有多个这样的子序列，则必须创建一个数组数组作为每个键的值条目。

使用此数据结构，您可以继续使用 from 的subseq方法Bio::Seq提取序列。

我希望这是一种适合您的方法，尽管我确信这对于纯 bash 也是可行的。

于 2013-05-30T13:11:24.100 回答
2

这不是答案，而是试图澄清您的问题；如果我正确理解了您的任务性质，请告诉我。

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

如果我的任务正确，那么您的编程语言 (bash) 和数据的特殊格式（跨行拆分的记录）会阻碍您。Perl 或 Python 会更适合这项任务；确实 Perl 的编写部分是因为当时的多个文件访问awk即使不是不可能也非常困难。

你已经用你知道的工具走了很远，但看起来你已经达到了它们方便表达的极限。

于 2013-05-30T13:11:31.790 回答
1

更新了答案：

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-30T13:31:44.427 回答
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 序列不超过一行。

于 2013-05-30T13:34:03.553 回答