获取启动子区域之后，我们可能想要提取这些区域的碱基序列，有两种常见做法：
- 直接使用Bioconductor发行的R包
- 将区域存储为类似BED格式的文件，使用
BEDTools命令行工具
这里着重介绍第1种方案，采用一个关键的R包BSgenome (BS代表biostrings)。此包与之前介绍的GenomicFeatures类似，预先存储了不同特种，不同版本的基因组序列信息（部分依赖数据如图1，所有依赖数据见官方说明），如果没有你感兴趣序列的话可以考虑BEDTools工具。
图1 BSgenome序列数据
通过以下命令安装BSgenome：
> BiocManager::install("BSgenome")
导入小鼠的参考基因组序列：
> library(BSgenome.Mmusculus.UCSC.mm10)
> mm_gm <- BSgenome.Mmusculus.UCSC.mm10
序列基本信息查看
通过metadata命令可以查看序列的物种，版本，来源等信息：
> metadata(mm_gm)
$organism
[1] "Mus musculus"
$common_name
[1] "Mouse"
$genome
[1] "mm10"
$provider
[1] "UCSC"
$release_date
[1] "Dec. 2011"
$source_url
[1] "http://hgdownload.cse.ucsc.edu/goldenPath/mm10/chromosomes/"
通过seqinfo命令查看序列信息：
> seqinfo(mm_gm)
Seqinfo object with 66 sequences (1 circular) from mm10 genome:
seqnames seqlengths isCircular genome
chr1 195471971 FALSE mm10
chr2 182113224 FALSE mm10
chr3 160039680 FALSE mm10
chr4 156508116 FALSE mm10
chr5 151834684 FALSE mm10
... ... ... ...
chrUn_GL456392 23629 FALSE mm10
chrUn_GL456393 55711 FALSE mm10
chrUn_GL456394 24323 FALSE mm10
chrUn_GL456396 21240 FALSE mm10
chrUn_JH584304 114452 FALSE mm10
直接查看某条染色体序列：
> mm_gm$chrM
16299-letter DNAString object
seq: GTTAATGTAGCTTAATAACAAAGCAAAGCACTGAAA...TCTAATCATACTCTATTACGCAATAAACATTAACAA
定位碱基序列
序列本质上字符串，那么我们就可以使用已有的一段序列来搜索其出现的位置(使用Biostrings::matchPattern函数)。比如说我们在1号染色体上搜索“TCGATCGA”序列：
> matchPattern("TCGATCGA", mm_gm$chr1)
Views on a 195471971-letter DNAString subject
subject: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
views:
start end width
[1] 5118747 5118754 8 [TCGATCGA]
[2] 12846411 12846418 8 [TCGATCGA]
[3] 20403153 20403160 8 [TCGATCGA]
[4] 24329147 24329154 8 [TCGATCGA]
[5] 28627400 28627407 8 [TCGATCGA]
... ... ... ... ...
[73] 181302459 181302466 8 [TCGATCGA]
[74] 184731611 184731618 8 [TCGATCGA]
[75] 184836336 184836343 8 [TCGATCGA]
[76] 185637438 185637445 8 [TCGATCGA]
[77] 189056519 189056526 8 [TCGATCGA]
注意：这种方式只限于小规模地查找，不可进行大规模的序列比对。
查看区域序列
上一节我们获取了启动子区域信息(命令如下)：
数据下载地址
> chr1_gtf <- import("Mus_musculus.GRCm38.75_chr1.gtf.gz")
> chr1_pcg <- mm_gtf[mm_gtf$type == "gene" & mm_gtf$gene_biotype == "protein_coding"]
> chr1_pcg_3kb_up <- flank(chr1_pcg, width = 3000)
在我们提取启动子区域序列之前，需要关注一个细节，启动子区域的染色体编号和参考基因组染色体编号方式并不一致（seqlevels提取染色体编号）：
> all(seqlevels(chr1_pcg_3kb_up) %in% seqlevels(mm_gm))
[1] FALSE
这是因为我们前面使用的注释数据来自NCBI，其采用纯数字来编号染色体（如“1”，“2”），而BSgenome采用来自UCSC的基因组，采用的染色体编号方式为“chr1”, "chr2"等，通过seqlevelsStyle函数确认：
> seqlevelsStyle(chr1_pcg_3kb_up)
[1] "NCBI" "Ensembl" "MSU6" "AGPvF"
> seqlevelsStyle(mm_gm)
[1] "UCSC"
那么，这里就需要先统一染色体命名方式，这里将NCBI的序列转变为UCSC的风格：
> seqlevelsStyle(chr1_pcg_3kb_up) <- "UCSC"
> all(seqlevels(chr1_pcg_3kb_up) %in% seqlevels(mm_gm))
[1] TRUE
接下来就可以进行启动子区域的序列提取了，采用getSeq函数：
> promoters_seq <- getSeq(mm_gm, chr1_pcg_3kb_up)
> promoters_seq
DNAStringSet object of length 1240:
width seq
[1] 3000 ATTCTGAGATGTGGTTACTAGATCAATGGGAT...CGGCTAGCCGGGCCCAGCGCCCAGCCCCGCGG
[2] 3000 GAAGTGGTATATCTGCCTAGTCTAGGTGTGCA...GCTGTACTTAATCTGTGAGCACACATGCTAGT
[3] 3000 CTTAAAAACCTAGATATTCTATTTTTTTTTTT...CTTTGATAACGTCGTGAGCTCGGCTTCCAACA
[4] 3000 GAATTGGCACAGTTTCACATGATTGGTCCATT...GTACGGCCGCTGCAGCGCGACAGGGGCCGGGC
[5] 3000 AAATATAAAGTTAACATACAAAAACTAGTCGC...TCGGGGCGCGAGCTCGGGGCCGAACGCGAGGA
... ... ...
[1236] 3000 CAACATGGGTAGTAGTGGGGGAGCTTTAGTTC...GAGGGGCTGGCCTCACCAAGACGCAACAGGGA
[1237] 3000 AGGTGTGTTATATAATAATTGGTTTGACACTG...CTTAAAACTTGCTCTCTGGCTTCCTGGCGCCC
[1238] 3000 TTGGCCAGGTGATTGATCTTGTCCAACTGGAA...GTAAGGCCGGGCTATATGCAAACCGAGTTCCC
[1239] 3000 GGCATTCCCCTATACTGGGGCATAGAACCTTC...ATTTAAGGGTCTGCTCCCCACTGCTTACAGCC
[1240] 3000 GTAAATTTTCAGGTATATTTCTTTCTACTCTT...CTTTGATATTTCTGTGGTCCTTATTTCTAGGT
getSeq函数的两个参数分别为存储基因组序列的BSgenome对象和存储范围的GRanges对象。
最后，我们可以将提取的碱基序列以fasta格式存储，采用writeXStringSet命令：
> writeXStringSet(promoters_seq, filepath= "Mmusculus.UCSC.mm10.promoters.fasta", format = "fasta")
网友评论