代码之家 › 专栏 › 技术社区 › C. John

整数编码VCF文件的最优解

vcf-variant-call-format bioinformatics awk bash

0

C. John · 技术社区 · 4 年前

CHROM    POS REF   ALT   Geno      value                                                                                       
Chr16 616504 T     C     X93.968   0|1:7,28:35:99:0|1:616504_T_C:787,0,177:616504   
Chr16 616504 T     C     BESC.1    0/0:48,0:48:99:.:.:0,114,1710:.                  
Chr16 616504 T     C     BESC.10   1|1:0,23:23:72:1|1:616504_T_C:1059,72,0:616504   
Chr16 616504 T     C     BESC.100  0/0:34,0:34:96:.:.:0,96,1440:.                   
Chr16 616504 T     C     BESC.1001 0/0:47,0:47:99:.:.:0,120,1800:.                  
Chr16 616504 T     C     BESC.1002 0/0:39,0:39:99:.:.:0,108,948:.

我们的目标是把第一个和第三个字符从 value 列并求和,然后输出一个类似的文件,其中值列替换为该和。前两行的输出示例:

CHROM    POS REF   ALT   Geno      value                                                                                       
Chr16 616504 T     C     X93.968   1   
Chr16 616504 T     C     BESC.1    0

以下是我当前的解决方案,其中STDIN 1是输入文件名,STDIN 2是输出文件名:

#!/bin/bash
i=0
len=$(cat $1 | wc -l)

touch $2
while read -r line; do
    let "i++"
    geno=$(echo "$line" | awk '{ print $6 }'| cut -c1,3 | fold -w1 | paste -sd+ - | bc)

    echo "$line" | awk -v g="$geno" '{ $6=""; print $0 " " g}' >> $2

    echo "Processed " "$i/$len"

done < $1

实际输入文件有1707993个条目。根据我的解决方案,这大约需要4.5小时来计算。理想情况下,我可以在一个小时内运行这个,但我不确定这有多现实。谢谢

1 回复 | 直到 4 年前

1

3

Nick ODell 4 年前

这个循环会很慢,因为它在输入文件的每行创建和销毁6个外部进程。

在Bash中,调用命令来处理整个文件通常比在该文件的每行调用一次命令要快得多。

例如,如果您想在Awk中进行文件处理,可以这样做:

#!/bin/bash
set -eu

awk -v out_file="$2" 'NR == 1 {
    print > out_file;
}
NR != 1 {
    $6 = substr($6, 1, 1) +  substr($6, 3, 1);
    print > out_file;
    print "Processed " NR;
}' "$1"

说明:

它取第6个字段,取第一个和第三个字符并求和。然后,它将它们重新分配回字段6。在awk中,这意味着 print 将在更改一个字段的情况下发出原始输出。
out_file .
发出一个进度指示器。