从单个字符串高效地计算统计信息(通过bowtie2实现bam文件)

Question

我的目标是有效地将包含短DNA序列读取的bowtie2映射的bam文件转换为包含映射开始和百分比标识的简单表。我即将完成这一任务，但我的解决方案非常缓慢，无法处理一个重要的例外。我会一步一步地解释情况：

我们从这样的字符串开始

    FCC5G2YACXX:5:1101:1224:2059#NNNNNNNN   97  genome  96003934    24  118M4D11M   =   96004135    0   GCA....ACG  P\..GW^EO   AS:i:-28    XN:i:0  XM:i:2  XO:i:1  XG:i:4  NM:i:6  MD:Z:54G53T9^TACA11 YT:Z:UP

我们只取第四列和MD:Z:54G53T9^TACA11部分，后者表示匹配(数字)和不匹配(字母)。除非字母前面有“^”，否则字母在引用序列中是删除的。

然后，我将它计算成一个类似于这个96003934 98.00的字符串

其中第一列与原始输入的第四列相同，第二列是匹配，除以和匹配和不匹配。

由于我更喜欢bash/zsh脚本，所以我做了以下工作

    if_sam2tab() {
    tab=$(echo $1 | rev | cut -d '.' -f 2- | rev)
    if      [ ! -e ./bowtie_results/$tab.tab ]
    then    echo -e "rstart\tmatch\tmismatch" > ./bowtie_results/$tab.tab
            while   read -r l
            do      rstart=$(echo $l | cut -f 4 -d " "      )
                    t=$(echo $l | grep -o 'MD:Z:[0-9A-Z]*' )
                    match=$(echo ${t: 5} | tr '[a-zA-Z]' '+' | bc )
                    mismatch=$(echo ${t: 5} | tr -d '[0-9]\n' | wc -c )
                    sum=$(echo "$match + $mismatch" | bc )
                    id=$(echo "scale=2; $match / $sum *100"         | bc  )
                    echo -e "$rstart\t$id" >> ./bowtie_results/$tab.tab
            done < <(grep 'MD:Z' ./bowtie_results/$1 )
    fi
    }

然而，这个解决方案是错误的，因为它认为删除(如^TCTAAG )是不匹配的。

其次，对我来说，这似乎非常缓慢，即使我在六个cpu上并行运行这些函数中的五个。

概括地说，我的目标是有效地将带有映射信息的字符串转换为标识百分比。

谢谢你的关心

Answer 1

何塞先前回答的一个补充: MD字符串不一定是第18列。因此，我在Jose的awk脚本中添加了一行，以查找MD字符串，而不管在行中的位置如何。毫无疑问，在这个过程中增加一个regex步骤将减缓整个过程。

awk '{
 match($0, /MD:Z:[0-9A-Z\^]*/,m );
 split(m[0],v,/[\^:]/);
 nmatch = split(v[3],vmatch, /[^0-9]/);
 cmatch=0;
 for(i=1; i<=nmatch; i++) cmatch+=vmatch[i];
 printf("%s"OFS"%.2f\n", $4, cmatch*100/(cmatch+nmatch-1));
}' 

我选择直接在这个awk脚本中输出管道，而不把它保存在磁盘上。在我有限的测试中，我没有发现磁盘i/o受到限制，也没有发现将bowtie2输出直接放入awk，在同时运行时会大大降低总体性能。