对 bed 文件中的每个区间计算平均覆盖深度:
bedtools coverage -a BED -b [BAM] -mean > MeanCoverageBED.bedgraph
如果想要计算总体的覆盖深度,可以使用下面的管道
bedtools coverage -a BED -b [BAM] -mean | awk '{total+=$4} END {print total/NR}'
一些其他有用的命令
输出 bam 文件所有区间的覆盖度
bedtools genomecov -ibam [BAM] -bga -split > CoverageTotal.bedgraph
输出 bed 文件区间中每个碱基的读段深度
bedtools coverage -a [BED] -b [BAM] -d > PerBaseDepthBED.bedgraph
获取覆盖的区间占比
zero=$(bedtools genomecov -ibam [BAM] -g [ref.fasta] -bga | awk '$4==0 {bpCountZero+=($3-$2)} {print bpCountZero}' | tail -1)
nonzero=$(bedtools genomecov -ibam [BAM] -g [ref.fasta] -bga | awk '$4>0 {bpCountNonZero+=($3-$2)} {print bpCountNonZero}' | tail -1)
percent=$(bc <<< "scale=6; ($nonzero / ($zero + $nonzero))*100")
网友评论