#! /bin/bash

#purpose:bam2bedGraph统计每个碱基的depth，并利用bedtools求overlapped region。
#set env
bam_dir=" "
cd $bam_dir  #bam文件所在位置
ls *.sorted.dedup.recal.bam > bam.list
probes_bed="truseq-exome-targeted-regions-manifest-v1-2.bed2" 

for sample in $(cat bam.list)
do
   bedtools genomecov  -bg -ibam   $bam_dir/$sample > $sample.bam2bed.bedGraph  # generate $sample.bam2bed.bedGraph file
  #chr start end depth #zero-based start and a one-based end
  
  cat $sample.bam2bed.bedGraph | awk '{print $1"\t"$2+1"\t"$3"\t"$4}' >$sample.bedGraph
  bedtools intersect -a $sample.bedGraph -b $probes_bed > $sample.intersect
  #chr overlapped_start overlapped_end depth 
  echo $sample ok
done

#cal the length and depth of overlapped regions
for sample in $(cat bam_list)
do
  perl -alne '{$len=$F[2]-$F[1]+1; $tmp+=$len}  END{print $tmp}'        $sample.intersect >> intersect.length
  perl -alne '{$len=$F[3]*($F[2]-$F[1]+1); $tmp+=$len}  END{print $tmp}' $sample.intersect >> intersect.depth
  perl -alne '{$len=$F[2]-$F[1]+1;$tmp+=$len}END{print $tmp}' $probes_bed  >> bed_length
done
paste bam.list intersect.length intersect.depth bed_length > intersect.cov_depth
#intersect.cov_depth: sampleID intersect.length intersect.depth bed_length 
#coverage=intersect.length/bed_length
#depth=intersect.depth/intersect.length

awk '{print $2/$4}' intersect.cov_depth #depth
awk '{print $3/$2}' intersect.cov_depth #coverage