☰

具有定义范围的组 ID

tag-icon • tag-icon text-processing bioinformatics
具有定义范围的组 ID

我有一个 ID 和数字（位置）的排序文件。我需要将第二列中的位置按 500 个间隔分组为一组。

如果该行的值与上一行相比小于500，则将它们分到同一组；如果该行的值大于 500，则将它们分到不同的组中。

输入文件：

snp00001    200
snp00002    300
snp00003    400
snp00004    500
snp00005    600
snp00006    900
snp00007    1500
snp00008    1800
snp00009    3000
snp00010    3500
snp00011    4000
snp00012    5000

所需输出

snp00001 200 Group1
snp00002 300 Group1
snp00003 400 Group1
snp00004 500 Group1
snp00005 600 Group1
snp00006 900 Group1
snp00007 1500 Group2
snp00008 1800 Group2
snp00009 3000 Group3
snp00010 3500 Group3
snp00011 4000 Group4
snp00012 5000 Group5

额外注意：snp00001到snp00006将被分为同一组，因为它们之间的范围（snp00002 - snp00001）或（snp00003 - snp00002）或（snp00004 - snp00003）...小于500。

snp00006 和 snp00007 被分到下一组，因为它们之间的范围 (snp00007 - snp00006) 大于 500。

我尝试过 awk，但没有成功。

awk -v step=500 -v OFS='\t' '{if(NR==1 || $2+limit){group++} file="Group"group; print file,$0}' input_file

答案1

您需要跟踪前一个值并将当前值与保存的值进行比较。如果差异超过 500，则增加组数。

例如

awk -v group=1 '{if ($2-prev>500) { group++ }} {prev=$2; $3="group" group; print}'
snp00001 200 group1
snp00002 300 group1
snp00003 400 group1
snp00004 500 group1
snp00005 600 group1
snp00006 900 group1
snp00007 1500 group2
snp00008 1800 group2
snp00009 3000 group3
snp00010 3500 group3
snp00011 4000 group3
snp00012 5000 group4

（FWIW，你的 9/10/11 输出不一致；9->10 是 500 但不增加组，但 10->11 也是 500 但确实增加组）。

答案2

使用乐（以前称为 Perl_6）

这是一种稍微不同的分组方案，可能会被证明是有用的。它使用 Raku 的~~smartmatch 运算符快速判断某个位置是否在某个范围内（或不在范围内）：

~$ raku -e 'my $i = 1; my $r = 1..500; for lines() {my $a = .words;  \
            if ($a.[1].Int ~~ $r) {say "$a Group", $i, " ", $r} else {  \
            repeat { $r+=500 } until ($a.[1].Int ~~ $r);  \
            say "$a Group", ++$i, " ", $r };}' file

输入示例：

snp00001    200
snp00002    300
snp00003    400
snp00004    500
snp00005    600
snp00006    900
snp00007    1500
snp00008    1800
snp00009    3000
snp00010    3500
snp00011    4000
snp00012    5000

样本输出（从核苷酸 1 开始，每 500 个核苷酸对 SNP 进行分组）：

snp00001 200 Group1 1..500
snp00002 300 Group1 1..500
snp00003 400 Group1 1..500
snp00004 500 Group1 1..500
snp00005 600 Group2 501..1000
snp00006 900 Group2 501..1000
snp00007 1500 Group3 1001..1500
snp00008 1800 Group4 1501..2000
snp00009 3000 Group5 2501..3000
snp00010 3500 Group6 3001..3500
snp00011 4000 Group7 3501..4000
snp00012 5000 Group8 4501..5000

上面的 Raku 代码声明了一个 Group# 迭代器$i和一个初始$r范围1..500.输入被视为lines，每行被分成 (空格分隔) words。运行 if/else 条件：if第二列在范围、行、组# 和范围~~内智能匹配，取范围和edly 递增，而 not（即）智能匹配成功。然后打印与之前相同的信息，但这次 Group# 正确递增 ( )。$rsayelse$rrepeat500until~~++$i

上述分组方案的优点是所得组的间隔长度均相等，在本例中约为 500 个核苷酸。该方案可防止组间隔长度的“扩张”，当多个 SNP 稍微位于一起时可能会发生这种情况（间隔“扩张”可能会造成“聚类”的错误印象）。

为了使其成为更通用的“分组”工具，您可以将 Range 的右侧抽象为变量 ( $m)，以进行快速分组：

~$ raku -e 'my $i=1; my $m=1000; my $r = 1..$m; for lines() {my $a = .words;   if ($a.[1].Int ~~ $r) {say "$a\tGroup$i\t", $r} else { repeat { $r+=$m } until ($a.[1].Int ~~ $r); say "$a\tGroup{++$i}\t", $r };}' file
snp00001 200    Group1  1..1000
snp00002 300    Group1  1..1000
snp00003 400    Group1  1..1000
snp00004 500    Group1  1..1000
snp00005 600    Group1  1..1000
snp00006 900    Group1  1..1000
snp00007 1500   Group2  1001..2000
snp00008 1800   Group2  1001..2000
snp00009 3000   Group3  2001..3000
snp00010 3500   Group4  3001..4000
snp00011 4000   Group4  3001..4000
snp00012 5000   Group5  4001..5000

https://docs.raku.org/type/Range
https://raku.org

相关内容