【ChIP-seq】对Bam 文件添加chr

参考：

Question: Bam File: Change Chromosome Notation

当我们分析ChIP-seq数据时候，经常会用到ENCODE数据，需要考虑所网站使用的参考基因组版本及其是否存在chr.等问题。如何对基因组添加或者删除chr呢？以bam 文件为例，对其进行添加chr 操作，将bam 文件中"1"变成"chr1".

image.png

实践

对bam 文件添加chr,画出profiler 图。

• Part1: 添加chr

 ##### 添加chr
$ samtools view -h ACT-07_rmdup_picard.bam | sed -e '/^@SQ/s/SN\:/SN\:chr/' -e '/^[^@]/s/\t/\tchr/2'|awk -F ' ' '$7=($7=="=" || $7=="*"?$7:sprintf("chr%s",$7))' |tr " " "\t" | samtools view -h -b -@ 10 -S - > ACT-07_rmdup_picard.chr.bam

• Part2: bam文件转成bw

i="ACT-07_rmdup_picard.chr.bam"

base=$(basename $i "_rmdup_picard.bam");
num1=10000000;
num2="$(samtools view -c  $i  2>&1 )";
res=$(printf "%.5f" `echo "scale=5;$num1/$num2"|bc`);
samtools index -@ 20 $i
bamCoverage -p 20 --scaleFactor  $res -b  $i   -o   ./${base}.rmdup.10M.bw;

• Part3：画profilter图

GeneRegion=/public/home/kcao/genome_human/gtf_human_hg38/Homo_sapiens.GRCh38.92.chr.bed6

for sample in *bw* ; do echo $sample;
computeMatrix reference-point -p 20  --referencePoint TSS -b 5000 -a 5000  -R $GeneRegion  -S  $sample --skipZeros -o ./2_TSS_profiler/${sample}.mat.gz && plotHeatmap -m ./2_TSS_profiler/${sample}.mat.gz -out ./2_TSS_profiler/${sample}.heatmap.png  &
done

## 画一起
computeMatrix reference-point -p 20  --referencePoint TSS -b 5000 -a 5000  -R $GeneRegion  -S ACT-07_rmdup_picard.chr.bam.rmdup.10M.bw ENCODE_merge.rmdup.10M.bw --skipZeros -o ./2_TSS_profiler/merge_TSS.mat.gz && plotHeatmap -m ./2_TSS_profiler/merge_TSS.mat.gz -out ./2_TSS_profiler/merge_TSS.heatmap.png &

image.png

思考：

• 使用deeptools 对多个文件进行一起绘图，推荐对bw 进行标准化，会影响最后画图效果。
  比如plotCorrelation 发现图中没有任何点，只有一个相关系数，这时候很可能是由于两个bw 文件，文库大小相差很大，需要先产生标准化后的bw ,再进行绘图。
• 需要注意不同工具是否要求chr.
  比如bedtools intersectBed 需要两个文件同时存在chr 或者没有chr.
  ChIPseeker 输入的peak 文件是否需要chr,需要考虑org.Hs.eg.db.

欢迎评论交流~