从包含多组基因坐标的文件中提取一组基因坐标 [英] Extract One Set of Gene Coordinates From a File Containing Several Sets of Gene Coordinates

问题描述

给出以下文件:

NW_022983499.1  RefSeq  CDS 6883    7503    .   +   0   ID=cds-XP_033376633.1
NW_022983500.1  RefSeq  CDS 5353    5898    .   +   0   ID=cds-XP_033376630.1
NW_022983500.1  RefSeq  CDS 6033    7994    .   +   0   ID=cds-XP_033376630.1
NW_022983502.1  RefSeq  CDS 5391    5543    .   +   0   ID=cds-XP_033376626.1
NW_022983502.1  RefSeq  CDS 5591    5673    .   +   0   ID=cds-XP_033376626.1
NW_022983502.1  RefSeq  CDS 5782    5895    .   +   1   ID=cds-XP_033376626.1
NW_022983502.1  RefSeq  CDS 5937    6424    .   +   1   ID=cds-XP_033376626.1
NW_022983502.1  RefSeq  CDS 6478    6680    .   +   2   ID=cds-XP_033376626.1
NW_022983502.1  RefSeq  CDS 6739    6858    .   +   0   ID=cds-XP_033376626.1
NW_022983502.1  RefSeq  CDS 6926    7408    .   +   0   ID=cds-XP_033376626.1
NW_022983504.1  RefSeq  CDS 5478    5513    .   -   0   ID=cds-XP_033376620.1
NW_022983504.1  RefSeq  CDS 5353    5419    .   -   0   ID=cds-XP_033376620.1
NW_022983504.1  RefSeq  CDS 5161    5297    .   -   2   ID=cds-XP_033376620.1
NW_022983504.1  RefSeq  CDS 5059    5115    .   -   0   ID=cds-XP_033376620.1
NW_022983508.1  RefSeq  CDS 4415    5392    .   -   1   ID=cds-XP_033376609.1
NW_022983508.1  RefSeq  CDS 4215    4344    .   -   1   ID=cds-XP_033376609.1
NW_022983512.1  RefSeq  CDS 2650    2831    .   +   0   ID=cds-XP_033376596.1
NW_022983512.1  RefSeq  CDS 2890    3112    .   +   1   ID=cds-XP_033376596.1
NW_022983512.1  RefSeq  CDS 3163    3267    .   +   0   ID=cds-XP_033376596.1

我想提取一组对应于第9列中ID的坐标(从低到高的数值)，以获得以下文件:

I would like to extract one set of coordinates (from lower to higher numerical value) corresponding to the IDs present in column 9, so as to obtain the following file:

NW_022983499.1  RefSeq  CDS 6883    7503    .   +   0   ID=cds-XP_033376633.1
NW_022983500.1  RefSeq  CDS 5353    7994    .   +   0   ID=cds-XP_033376630.1
NW_022983502.1  RefSeq  CDS 5391    7408    .   +   0   ID=cds-XP_033376626.1
NW_022983504.1  RefSeq  CDS 5059    5513    .   -   0   ID=cds-XP_033376620.1
NW_022983508.1  RefSeq  CDS 4215    5392    .   -   0   ID=cds-XP_033376609.1
NW_022983512.1  RefSeq  CDS 2650    3267    .   +   0   ID=cds-XP_033376596.1

请注意，对于在第7列中具有正值的ID=cds-XP_033376630.1，我需要选择第2行第4列5353和第3行第5列7994的值.

Note that in the case of ID=cds-XP_033376630.1 who has a positive value in column 7, I need to select the value of line 2 column 4 5353 and of line 3 column 5 7994.

相反，如果第7列的值为负，如ID=cds-XP_033376620.1所示，则逻辑取反，我需要选择第14行，第4列5059和第11行，第5列

In contrast, if the value of column 7 is negative, as in ID=cds-XP_033376620.1, the logic gets inverted, I need to select the value of line 14, column 4 5059 and of line 11, column 5 5513

我对使用AWK(不是Perl或Python)来解决这种经典的生物信息学问题特别感兴趣，如果有人能指出正确的方向，我将不胜感激.

I am specially interested in using AWK (not Perl or Python), to solve this classical bioinformatics problem, and I would appreciate if anyone could point me in the right direction.

推荐答案

$ cat tst.awk
$NF != prevKey {
    if ( NR > 1 ) {
        prt()
    }
    min  = $4
    max  = $5
    line = $0
    prevKey = $NF
}
{
    min = ($4 <= min ? $4 : min)
    max = ($4 >= max ? $5 : max)
}
END { prt() }

function prt(   orig) {
    orig = $0
    $0 = line
    $4 = min
    $5 = max
    $8 = 0
    print
    $0 = orig
}

.

$ awk -f tst.awk file
NW_022983499.1 RefSeq CDS 6883 7503 . + 0 ID=cds-XP_033376633.1
NW_022983500.1 RefSeq CDS 5353 7994 . + 0 ID=cds-XP_033376630.1
NW_022983502.1 RefSeq CDS 5391 7408 . + 0 ID=cds-XP_033376626.1
NW_022983504.1 RefSeq CDS 5059 5513 . - 0 ID=cds-XP_033376620.1
NW_022983508.1 RefSeq CDS 4215 5392 . - 0 ID=cds-XP_033376609.1
NW_022983512.1 RefSeq CDS 2650 3267 . + 0 ID=cds-XP_033376596.1

这篇关于从包含多组基因坐标的文件中提取一组基因坐标的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！