从BAM文件中提取指定区域

发表于 | 分类于 | | 热度 °C

首先，看一下BAM文件的内容，BAM文件名为heart.bam

1
2
3
4
5
6
7
8
9
10
11
> samtools view heart.bam | head
HWI-ST1113:280:H7KE6ADXX:2:1201:2203:75105      16      chr19   3079081 35      36M     *       0       0       AGTCCATTTATTCTTTGTGTAGAGAAAGACTACTGA     JIIJJJJJJJJJIIGJJIJJJIJGHGHHFFFDFCC@     AS:i:0  XS:i:-12        XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:2206:14516:57404     0       chr19   3079310 31      36M     *       0       0       CCACAAGATGACTCCAGACCTGAGTGATGCAAACAG     @@@FFFFFGHFHHJJGIJJEHJJJEEHIJJJJJJJI     AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:2214:19987:44958     16      chr19   3080044 30      36M     *       0       0       ACCAGAGAAACAGACAGCTTCTGGGACAGGTGGAAG     ####AC>C?<+,BC?7=7CBABBB@7BA7AA+A;==     AS:i:0  XS:i:-4 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:2:2214:13543:99477     16      chr19   3082639 31      36M     *       0       0       GCTCTAGAAAAAATGGAAGCAAATTCAGCCAAGAGG     IJIJJJIJJJJJJIJIJJJJJJJHHHHHFFFFFCCC     AS:i:0  XS:i:-5 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:2:1115:6716:37571      0       chr19   3083875 31      36M     *       0       0       TACCCCAACTCGGTGCTGGTTTGGGTTTCTTTTTGG     CCCFFFFFHHHHHFGHJJJFHIGJJEHIIJJJJJJJ     AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:1216:1670:96735      0       chr19   3085742 31      36M     *       0       0       TCACTACTGCCTGAGACTTCGCCTACTCATCATTGT     CCCFFFFFHHHHHIIJJIJJJJIIJJJIJJJIJJIG     AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:1202:4496:44872      16      chr19   3087466 42      36M     *       0       0       TCCCTTTGGGGATCTTTCATTAGAAGCAGTTATAGT     JJJJIJJJJJJIJHJJJJJJJJJHHHHHFFFFFCCC     AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:1111:11809:31419     16      chr19   3088647 31      36M     *       0       0       ATTCATATGGAATAACAAAAAACTGAGGATAGCAAA     JIIHJJIIGIJJIJJIJJIIIEIHHGHHFFFFD@CC     AS:i:0  XS:i:-5 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:2:1113:1925:16628      16      chr19   3090618 42      36M     *       0       0       TATTGTCAAGTATATTTTAGGACAATTTTCTTGATT     JIIIIGJJIIIGIIIJIIGHEEJGGHGHFFFFF@CC     AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:2205:2262:65371      0       chr19   3090666 35      36M     *       0       0       CAATCATTTTGATAACAGGAATCAGGAGAACTAACT     CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJJJJ     AS:i:0  XS:i:-11        XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU

下面，准备一个区域信息文件，即BED文件，命名为chr19.bed，查看一下内容

1
2
> less chr19.bed
chr19   3079000 3080000

使用bedtools intersect取交集，并查看生成的target.bam文件

1
2
3
4
bedtools intersect -a heart.bam -b chr19.bed > target.bam
> samtools view target.bam | head
HWI-ST1113:280:H7KE6ADXX:2:1201:2203:75105      16      chr19   3079081 35      36M     *       0       0       AGTCCATTTATTCTTTGTGTAGAGAAAGACTACTGA    JIIJJJJJJJJJIIGJJIJJJIJGHGHHFFFDFCC@    AS:i:0  XS:i:-12        XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:2206:14516:57404     0       chr19   3079310 31      36M     *       0       0       CCACAAGATGACTCCAGACCTGAGTGATGCAAACAG    @@@FFFFFGHFHHJJGIJJEHJJJEEHIJJJJJJJI    AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU

在目标范围chr19:3079000-3080000之间，共有两个reads，即3079081-3079116和3079310-3079345，读长皆为36 bp，这是最简单的情形，容易理解，因为两段读长都完全位于目标范围之内。

假设，读长只有一部分位于目标区域之内，结果如何？修改chr19.bed为chr19b.bed，并更改目标范围，看看结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
chr19   3079081 3079310 # [1]第一段完全覆盖，第二段只覆盖第一个碱基

HWI-ST1113:280:H7KE6ADXX:2:1201:2203:75105      16      chr19   3079081 35      36M     *       0       0       AGTCCATTTATTCTTTGTGTAGAGAAAGACTACTGA     JIIJJJJJJJJJIIGJJIJJJIJGHGHHFFFDFCC@     AS:i:0  XS:i:-12        XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:2206:14516:57404     0       chr19   3079310 31      36M     *       0       0       CCACAAGATGACTCCAGACCTGAGTGATGCAAACAG     @@@FFFFFGHFHHJJGIJJEHJJJEEHIJJJJJJJI     AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU

chr19   3079081 3079345 # [2]两段完全覆盖

HWI-ST1113:280:H7KE6ADXX:2:1201:2203:75105      16      chr19   3079081 35      36M     *       0       0       AGTCCATTTATTCTTTGTGTAGAGAAAGACTACTGA     JIIJJJJJJJJJIIGJJIJJJIJGHGHHFFFDFCC@     AS:i:0  XS:i:-12        XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:2206:14516:57404     0       chr19   3079310 31      36M     *       0       0       CCACAAGATGACTCCAGACCTGAGTGATGCAAACAG     @@@FFFFFGHFHHJJGIJJEHJJJEEHIJJJJJJJI     AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU

chr19   3079116 3079310 # [3]第一段只覆盖最后一个碱基，第二段只覆盖第一个碱基

HWI-ST1113:280:H7KE6ADXX:1:2206:14516:57404     0       chr19   3079310 31      36M     *       0       0       CCACAAGATGACTCCAGACCTGAGTGATGCAAACAG     @@@FFFFFGHFHHJJGIJJEHJJJEEHIJJJJJJJI     AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU

chr19   3079116 3079345 # [4]第一段只覆盖最后一个碱基，第二段完全覆盖

HWI-ST1113:280:H7KE6ADXX:1:2206:14516:57404     0       chr19   3079310 31      36M     *       0       0       CCACAAGATGACTCCAGACCTGAGTGATGCAAACAG     @@@FFFFFGHFHHJJGIJJEHJJJEEHIJJJJJJJI     AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU

chr19   3079081 3079309 # [5]第一段完全覆盖，第二段无覆盖

HWI-ST1113:280:H7KE6ADXX:2:1201:2203:75105      16      chr19   3079081 35      36M     *       0       0       AGTCCATTTATTCTTTGTGTAGAGAAAGACTACTGA     JIIJJJJJJJJJIIGJJIJJJIJGHGHHFFFDFCC@     AS:i:0  XS:i:-12        XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU

chr19   3079116 3079344 # [6]第一段只覆盖最后一个碱基，第二段覆盖绝大部分

HWI-ST1113:280:H7KE6ADXX:1:2206:14516:57404     0       chr19   3079310 31      36M     *       0       0       CCACAAGATGACTCCAGACCTGAGTGATGCAAACAG     @@@FFFFFGHFHHJJGIJJEHJJJEEHIJJJJJJJI     AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU

chr19   3079090 3079330 # [7]两段部分覆盖

HWI-ST1113:280:H7KE6ADXX:2:1201:2203:75105      16      chr19   3079081 35      36M     *       0       0       AGTCCATTTATTCTTTGTGTAGAGAAAGACTACTGA     JIIJJJJJJJJJIIGJJIJJJIJGHGHHFFFDFCC@     AS:i:0  XS:i:-12        XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU
HWI-ST1113:280:H7KE6ADXX:1:2206:14516:57404     0       chr19   3079310 31      36M     *       0       0       CCACAAGATGACTCCAGACCTGAGTGATGCAAACAG     @@@FFFFFGHFHHJJGIJJEHJJJEEHIJJJJJJJI     AS:i:0  XS:i:-6 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:36 YT:Z:UU

综上可见，对于一个read而言，只要其中一个碱基（最后一个例外）落入目的范围，则会被认为是与此目的区域有交集并保留。

  • 本文作者：括囊无誉
  • 本文链接： CHIPSEQ/intersectbam/
  • 版权声明： 本博客所有文章均为原创作品，转载请注明出处！
------ 本文结束 ------
坚持原创文章分享，您的支持将鼓励我继续创作！