使用 awk 从 vcf 文件中提取字符串答案

【问题标题】：pull string from a vcf file using awk使用 awk 从 vcf 文件中提取字符串
【发布时间】：2019-12-04 00:05:49
【问题描述】：

我正在运行以下代码来操作 vcf 表中的数字数据。

 cat inputfile | while read row; do
                echo $row > tmp
                originalProb= `awk '{print $1}' tmp`
                probabilityHom1=`awk '{print $2}' tmp`
                probabilityHom2=`awk '{print $4}' tmp`
                numCols=`awk '{print NF}' tmp`

                if [ $numCols -gt 4 ]; then
                        echo "${originalProb}" >> currentRowGenotypes
                elif [ "$probabilityHom1" -gt "$probabilityHom2" ]; then
                        echo "1/1" >> currentRowGenotypes
                elif [ "$probabilityHom1" -lt "$probabilityHom2" ]; then
                        echo "0/0" >> currentRowGenotypes
                elif [ "$probabilityHom1" -eq "$probabilityHom2" ] && [ "$probabilityHom1" -eq 0 ]; then
                        echo "${originalProb}" >> currentRowGenotypes
                else                    
                        echo "het" >> currentRowGenotypes
                fi

        done

        cat tmpHeaders currentRowGenotypes > currentFullCol

输入文件如下所示

1/1     255     231     0
0/1     255     0       152
0/1     255     0       82
0/1     255     0       151
0/1     239     0       31
0/1     255     0       255

由于某种原因，awk 命令无法识别第一列。有什么建议吗？

【问题讨论】：

"originalProb= `awk '{print $1}' tmp`" - bash 可以识别空格，= 后面有一个空格。请使用$(...) 而不是反引号`
同上，但所有这些都可以是一个 awk 脚本。见grymoire.com/Unix/Awk.html。祝你好运。

标签： bash unix bioinformatics vcf-variant-call-format

【解决方案1】：

创建一个临时文件只是为了使awk 不是一个好主意将行拆分为列，因为：

逐行创建临时文件会产生开销。
它会多次生成子进程来调用awk。
bash 和awk 之间的语法差异可能是导致错误的原因。

您可以不使用awk。请尝试以下方法：

while read -ra row; do
    originalProb="${row[0]}"
    probabilityHom1="${row[1]}"
    probabilityHom2="${row[3]}"
    numCols="${#row}"

    if (( numCols > 4 )); then
        echo "$originalProb" >> currentRowGenotypes
    elif (( probabilityHom1 > probabilityHom2 )); then
        echo "1/1" >> currentRowGenotypes
    elif (( probabilityHom1 < probabilityHom2 )); then
        echo "0/0" >> currentRowGenotypes
    elif (( probabilityHom1 == probabilityHom2 &&  probabilityHom1 == 0 )); then
        echo "$originalProb" >> currentRowGenotypes
    else
        echo "het" >> currentRowGenotypes
    fi
done < inputfile

cat tmpHeaders currentRowGenotypes > currentFullCol

正如其他人反复建议的那样，更好的方法是使用awk：

awk '{
    originalProb = $1
    probabilityHom1 = $2
    probabilityHom2 = $4
    numCols = NF

    if ( numCols > 4 )
        print originalProb >> "currentRowGenotypes"
    else if ( probabilityHom1 > probabilityHom2 )
        print "1/1" >> "currentRowGenotypes"
    else if ( probabilityHom1 < probabilityHom2 )
        print "0/0" >> "currentRowGenotypes"
    else if ( probabilityHom1 == probabilityHom2 && probabilityHom1 == 0 )
        print originalProb >> "currentRowGenotypes"
    else
        print "het" >> "currentRowGenotypes"
}' inputfile

cat tmpHeaders currentRowGenotypes > currentFullCol

希望这会有所帮助。

【讨论】：

【解决方案2】：

为什么不使用Pysam？非常适合解析 BCF/VCF。

【讨论】：