獲取啟動子區(qū)域之后，我們可能想要提取這些區(qū)域的堿基序列，有兩種常見做法：

1. 直接使用Bioconductor發(fā)行的R包
2. 將區(qū)域存儲為類似BED格式的文件，使用BEDTools命令行工具

這里著重介紹第1種方案，采用一個關(guān)鍵的R包BSgenome (BS代表biostrings)。此包與之前介紹的GenomicFeatures類似，預(yù)先存儲了不同特種，不同版本的基因組序列信息（部分依賴數(shù)據(jù)如圖1，所有依賴數(shù)據(jù)見官方說明），如果沒有你感興趣序列的話可以考慮BEDTools工具。

圖1 BSgenome序列數(shù)據(jù)

通過以下命令安裝BSgenome：

> BiocManager::install("BSgenome")

導(dǎo)入小鼠的參考基因組序列：

> 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

定位堿基序列

序列本質(zhì)上字符串，那么我們就可以使用已有的一段序列來搜索其出現(xiàn)的位置(使用Biostrings::matchPattern函數(shù))。比如說我們在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]

注意：這種方式只限于小規(guī)模地查找，不可進(jìn)行大規(guī)模的序列比對。

查看區(qū)域序列

上一節(jié)我們獲取了啟動子區(qū)域信息(命令如下)：

數(shù)據(jù)下載地址

> 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)

在我們提取啟動子區(qū)域序列之前，需要關(guān)注一個細(xì)節(jié)，啟動子區(qū)域的染色體編號和參考基因組染色體編號方式并不一致（seqlevels提取染色體編號）：

> all(seqlevels(chr1_pcg_3kb_up) %in% seqlevels(mm_gm))
[1] FALSE

這是因?yàn)槲覀兦懊媸褂玫淖⑨寯?shù)據(jù)來自NCBI，其采用純數(shù)字來編號染色體（如“1”，“2”），而BSgenome采用來自UCSC的基因組，采用的染色體編號方式為“chr1”, "chr2"等，通過seqlevelsStyle函數(shù)確認(rèn)：

> seqlevelsStyle(chr1_pcg_3kb_up)
[1] "NCBI"    "Ensembl" "MSU6"    "AGPvF"
> seqlevelsStyle(mm_gm)
[1] "UCSC"

那么，這里就需要先統(tǒng)一染色體命名方式，這里將NCBI的序列轉(zhuǎn)變?yōu)閁CSC的風(fēng)格：

> seqlevelsStyle(chr1_pcg_3kb_up) <- "UCSC"
> all(seqlevels(chr1_pcg_3kb_up) %in% seqlevels(mm_gm))
[1] TRUE

接下來就可以進(jìn)行啟動子區(qū)域的序列提取了，采用getSeq函數(shù)：

> 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函數(shù)的兩個參數(shù)分別為存儲基因組序列的BSgenome對象和存儲范圍的GRanges對象。

最后，我們可以將提取的堿基序列以fasta格式存儲，采用writeXStringSet命令：

> writeXStringSet(promoters_seq, filepath= "Mmusculus.UCSC.mm10.promoters.fasta", format = "fasta")