DamID-seq分析（三）---liu_ll

质控分析： https://www.jianshu.com/p/e0bfa7b75ccb

mapping: https://www.jianshu.com/writer#/notebooks/30869589/notes/66543025

接下来，我需要去提取mapping完后的bam文件里的GATC的具体位置信息。

想法:

1：把bam的一对reads里面的chr start end提取出来 2：因为DpnI这个酶切的GATC中间的位点，为了还原其原来的信息，所以对start-2 end+2处理

bedtools bamtobed -i CTCF2_q10.tmpp.bam > CTCF2_q10.pos.bed #get chr_start chr_end cat CTCF2_q10.pos.bed|awk '{print"chr"$1"t"$2-2"t"$2+2"n"}' > chr_start.bed cat CTCF2_q10.pos.bed|awk '{print"chr"$1"t"$3-2"t"$3+2"n"}' > chr_end.bed cat chr_end.bed chr_start.bed > chr_star_end.bed

其实还可以用命令行输出bedpe格式，然后再取bedpe里面的坐标（老师推荐~）

全基因组GATC的坐标，我从网上找到了一个脚本，具体的下载地址和用法见下面链接。我们可以直接输入参考基因组序列，然后就可以直接得到结果。 github:GATC_COUNT

（PS：这个分析 GATC 位点的脚本属于这个作者自己写的 DamID 的分析脚本其中之一，这个脚本非常友好，对小白很 nice ，前提是如果把 DamID 当做 ChIP-seq 的替代品，看一下自己的实验目的）

计算 DamID per GATC 位点的值 一：把 bam 转成 bigwig ，然后得到的结果就每一个 motif 上面的值

bigWigAverageOverBed CTCF.bigWig GATC.wholegenome.bed -bedOut=CTCF.bed out.tab

二：用 bedtools intersect + groupby

bedtools intersect -a GATC.wholegenome.bed -b CTCF.chr_start_end.bed -wao |bedtools sort -i > out1.bed bedtools groupby -i out.bed -g 1,2,3 -c 5 -o

得到的结果是每个 GATC 位点上有多少个 reads 落在上面

接下来就去 normalize： 用每个位点的 reads 出一下 reads 的总和（因为 GATC 都是 4bp 所以不用考虑长度的问题） 最后得到的是 DamID 的结果！因为我的 DamID 的测序量不够，所以我画 bin 来计算，计算一个 bin 里面的 Dam-CTCF/Dam 的信号：

10k_bin

Ref: https://github.com/Vift/DamID-Seq

小结结束，给自己撒花~

最近算东西真的是算的头发昏。。。。。继续加油呀~

---来自腾讯云社区的---liu_ll