這里是佳奧!
獲得了比對(duì)的bam文件先誉,生成了index的bai文件,生成了比對(duì)結(jié)果stat文件褐耳。
##可以先multiqc看一下結(jié)果
multiqc ./
合并bam文件和去除PCR重復(fù)沒(méi)有嚴(yán)格先后順序铃芦,是否必要根據(jù)文章決定
1 合并bam文件
需要回文章看
是同一個(gè)樣本的多個(gè)測(cè)序、還是多個(gè)生物學(xué)重復(fù)刃滓。(一個(gè)樣本測(cè)了兩次,還是兩次實(shí)驗(yàn)同一個(gè)樣本)
因?yàn)?strong>一個(gè)樣品分成了多個(gè)lane進(jìn)行測(cè)序晃危,所以在進(jìn)行peaks calling的時(shí)候老客,需要把bam進(jìn)行合并震叮。
## 如果不用循環(huán) 格式:輸出 兩個(gè)輸入
samtools merge control.merge.bam Control_1_trimmed.bam Control_2_trimmed.bam
## 通常我們用批處理
mkdir mergeBam
cd /home/kaoku/chipseq/mouse_project/align
ls *.bam| sed 's/_[0-9]_trimmed.bam//g' |sort -u | while read id; do samtools merge ../mergeBam/$id.merge.bam $id*.bam ; done
##合并以后
(chipseq) root 16:05:07 /home/kaoku/chipseq/mouse_project/mergeBam
$ ls -lh
總用量 9.6G
-rw-r--r-- 1 root root 836M 8月 11 15:57 Control.merge.bam
-rw-r--r-- 1 root root 1.3G 8月 11 15:58 H2Aub1.merge.bam
-rw-r--r-- 1 root root 1.6G 8月 11 15:59 H3K36me3.merge.bam
-rw-r--r-- 1 root root 1.3G 8月 11 16:00 Ring1B.merge.bam
-rw-r--r-- 1 root root 1.1G 8月 11 16:01 RNAPII_8WG16.merge.bam
-rw-r--r-- 1 root root 1.5G 8月 11 16:02 RNAPII_S2P.merge.bam
-rw-r--r-- 1 root root 1.4G 8月 11 16:03 RNAPII_S5P.merge.bam
-rw-r--r-- 1 root root 215M 8月 11 16:04 RNAPII_S5PRepeat.merge.bam
-rw-r--r-- 1 root root 713M 8月 11 16:04 RNAPII_S7P.merge.bam
2 去除PCR重復(fù)
##使用軟件samtools(后臺(tái)運(yùn)行苇瓣,top查看后臺(tái))
ls *merge.bam | while read id ; do ( samtools markdup -r $id $(basename $id ".bam").rmdup.bam & );done
##把去除重復(fù)的bam文件再次建立索引,查看比對(duì)效果
ls *.rmdup.bam | xargs -i samtools index {}
ls *.rmdup.bam | while read id ; do ( samtools flagstat $id > $(basename $id ".bam").stat & );done
##結(jié)果如下哲嘲,去除重復(fù)前后比較
(chipseq) root 20:58:44 /home/kaoku/chipseq/mouse_project/mergeBam
$ ls -lh
總用量 18G
-rw-r--r-- 1 root root 836M 8月 11 15:57 Control.merge.bam
-rw-r--r-- 1 root root 743M 8月 11 20:57 Control.merge.rmdup.bam
-rw-r--r-- 1 root root 1.3G 8月 11 15:58 H2Aub1.merge.bam
-rw-r--r-- 1 root root 1.1G 8月 11 20:58 H2Aub1.merge.rmdup.bam
-rw-r--r-- 1 root root 1.6G 8月 11 15:59 H3K36me3.merge.bam
-rw-r--r-- 1 root root 1.5G 8月 11 20:58 H3K36me3.merge.rmdup.bam
-rw-r--r-- 1 root root 1.3G 8月 11 16:00 Ring1B.merge.bam
-rw-r--r-- 1 root root 1006M 8月 11 20:58 Ring1B.merge.rmdup.bam
-rw-r--r-- 1 root root 1.1G 8月 11 16:01 RNAPII_8WG16.merge.bam
-rw-r--r-- 1 root root 984M 8月 11 20:58 RNAPII_8WG16.merge.rmdup.bam
-rw-r--r-- 1 root root 1.5G 8月 11 16:02 RNAPII_S2P.merge.bam
-rw-r--r-- 1 root root 1.2G 8月 11 20:58 RNAPII_S2P.merge.rmdup.bam
-rw-r--r-- 1 root root 1.4G 8月 11 16:03 RNAPII_S5P.merge.bam
-rw-r--r-- 1 root root 775M 8月 11 20:58 RNAPII_S5P.merge.rmdup.bam
-rw-r--r-- 1 root root 215M 8月 11 16:04 RNAPII_S5PRepeat.merge.bam
-rw-r--r-- 1 root root 210M 8月 11 20:57 RNAPII_S5PRepeat.merge.rmdup.bam
-rw-r--r-- 1 root root 713M 8月 11 16:04 RNAPII_S7P.merge.bam
-rw-r--r-- 1 root root 610M 8月 11 20:57 RNAPII_S7P.merge.rmdup.bam
##查看一下比對(duì)成功率
$ grep 'N/A' *.stat | grep '%'
Control.merge.rmdup.stat:12330969 + 0 mapped (85.16% : N/A)
H2Aub1.merge.rmdup.stat:17516222 + 0 mapped (96.82% : N/A)
H3K36me3.merge.rmdup.stat:22685679 + 0 mapped (98.51% : N/A)
Ring1B.merge.rmdup.stat:24901367 + 0 mapped (93.46% : N/A)
RNAPII_8WG16.merge.rmdup.stat:23397509 + 0 mapped (94.84% : N/A)
RNAPII_S2P.merge.rmdup.stat:26655659 + 0 mapped (95.36% : N/A)
RNAPII_S5P.merge.rmdup.stat:13680963 + 0 mapped (90.78% : N/A)
RNAPII_S5PRepeat.merge.rmdup.stat:3997567 + 0 mapped (82.22% : N/A)
RNAPII_S7P.merge.rmdup.stat:9759486 + 0 mapped (77.96% : N/A)
RNAPII_S7P.merge.rmdup.stat:9759486 + 0 primary mapped (77.96% : N/A)
去除PCR重復(fù)和不去除PCR重復(fù)的樣本都找一次peaks看一下眠副。
我們拿到兩批文件:合并的bam文件竣稽,去除PCR重復(fù)的合并bam文件。
下一步就是使用macs2尋找peaks了娃弓!
我們下一篇再見(jiàn)岛宦!