我是一个生物信息学专业的学生,刚接触到bash和编程的东西。我想计算基因组覆盖率。
这是我的脚本。我用xx
替换了真实的参数,但我确信xx
没有问题。其他学生已经执行了这个脚本,没有出现错误。
filename=$1
reference=/xx
filebase=$(basename $filename .bam)
samtools view ${filename} -F 4 -q 30 -b > ${filebase}.f.bam
genomeCoverageBed -ibam -g ${filebase}.f.bam ${reference} > /mnt/ABC/projects/abc/def/${filebase}.cov
coverage=$(grep genome /mnt/ABC/projects/abc/def/${filebase}.cov | awk '{NUM+=$2*$3; DEN+=$3} END {print NUM/DEN}')
echo -e "${filename},${coverage}" >> coverages.txt
当我用sh ./coverage.sh /mnt/XYZ/share/sdf_rawdata/hsa/mergedbams/ghj_merged_200203.hs37d5.cons.90perc.bam
它不起作用,并给出:awk: cmd. line:1: fatal: division by zero attempted and unrecognized parameter:/mnt/XYZ/share/sdf_rawdata/hsa/mergedbams/ghj_merged_200203.hs37d5.cons.90perc.bam
错误,并且在coverages.txt文件中只有这一行:-e /mnt/NEOGENE2/share/compevo_rawdata/hsa/mergedbams/Ash128_all.merged.hs37d5.fa.cons.90perc.bam
,仅此而已。
谢谢你的帮助。
1条答案
按热度按时间djp7away1#
您需要设置一个条件来检查
DEN
变量是否为NOT NULL,然后仅在awk
代码的END
块中执行除法(尝试在此处修复OP的尝试)。coverage=$(grep genome /mnt/ABC/projects/abc/def/${filebase}.cov | awk '{NUM+=$2*$3; DEN+=$3} END {if(DEN){print NUM/DEN}}')
你不需要把
grep
命令和awk
沿着使用,我们可以在awk
本身中搜索字符串,可能是这样的:coverage=$(awk '/genome/{NUM+=$2*$3; DEN+=$3} END {if(DEN){print NUM/DEN}}' "/mnt/ABC/projects/abc/def/${filebase}.cov")
***为什么会出现错误:***因为有时候变量
DEN
中的值为零。让我们以更短的方式举一个例子(只是一些例子来理解这里的错误):当变量
a
为NULL时,我们也会得到相同错误当变量
a
为零时,我们也会得到相同的错误: