批量处理上千个fastq的时候会产生上千个bam，因为各种原因会造成bam文件损毁或不完整，因此必须要能快速检测出错误的bam文件，重跑或找出根本问题。
处理方法：检测bam文件的完整度
for i in *.bam ;do (samtools quickcheck $i || echo $i error);done
根据基因的名字，找到对应的染色体上的位置。【简单】
cat gencode.v27.annotation.gtf | grep -v '#' | awk '{if($3=="gene")print $0}' | cut -f1,4,5,9 | sed 's/ /;/g' | sed 's/ /;/g' | sed 's/;;/;/g' | sed 's/"//g' | cut -d; -f1,2,3,9 | sed 's/;/ /g' > gene.chr.pos.hg19.bed
参考：hg19 | GRCh38 | mm10 | 基因组 | 功能区域 | 位置提取
How to convert chromosomal coordinates (Bedgraph) into gene symbols?
根据染色体的位置，找到对应的基因。【有点繁琐，核心就是用bedtools的交集功能】
cat chr.pos.bed | sed 's/ / /g' > chr.pos.feature.bed bedtools intersect -wa -wb -a gene.chr.pos.hg19.bed -b chr.pos.feature.bed | cut -f4,8 > chr.pos.to.gene.name.bed
查看reads比对BAM文件中指定区域的reads的数量【涉及到exon或非基因区域时非要重要】
# index first samtools index bam.link.2650/23C9-E16-F8_STAR_Aligned.sortedByCoord.out.bam # samtools samtools depth -a -r chr19:803565-803565 bam.link.2650/23C9-E16-F8_STAR_Aligned.sortedByCoord.out.bam
批量运行bam文件
# mv bam.link.2650/10c2-E11-A11_STAR_Aligned.sortedByCoord.out.bam error_bam/ ls bam.link.2650/10c2*.bam > 10c2.bam.list samtools depth -a -r chr19:803565-803565 -f 10c2.bam.list
w完整版
samtools depth -a -r chr19:803561-803561 -f IMR.bam.list >> count.txt && samtools depth -a -r chr19:803561-803561 -f UE.bam.list >> count.txt && samtools depth -a -r chr19:803561-803561 -f 1c11.bam.list >> count.txt && samtools depth -a -r chr19:803561-803561 -f 5c3.bam.list >> count.txt && samtools depth -a -r chr19:803561-803561 -f 10c2.bam.list >> count.txt && samtools depth -a -r chr19:803561-803561 -f 17c8.bam.list >> count.txt && samtools depth -a -r chr19:803561-803561 -f 20c7.bam.list >> count.txt && samtools depth -a -r chr19:803561-803561 -f 23c9.bam.list >> count.txt && echo "done!!!"
参考：samtools depth 用于外显子未覆盖区域的统计及统计未覆盖区域的意义
待续~