从包含多组基因坐标的文件中提取一组基因坐标

问题描述：

给出以下文件:

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

从包含多组基因坐标的文件中提取一组基因坐标

相关推荐