☰

保留范围最大的行

tag-icon • tag-icon text-processing bioinformatics
保留范围最大的行

我有一个文件（fileA），如下所示：

>ENST00000614578.1 gene=WASP12 CDS=1-526
>ENST00000581055.1 gene=PTP2 CDS=138-579
>ENST00000577541.1 gene=PTP2 CDS=1-81
>ENST00000423485.5 gene=PTP2 CDS=160-4752
>ENST00000367142.4 gene=PRPF40B CDS=304-1032
>ENST00000366955.7 gene=WASP12 CDS=169-9511

我想只保留 CDS 值范围最大的行。该范围由第三列中的 CDS 值给出。例如，第一行的范围是 525 (526-1)，第二行的范围是 441 (579-138)

期望的输出：

>ENST00000423485.5 gene=PTP2 CDS=160-4752
>ENST00000367142.4 gene=PRPF40B CDS=304-1032
>ENST00000366955.7 gene=WASP12 CDS=169-9511

我尝试通过排序

sort -nrk3,3  fileA

但我不认为这是可行的方法，有什么建议吗？过滤标准是对具有相同基因名称的行取范围的最大值

答案1

鉴于新信息，我来到了这里

awk -F'[ =-]' '{ print $0" "$6-$5 | "sort -k4nr" }' fileA | \
cut -d' ' -f1-3 | \
awk -F'[ =]' '!seen[$3]++'

输出：

>ENST00000366955.7 gene=WASP12 CDS=169-9511
>ENST00000423485.5 gene=PTP2 CDS=160-4752
>ENST00000367142.4 gene=PRPF40B CDS=304-1032

这当然可以用一个命令来完成awk，但我仍在学习如何使用它。

OP编辑之前的先前答案：

给定一个最小值，假设为 700（以匹配您的示例输出），您可以尝试以下操作：

awk -F'[=-]' '$4-$3 > 700' fileA

输出：

>ENST00000423485.5 gene=PTP2 CDS=160-4752
>ENST00000367142.4 gene=PRPF40B CDS=304-1032
>ENST00000366955.7 gene=WASP12 CDS=169-9511

-F'[=-]'：使用=和-作为列分隔符，这会导致第 3 列和第 4 列：

1 526
138 579
1 81
160 4752
304 1032
169 9511

$4-$3 > 700：选择第 4 列减去第 3 列大于 700 的行。

答案2

假设范围始终呈现在正链上并且从不以相反方向呈现，并且假设原始文件中的列由单个空格分隔，

$ awk -F '[ =-]' '{ k = $3; r = $6 - $5 } (m[k] == "" || m[k] < r) { d[k] = $0; m[k] = r } END { for (k in d) print d[k] }' file
>ENST00000366955.7 gene=WASP12 CDS=169-9511
>ENST00000367142.4 gene=PRPF40B CDS=304-1032
>ENST00000423485.5 gene=PTP2 CDS=160-4752

这将数据视为具有由空格=和分隔的字段的行-。因此，范围的长度可以通过第 6 个字段减去第 5 个字段来找到。任何基因名称的最大范围保存在m（“最大”）数组中，相应的行保存在d（“数据”）数组中。

最后将采集到的数据d输出。

awk格式更好的代码：

BEGIN { FS = "[ =-]" }

{
    k = $3
    r = $6 - $5
}

(m[k] == "" || m[k] < r) {
    d[k] = $0
    m[k] = r
}

END {
    for (k in d) print d[k]
}

相关内容