49.《Bioinformatics Data Skills》之

上文介绍的GenomicFeatures包在统一方便的同时，也牺牲掉了灵活性。这里介绍另外一个R包rtracklayer，提供GTF/GFF，BED，BED Graph文件的读取，处理与导出。

本节数据下载地址： https://github.com/vsbuffalo/bds-files/tree/master/chapter-09-working-with-range-data

数据读取

读取采用import函数，可将上述文件类型读取为GenomicRanges对象，而且自动地分配meta-colmns列：

> library(rtracklayer)
> mm_gtf <- import("Mus_musculus.GRCm38.75_chr1.gtf.gz")
> names(mcols(mm_gtf))
 [1] "source"            "type"              "score"
 [4] "phase"             "gene_id"           "gene_name"
 [7] "gene_source"       "gene_biotype"      "transcript_id"
[10] "transcript_name"   "transcript_source" "tag"
[13] "exon_number"       "exon_id"           "ccds_id"
[16] "protein_id"

数据导出

采用export函数导出，同样可以导出为上述的文件类型。例如提取5个pseudogene并导出为gtf文件，可以：

> set.seed(0)
> pseudogene_i <- which(mm_gtf$gene_biotype == "pseudogene")
> sample_i <- sample(pseudogene_i, 5)
> pseudo_gene_sample <- mm_gtf[sample_i, ]
> export(pseudo_gene_sample, con = "pseudo_gene_sample.gtf", format = "GTF")

假如我们并不关心meta-column列的数据，可以将它们删除，保存为BED格式（BED文件仅包含染色体编号，起始与结束位点三列数据）：

> pseudo_gene_bed <- pseudo_gene_sample
> mcols(pseudo_gene_bed) <- NULL
> pseudo_gene_bed
GRanges object with 5 ranges and 0 metadata columns:
      seqnames              ranges strand
         <Rle>           <IRanges>  <Rle>
  [1]        1   88153788-88154571      +
  [2]        1   55443733-55445064      -
  [3]        1 121969217-121970524      -
  [4]        1   85332716-85342086      -
  [5]        1 119787769-119788086      +
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
> export(pseudo_gene_bed, con = "pseudo_gene_bed.bed", format="BED")

最后补充一下，rtracklayer提供和UCSC浏览器的交互，可以将GenomicRanges对象生成为UCSC的track并上传到UCSC网页，需要的话可以阅读一下官方文档了解一下。