今天我们通过2个实战来掌握函数gaps，setdiff与reduce在GenomicRanges中的使用：
- 获取基因间区域；
- 获取基因的内含子区域。
Gaps函数了解
首先了解一下提取基因间区域需要用到的gaps函数，此函数用于返回范围间的空白区域，之前在IRanges包中介绍过。但GRanges对象增加了strand信息，结果有所不同，通过一个简单的例子说明：
首先创建一个GRanges对象gr2：
> library(GenomicRanges)
> gr2 <- GRanges(c("chr1", "chr2"), IRanges(start = c(5, 10), width = 6), strand = c("+", "-"), seqlengths = c(chr1 = 20, chr2 = 40))
> gr2
GRanges object with 2 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 5-10 +
[2] chr2 10-15 -
-------
seqinfo: 2 sequences from an unspecified genome
对gr2使用gaps函数：
> gaps(gr2)
GRanges object with 8 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 1-4 +
[2] chr1 11-20 +
[3] chr1 1-20 -
[4] chr1 1-20 *
[5] chr2 1-40 +
[6] chr2 1-9 -
[7] chr2 16-40 -
[8] chr2 1-40 *
-------
seqinfo: 2 sequences from an unspecified genome
观察以上结果，发现gaps并没有如我们想象那样返回一条区域的空白，而是根据链信息保留所有的可能性（8条）。这是严谨正确的做法，然而有时我们不考虑链信息，可以将GRanges对象的链设置为“*”，通过如下方式获取gaps区域：
> gr3 <- gr2
> gaps(gr3)[strand(gaps(gr3)) == "*"]
GRanges object with 4 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 1-4 *
[2] chr1 11-20 *
[3] chr2 1-9 *
[4] chr2 16-40 *
-------
seqinfo: 2 sequences from an unspecified genome
提取基因间区域
以小鼠的所有基因转录本数据为例：
> library(TxDb.Mmusculus.UCSC.mm10.ensGene)
> txdb <- TxDb.Mmusculus.UCSC.mm10.ensGene
通过seqinfo函数获取所有染色体区域：
> chr_rngs <- GRanges(seqinfo(txdb))
> head(chr_rngs, 2)
GRanges object with 2 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
chr1 chr1 1-195471971 *
chr2 chr2 1-182113224 *
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
使用reduce函数融合所有重叠的转录本区域：
> collapsed_tx <- reduce(transcripts(txdb))
避免考虑链信息：
> strand(collapsed_tx) <- "*"
最终采用setdiff命令得到基因间区域：
> intergenic <- setdiff(chr_rngs, collapsed_tx)
注：这里转录本间的范围都当作基因间区域，相对只考虑蛋白编码基因来说范围更少。
提取内含子区域
这里介绍两种方式获取基因的内含子区域：
-
直接使用封装好的函数，非常简单，如下：
> mm_introns <- intronsByTranscript(txdb) 为了练习，采用手动方式提取，内含子区域通过转录本区域减去外显子区域获得。这部分的内容稍有难度，掌握之后可以加深我们对基因组操作的理解。简单起见，这里只提取淀粉酶基因（amylase 1，Ensembl id：ENSMUSG00000074264）的内含子区域，具体过程如下：
首先根据基因id获取amy1基因的转录本：
> amy1 <- transcriptsBy(txdb, "gene")$ENSMUSG00000074264
> amy1
GRanges object with 5 ranges and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
<Rle> <IRanges> <Rle> | <integer> <character>
[1] chr3 113555710-113577830 - | 18879 ENSMUST00000067980
[2] chr3 113555953-113574762 - | 18880 ENSMUST00000106540
[3] chr3 113556149-113562018 - | 18881 ENSMUST00000172885
[4] chr3 113562690-113574272 - | 18882 ENSMUST00000142505
[5] chr3 113564987-113606699 - | 18883 ENSMUST00000174147
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
提取小鼠所有转录本的外显子区域：
mm_exons <- exonsBy(txdb, 'tx')
> head(mm_exons,2)
GRangesList object of length 2:
$`1`
GRanges object with 1 range and 3 metadata columns:
seqnames ranges strand | exon_id exon_name exon_rank
<Rle> <IRanges> <Rle> | <integer> <character> <integer>
[1] chr1 3054233-3054733 + | 1 <NA> 1
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
$`2`
GRanges object with 1 range and 3 metadata columns:
seqnames ranges strand | exon_id exon_name exon_rank
<Rle> <IRanges> <Rle> | <integer> <character> <integer>
[1] chr1 3102016-3102125 + | 2 <NA> 1
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
由于外显子区域是以tx_id（转录本ID）为名字的GRangesList，这里我们将amy1的转录本也转换为以tx_id命名的GRangesList对象：
> amy1_tx <- split(amy1, amy1$tx_id)
> head(amy1_tx, 2)
GRangesList object of length 2:
$`18879`
GRanges object with 1 range and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
<Rle> <IRanges> <Rle> | <integer> <character>
[1] chr3 113555710-113577830 - | 18879 ENSMUST00000067980
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
$`18880`
GRanges object with 1 range and 2 metadata columns:
seqnames ranges strand | tx_id tx_name
<Rle> <IRanges> <Rle> | <integer> <character>
[1] chr3 113555953-113574762 - | 18880 ENSMUST00000106540
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
接下来，我们通过匹配Amy1基因的tx_id来获取Amy1的外显子区域：
> amy1_exons <- mm_exons[match(names(amy1_tx), names(mm_exons))]
在确定转录本ID与外显子区域完全匹配的情况下，使用psetdiff命令获取该基因内含子区域(注：psetdiff用于配对数据，效果类似setdiff)：
> all(names(amy1_tx) == names(amy1_exons))
[1] TRUE
> amy1_intron <- psetdiff(unlist(amy1_tx), amy1_exons)
通过与方法一获取的内含子区域对比，检验结果是否正确：
> identical(mm_introns[names(amy1_tx)], amy1_intron)
[1] TRUE
