阅读:5805回复:0
全外显子分析(WES)的一般步骤-7 初级统计篇
这步主要目的是对比对情况,做一个统计。
#Make a temp directory mkdir -p ${tmp_folder}_dept #Running GATK depth of coverage java -Xmx${heap}m -Djava.io.tmpdir\=${tmp_folder}_dept \ -jar $gatk \ -T DepthOfCoverage \ -L $ExonFile \ -l INFO \ -R $REF \ -I $PWDS/${subjectID}.realigned.recal.bam \ -o $PWDS/${subjectID}.coverage.dept \ -L $ExonFile #Remove the temp directory rm -rf ${tmp_folder}_dept #Make a coverage directory and move all output to this folder mkdir -p $PWDS/coverage mv $PWDS/${subjectID}.coverage.dept* $PWDS/coverage/ #Calculate coverage for all subjects for i in `cat list.txt` do echo -ne "$i\t" >> coverage.txt cat ../$i/analysis/coverage_normdup/$i.coverage_normdup.dept.sample_cumulative_coverage_proportions \ | cut -f 8,10,12,22 >> coverage.txt cat ../$i/analysis/coverage/$i.coverage.dept.sample_cumulative_coverage_proportions \ | cut -f 8,10,12,22 >> coverage.txt done perl -0777 -i -pe "s/\ngte_6\tgte_8\tgte_10\tgte_20\n/\t/g" coverage.txt perl -0777 -i -pe "s/gte_6\tgte_8\tgte_10\tgte_20\n//g" coverage.txt mv coverage.txt tmp1 cat tmp1 | sort -k 2 -n -r > coverage.txt rm -f tmp1 最终生成这样子的一个统计结果:(xls表格,可以下载附件看一下,放在这里不美) 以及region.tsz,打开看就能知道了,是具体哪些区域覆盖到了,哪些没有。覆盖率多少左右 生信人员可以根据这几个文件,写脚本,画图,或者做更多的事,比如统计20X以下的区域有多少 同时未覆盖到的区域也会生成一个uncover.bed文件,告诉你哪些区域没有覆盖到。 [Total] Raw Reads (All reads) 830433 [Total] QC Fail reads 0 [Total] Raw Data(Mb) 234.3 [Total] Paired Reads 830433 [Total] Mapped Reads 745054 [Total] Fraction of Mapped Reads 89.72% [Total] Mapped Data(Mb) 218.77 [Total] Fraction of Mapped Data(Mb) 93.37% [Total] Properly paired 707868 [Total] Fraction of Properly paired 85.24% [Total] Read and mate paired 728720 [Total] Fraction of Read and mate paired 87.75% [Total] Singletons 16334 [Total] Read and mate map to diff chr 12783 [Total] Read1 412726 [Total] Read2 417707 [Total] Read1(rmdup) 373502 [Total] Read2(rmdup) 358845 [Total] forward strand reads 378345 [Total] backward strand reads 366709 [Total] PCR duplicate reads 0 [Total] Fraction of PCR duplicate reads 0.00% [Total] Map quality cutoff value 20 [Total] MapQuality above cutoff reads 727973 [Total] Fraction of MapQ reads in all reads 87.66% [Total] Fraction of MapQ reads in mapped reads 97.71% [Target] Target Reads 448085 [Target] Fraction of Target Reads in all reads 53.96% [Target] Fraction of Target Reads in mapped reads 60.14% [Target] Target Data(Mb) 115.14 [Target] Target Data Rmdup(Mb) 115.11 [Target] Fraction of Target Data in all data 49.14% [Target] Fraction of Target Data in mapped data 52.63% [Target] Len of region 30297 [Target] Average depth 3800.31 [Target] Average depth(rmdup) 3799.28 [Target] Coverage (>0x) 100.00% [Target] Coverage (>=4x) 100.00% [Target] Coverage (>=10x) 100.00% [Target] Coverage (>=30x) 100.00% [Target] Coverage (>=100x) 100.00% [Target] Target Region Count 86 [Target] Region covered > 0x 86 [Target] Fraction Region covered > 0x 100.00% [Target] Fraction Region covered >= 4x 100.00% [Target] Fraction Region covered >= 10x 100.00% [Target] Fraction Region covered >= 30x 100.00% [Target] Fraction Region covered >= 100x 100.00% [flank] flank size 200 [flank] Len of region (not include target region) 33191 [flank] Average depth 363.25 [flank] flank Reads 405120 [flank] Fraction of flank Reads in all reads 48.78% [flank] Fraction of flank Reads in mapped reads 54.37% [flank] flank Data(Mb) 12.06 [flank] Fraction of flank Data in all data 5.15% [flank] Fraction of flank Data in mapped data 5.51% [flank] Coverage (>0x) 14.60% [flank] Coverage (>=4x) 14.35% [flank] Coverage (>=10x) 14.15% [flank] Coverage (>=30x) 14.11% [flank] Coverage (>=100x) 14.11% |
|