||
#! /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
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 07:17
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社