117-Bioconductor没想象的那么简单-part7

刘小泽写于19.5.20 最近打算从头看一下Bioconductor这本书,目前最新版是2018的,https://bioconductor.github.io/BiocWorkshops/ ,感兴趣的小伙伴可以一起

先写写R的常用基础

  • read.csv函数的colClasses参数,可以在读入数据时就设置好各列的数据类型

    # 比如这样是读进来fname这个4列的数据框
    pdata <- read.csv( fname,
    colClasses = c("character", "factor", "integer", "factor") )
    
  • %in% 判断左边的数据是否在右边出现,返回逻辑值

  • 取子集:subset函数,利用逻辑值对行进行筛选,比如:

    subset(pdata, mol.biol %in% c("BCR/ABL", "NEG"))
    # pdata有一列叫mol.biol,这里选出mol.biol为BCR/ABL或NEG的行
    
  • formula的使用:基本格式就是y ~ x,左边是因变量(作图的纵坐标),右边是自变量(作图的横坐标)。比如有个数据x,要画mol.biol与age的关系,就是这样:

    boxplot(age ~ mol.biol, bcrabl)
    
  • R目前包含了超过15000个包;Bioconductor包含了超过1500个关于高通量数据分析与理解的包(RNA-seq,ChIP-seq,CNV,microarray等),包括软件包、注释包、实验包、流程包

  • 要查找某个R包的名称,支持部分匹配

    BiocManager::available("TxDb.Hsapiens")
    #> [1] "TxDb.Hsapiens.BioMart.igis"
    #> [2] "TxDb.Hsapiens.UCSC.hg18.knownGene"
    #> [3] "TxDb.Hsapiens.UCSC.hg19.knownGene"
    #> [4] "TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts" #> [5] "TxDb.Hsapiens.UCSC.hg38.knownGene"
    
  • 要看帮助文档:

    browseVignettes("simpleSingleCell")
    

开始练习

首先是准备工作

下载地址:https://github.com/Bioconductor/BiocWorkshops/blob/master/100_Morgan_RBiocForAll/CpGislands.Hsapiens.hg38.UCSC.bed

主要利用两个包:rtracklayerGenomicRanges ,其中rtracklayer 提供import函数,可以读取多种基因组文件作为对象(比如:BED、GTF、VCF、FASTA) ;GenomicRanges 又可以处理基因区域(比如:外显子、基因、ChIP peaks、called variants),这些信息都以坐标形式存储

library("rtracklayer") 
library("GenomicRanges")
导入数据

读入的数据是UCSC的BED文件,其中存储了人类基因组所有的CpG岛位置信息

文本大概长这样:包括了染色体、起始终止信息、每个CpG岛的ID

chr1    155188536   155192004   CpG:_361
chr1    2226773 2229734 CpG:_366
chr1    36306229    36307408    CpG:_110
chr1    47708822    47710847    CpG:_164
chr1    53737729    53739637    CpG:_221
chr1    144179071   144179313   CpG:_20

然后用函数读取(读取后看看和原来哪里不同)

# 选择文件
fname <- file.choose()
fname
file.exists(fname)

# 读入数据
cpg <- import(fname)
cpg

> head(cpg)
GRanges object with 6 ranges and 1 metadata column:
      seqnames              ranges strand |        name
         <Rle>           <IRanges>  <Rle> | <character>
  [1]     chr1 155188537-155192004      * |    CpG:_361
  [2]     chr1     2226774-2229734      * |    CpG:_366
  [3]     chr1   36306230-36307408      * |    CpG:_110
  [4]     chr1   47708823-47710847      * |    CpG:_164
  [5]     chr1   53737730-53739637      * |    CpG:_221
  [6]     chr1 144179072-144179313      * |     CpG:_20
  -------
  seqinfo: 254 sequences from an unspecified genome; no seqlengths

发现坐标出现了一点变化,比如第一行,源文件是155188536 155192004,现在是155188537-155192004

这是因为:原来的BED文件是0-based,而且是半开放区间(也就是包含左侧,不包含右侧坐标);Bioconductor转换得到的是1-based,而且是闭区间(左右两侧坐标都包含在内)

操作Genomic ranges

从上面的cpg来看,染色体名是存在像chr22_KI270734v1_random这样的形式,但是我们只想保留标准的常染色体和性染色体,因此先过滤一下

cpg <- keepStandardChromosomes(cpg, pruning.mode = "coarse") 
# 检查一下过滤结果
> levels(seqnames(cpg))
 [1] "chr1"  "chr2"  "chr3"  "chr4" 
 [5] "chr5"  "chr6"  "chr7"  "chr8" 
 [9] "chr9"  "chrX"  "chrY"  "chr10"
[13] "chr11" "chr12" "chr13" "chr14"
[17] "chr15" "chr16" "chr17" "chr18"
[21] "chr19" "chr20" "chr21" "chr22"

标准的GRange对象取其中的某一列使用$ ,比如:

head( cpg$name )
[1] "CpG:_361" "CpG:_366" "CpG:_110"
[4] "CpG:_164" "CpG:_221" "CpG:_20"

然后可以利用width()函数计算每个CpG岛的长度,因为这个值之间差别还较大,因此为了可视化可以利用log10处理一下:

hist(log10(width(cpg)))

取个子集试试,比如取出1、2染色体的cpg位点

subset(cpg, seqnames %in% c("chr1", "chr2"))
基因组注释

之前说到Bioconductor有一类叫注释包,举个例子,就是TxDb.*家族的包,这些包里一般会包含外显子、基因、转录本等,比如:TxDb.Hsapiens.UCSC.hg38.knownGene就是UCSC 存储的人类已知基因的注释

library("TxDb.Hsapiens.UCSC.hg38.knownGene")
# 比如提取所有转录本的坐标
tx <- transcripts(TxDb.Hsapiens.UCSC.hg38.knownGene)
tx
# 然后同样过滤不常用的染色体名称
tx <- keepStandardChromosomes(tx, pruning.mode="coarse")
tx
找overlap

这是一个较为常用的操作,就是在两个不同的Grange对象中找重叠区域,比如现在要找每个转录本可以与多少个CpG有交叉,那么query就是转录本对象txsubject就是CpG对象cpg

olaps <- countOverlaps(tx, cpg)
length(olaps) #  182435,表示一共有 182435个转录本
# 然后用table看看各个交叉范围(上面一行记录是有交叉的CpG的数量;下面一行是对应转录本的数量)
table(olaps) 
#> olaps

#> 0 1 2 3 4 5 6 7 8 9 10 11 
#> 94621 70551 10983 3228 1317 595 351 214 153 93 64 39 

#> 12 13 14 15 16 17 18 19 20 21 22 23 
#> 41 31 21 20 17 8 14 8 7 3 6 6 

#> 24 25 26 27 28 29 30 31 32 33 34 35 
#> 3 2 4 1 3 2 6 5 3 3 1 2 

#> 36 37 38 63
#> 3 1 1 4

可以把这个结果保存到转录本tx对象中

tx$cpgOverlaps <- countOverlaps(tx, cpg)
# 然后最后就会增加一列cpgOverlaps

那么有了这个结果,就很容易过滤了,比如想要看与CpG位点有10个以上交叉的转录本坐标

subset(tx, cpgOverlaps > 10)
Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related