基于blastx/diamond结果计算参考序列水平覆盖的方案

将DNA序列reads比对至蛋白质数据库时，多数比对工具无法提供参考序列层面的覆盖度信息，但blastx/diamond可输出每条query reads在参考基因上的比对起始（Start）与终止（End）位置。因此需按参考基因ID统计多个[Start..End]区间内的唯一整数数量，以此得到比对序列的有效长度。

预处理步骤 #

先用以下命令处理blastx/diamond的输出结果，保留参考基因ID、Start、End三列，同时修正反向序列（确保Start ≤ End）：

awk -F"\t" '{if ($2 <= $3) {print $0} else {print $1"\t"$3"\t"$2}}'

注：输入数据中同一ID对应的区间可能存在重叠或重复，重复的整数仅需计数一次。

输入示例 #

ID  Start   End
A   1   50
A   2   45
A   25  150
A   50  150
A   155 200
A   205 300
B   5   50
B   61  70
B   81  100
C   1   500

期望输出示例 #

ID  count
A   292
B   76
C   500

Ruby 实现（@dawg） #

ruby -lane 'BEGIN{h=Hash.new { |hash, key| hash[key] = Set.new() }}
    h[$F[0]].merge(($F[1].to_i..$F[2].to_i)) if $.>1
END{
    puts "ID\tCount"
    h.each{|k,v| puts "#{k}\t#{v.length}"}
}
' Input.file > Output.file

Perl 实现（@zdim） #

perl -MData::Dumper -MList::Util=uniq -wnE'
    ($id, $beg, $end) = split; 
    next if not $beg or $beg =~ /[^0-9]/ or not $end or $end =~ /[^0-9]/;
    push @{$res{$id}}, $beg..$end; 
    }{ 
    for (keys %res) { $res{$_} = uniq sort { $a <=> $b } @{$res{$_}} };
    say Dumper \%res
' Input.file > Output.file

Awk 实现（@pmf） #

awk 'NR>1 {s[$1] += $3 - ($2 <= b[$1] ? ($3 <= b[$1] ? $3 : b[$1]) + 1 : $2) + 1; b[$1] = b[$1] <= $3 ? $3 : b[$1]} END {OFS="\t"; print "ID", "count"; for (i in s) {print i, s[i]}}' Input.file > Output.file

Awk 实现（@dawg） #

awk '
FNR>1{for(i=$2;i<=$3;i++) ss[$1 "|" i]}
END{
    print "ID", "Count"
    for (e in ss) {
        split(e,idx,"|")
        cnt[idx[1]]++
    }
    for (e in cnt) print e, cnt[e]
}
' OFS="\t" Input.file > Output.file

内容的提问来源于stack exchange，提问作者SWK