DamID-seq分析（三）

质控分析：
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

image.png

（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

---

小结结束，给自己撒花~

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