问题描述
在分析的过程中,有些数据的染色体命名为“chr1、chr2、…、chrX、chrY”,而有些数据的染色体命名则为“1、2、…、X、Y” (也就是不包含 chr 字符)。这里,通过代码对 bam 文件作为修改,实现染色体名的统一。
代码实现
假设我们有一个名为
test.bam
的文件,其中染色体名不包含chr
字符,需要在染色体名前加上chr
字符。
通过 samtools
和 shell
实现 (注:samtools reheader 需要给一个-
的参数,不给会报错):
samtools view -H test.bam | sed -e 's/SN:\([0-9XY]\)/SN:chr\1/' -e 's/SN:MT/SN:chrM/' | samtools reheader - test.bam > test.CHR.bam
代码封装
因为会经常碰到这样的情况,因此就将上面的这段代码封装到一个名为 bam_add_chr.sh
的脚本,放在 bin
目录下面,方便调用。
#! /usr/bin/bash
samtools view -H $1 | sed -e 's/SN:\([0-9XY]\)/SN:chr\1/' -e 's/SN:MT/SN:chrM/' | samtools reheader - $1 > $2
echo "Finished!"
调用方法:
bam_add_chr.sh test.bam test.CHR.bam
其它方法
可以通过 Python 的 Pysam 模块进行修改,但计算速度相对更低。