收藏链接

批量处理上千个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 用于外显子未覆盖区域的统计及统计未覆盖区域的意义

---

待续~