合并基因列表与统计文件时首个基因丢失问题及代码修改需求

现有文件内容 #

基因列表文件（gene_list.txt） #

Gene
ABC
ASF
AGH

统计数据文件（{sample}.stats） #

Gene    total_coverage  average_coverage    S12_total_cvg   S12_mean_cvg    S12_granular_Q1 S12_granular_median S12_granular_Q3 S12_%_above_10  S12_%_above_20  S12_%_above_50
ABC 68264   74.36   68264   74.36   51  68  92  100.0   99.8    76.0
ASF 653934  232.14  653934  232.14  178 216 265 100.0   100.0   100.0
AGH 653934  232.14  653934  232.14  178 216 265 100.0   100.0   100.0
ORD 653934  232.14  653934  232.14  178 216 265 100.0   100.0   100.0
NOC 559425  248.63  559425  248.63  174 246 316 100.0   100.0   100.0

问题现象 #

执行代码后，输出结果中首字母排序的第一个基因（ABC）脱离表头单独显示：

ABC
Gene    total_coverage  average_coverage    S12_%_above_20  S12_%_above_50
ASF 653934  232.14  100.0   100.0
AGH 653934  232.14  100.0   100.0

使用的代码 #

rule join:
    input:
        stat="{sample}.stats",
        genes="gene_list.txt"
    output: 
        "{sample}.stats.output"
    shell:
        """
        join  --header -1 1 -2 1 -t $'\t' <(awk 'NR==1; NR > 1 {{print $0 | "sort -k1,1"}}' {input.stat}) <(sort -k1,1 {input.genes}) | cut -d$'\t' -f1-3,10-11 >> {output}
        """

期望正确输出 #

Gene    total_coverage  average_coverage    S12_%_above_20  S12_%_above_50
ABC 68264   74.36   99.8    76.0
ASF 653934  232.14  100.0   100.0
AGH 653934  232.14  100.0   100.0

原因分析 #

问题出在awk处理统计文件的逻辑：NR==1; NR > 1 {{print $0 | "sort -k1,1"}}中，表头行（NR==1）会直接输出，而后续行通过管道交给sort命令处理。但sort是缓冲输出的，导致表头先被join读取，排序后的内容（含ABC行）后输出。同时基因列表排序后，表头Gene与统计文件表头匹配，ABC作为第一个排序后的基因，会在表头匹配完成前被join处理，最终脱离表头单独显示。

解决方法 #

修改awk逻辑，先收集所有行再统一处理，确保表头和排序后的内容输出时序一致：

rule join:
    input:
        stat="{sample}.stats",
        genes="gene_list.txt"
    output: 
        "{sample}.stats.output"
    shell:
        """
        join --header -1 1 -2 1 -t $'\t' <(awk '{{lines[NR]=$0}} END {{print lines[1]; for(i=2;i<=NR;i++) print lines[i] | "sort -k1,1"}}' {input.stat}) <(sort -k1,1 {input.genes}) | cut -d$'\t' -f1-3,10-11 > {output}
        """

修改说明 #

1. 用awk数组lines存储所有行，避免表头与内容输出时序不一致的问题；
2. 先输出表头行，再将剩余行排序后输出，保证统计文件是「表头+排序后基因数据」的结构；
3. 将>>改为>，避免重复运行时输出文件内容累加。

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