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")
开始练习
首先是准备工作
主要利用两个包：rtracklayer和GenomicRanges ，其中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就是转录本对象tx，subject就是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)
