寻找区域间的overlap算得上是基因组最常见的操作，貌似简单的寻找overlap追究细节之后也会变得非常复杂，比如为了解决RNA-seq表达定量的问题就衍生出一系列专业的工具或者R包（例如RSEM，TopHat，HTSeq，GenomicAlignments等等）。
今天我们通过一个简单的实例来学习寻找overlaps的操作：统计落在外显子区域的dbSNP数量。
本节数据下载
数据准备
首先，导入小鼠1号染色体上的dbSNP文件（dbSNP包括了单核苷酸多态以及碱基的插入缺失，短的串连重复，多核苷酸多态）：
> library(rtracklayer)
> dbsnp <- import("mm10_snp137_chr1_trunc.bed.gz")
首先查看长度的分布（永远对自己的数据保持怀疑）：
> summary(width(dbsnp))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 1.000 1.000 1.142 1.000 732.000
最长的dbSNP居然长达732kb，貌似不太合理，它的名字为：
> dbsnp$name[which.max(width(dbsnp))]
[1] "rs232497063"
不过我们通过UCSC浏览器可以确定它是一个真实的的dbsnp（图1）：
图1
此外，大量dbSNP长度为0，这是因为相对于参考基因组来说它们是插入序列。不过其长度为0的话，我们无法确定其和外显子区域的overlap，为了方便统计，使用resize函数将它们的size定义为1：
> zw_i <- which(width(dbsnp) == 0)
> dbsnp_resize <- dbsnp
> dbsnp_resize[zw_i] <- resize(dbsnp[zw_i], width = 1)
最后，我们导入小鼠1号染色体上所有外显子区域（忽略链位置信息）：
> library(TxDb.Mmusculus.UCSC.mm10.ensGene)
> txdb <- TxDb.Mmusculus.UCSC.mm10.ensGene
> collapsed_exons <- reduce(exons(txdb), ignore.strand = T)
> chr1_collapsed_exons <- collapsed_exons[seqnames(collapsed_exons) == "chr1"]
寻找overlaps
采用findOverlaps函数确定dbSNP与外显子区域的交集：
> hits <- findOverlaps(dbsnp_resize, chr1_collapsed_exons, ignore.strand = T)
外显子区域dbSNP数目与其占所有dbSNP的比例为：
> length(unique(queryHits(hits)))
[1] 57623
> length(unique(queryHits(hits)))/length(dbsnp_resize)
[1] 0.02134185
可见只有约2%的dbSNP落在外显子区域。
我们还可以筛选出落在外显子上面的dbSNP，一种方式为根据hits检索，另外一种方式为采用函数subsetByOverlaps:
> subsetByOverlaps(dbsnp_resize, chr1_collapsed_exons, ignore.strand = T)
GRanges object with 57623 ranges and 2 metadata columns:
seqnames ranges strand | name score
<Rle> <IRanges> <Rle> | <character> <numeric>
[1] chr1 43032144 + | rs250123171 0
[2] chr1 36713805 + | rs50487270 0
[3] chr1 132567494 + | rs247294715 0
[4] chr1 160995431 + | rs47617081 0
[5] chr1 84036552-84036553 + | rs216202117 0
... ... ... ... . ... ...
[57619] chr1 188263219 + | rs13476293 0
[57620] chr1 134780954 + | rs218301913 0
[57621] chr1 130270464 + | rs266050681 0
[57622] chr1 107380295 + | rs224267626 0
[57623] chr1 98421207 + | rs224196900 0
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
最后，使用countOverlaps函数统计各个外显子区域的dbSNP数目：
var_count <- countOverlaps(chr1_collapsed_exons, dbsnp_resize, ignore.strand = T)
为了方便追踪，我们将这个数目作为外显子区域的meta-data：
> chr1_collapsed_exons$var_count <- var_count
> chr1_collapsed_exons
GRanges object with 15048 ranges and 1 metadata column:
seqnames ranges strand | var_count
<Rle> <IRanges> <Rle> | <integer>
[1] chr1 3054233-3054733 * | 1
[2] chr1 3102016-3102125 * | 0
[3] chr1 3205901-3207317 * | 17
[4] chr1 3213439-3216968 * | 21
[5] chr1 3421702-3421901 * | 1
... ... ... ... . ...
[15044] chr1 195169702-195169801 * | 0
[15045] chr1 195170991-195171168 * | 0
[15046] chr1 195176553-195176715 * | 1
[15047] chr1 195228278-195228398 * | 0
[15048] chr1 195240910-195241007 * | 0
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
第一个区域唯一的dbSNP叫做rs233033126：
> dbsnp_resize$name[queryHits(hits)[which(subjectHits(hits) == 1)]]
[1] "rs233033126"
根据UCSC浏览器，该dbSNP确实在区域chr1:3054233-3054733内出现（图2），不过该区域包含的dbSNP明显不止1个，说明我们的数据可能存在问题。
图2
网友评论