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

问题描述

给出以下文件：

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的坐标（从低到高的数值），以获得以下文件：

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

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

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

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

解决方法

$ 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

,

$ awk 'p9!=$9{if(p0) print p0} !a[$9]++; {p9=$9; p0=$0} END{print p0}' file | 
  awk 'NR%2{k=($7=="+")?4:5; v=$k; next} {$k=v}1'

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 . - 1 ID=cds-XP_033376609.1
NW_022983512.1 RefSeq CDS 2650 3267 . + 0 ID=cds-XP_033376596.1

两个单独的脚本将简化逻辑，第一个脚本打印每个键的第一行和最后一行（如果存在一行，则重复）。第二个脚本根据符号选择正确的值。

,

另一个awk（还添加了$ 0的打印零，如注释所示）

> cat tst.awk
$9 == prev {
    $keep = val
    $8 = 0
    row = $0
    next
}

{ 
    print row
    prev = $9
    $8 = 0
    row = $0
    keep = ( $7=="+"? 4: 5 )
    val = $keep
}

END {
    print row
}

输出：

> 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