50.《Bioinformatics Data Skills》之

本节数据下载地址

上一节使用rtracklayer包加载了小鼠的1号染色体的基因注释数据（如下）。

> library(rtracklayer)
> mm_gtf <- import("Mus_musculus.GRCm38.75_chr1.gtf.gz")

此文件的type列和gene_biotype列内容分别如下：

> table(mm_gtf$type)
       gene  transcript        exon         CDS start_codon  stop_codon
       2027        4993       36128       25901        2290        2299
        UTR
       7588

> table(mm_gtf$gene_biotype)
             antisense                lincRNA                  miRNA
                   480                    551                    354
              misc_RNA polymorphic_pseudogene   processed_transcript
                    93                     61                    400
        protein_coding             pseudogene                   rRNA
                 77603                    978                     69
        sense_intronic      sense_overlapping                 snoRNA
                    21                      4                    297
                 snRNA
                   315

我们可以通过以下条件获取蛋白编码基因区域：

> chr1_pcg <- mm_gtf[mm_gtf$type == "gene" & mm_gtf$gene_biotype == "protein_coding"]

确认所选区域长度与数量是否正常：

> summary(width(chr1_pcg))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
     78    9429   25754   60640   62422 1075874
> length(chr1_pcg)
[1] 1240

启动子位于基因上游100到1kbp左右，这里定义其为蛋白编码区域上游3Kbp的区域。可以通过以下两种方式提取启动子区域：

1. 使用flank函数（详见IRanges操作）

> chr1_pcg_3kb_up <- flank(chr1_pcg, width = 3000)

注：使用help(flank)可以看到参数ignore.strand默认为FALSE，说明此函数还考虑了链信息。

2. 使用promoters函数
   此函数可以直接回溯启动子区域，默认为上游2Kbp与下游200bp作为启动子区域。设置相同的参数可以看出结果是一致的：

> chr1_pcg_3kb_up2 <- promoters(chr1_pcg, upstream=3000, downstream=0)
> identical(chr1_pcg_3kb_up, chr1_pcg_3kb_up2)
[1] TRUE

【完】