【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.

欢迎评论交流~