R语言教程之序列处理
Biostrings最后一节，介绍模式匹配和序列比对的相关函数和操作。
下面我们使用拟南芥基因转录起始点上游1kb的序列进行分析。序列文件可以从TAIR网站（http://www.arabidopsis.org）下载。先用readDNAStringSet函数从FASTA文件中读取序列并查看头2个序列的信息：
library(Biostrings) upstream.1k <- readDNAStringSet(file.choose(), "fasta")head(upstream.1k, 2)
## A DNAStringSet instance of length 2 ## width seq names ## [1] 1000 ACAAGCATGCTTAGCATACGT...ATCTTAATGTGGATAGTGCT AT1G08520 | chr1:... ## [2] 1000 CGTAACAGTAACAACTATATT...CTCCTGGAGAAGAGAAGACT AT1G08530 | chr1:...
# 获取第29条序列DNA.str <- upstream.1k[[29]]
一、模式匹配（pattern match）
“模式”是一段序列，可以是我们通常说的DNA motif，当然也可以是高通量测序获得的reads。模式匹配就是将这些序列比对或者map到目标序列如染色体序列上。模式匹配可用于统计和分析motif/reads在目标序列上的分布情况。
1、单模式匹配
在一般序列模式匹配的应用中，无论是查询模式还是目标序列都比较少，使用Biostrings的matchPattern和vmatchPattern函数完全可以胜任这方面的数据处理。这一系列的函数有四个，两个函数返回Views对象，另外两个函数统计匹配的数量：
matchPattern()：1个查询模式1条序列countPattern()：1个查询模式1条序列，仅计数vmatchPattern()：1个查询模式n条序列vcountPattern()：1个查询模式n条序列，仅计数
这四个函数的参数均为： pattern, 匹配模式 subject, 参考序列 max.mismatch=0, 最大错配数 min.mismatch=0, 最小错配数 with.indels=FALSE/TRUE, 是否允许Indels（插入/缺失）。如果设置为TRUE，min.mismatch必需为0，而max.mismatch解释为编辑位点与匹配区域间的距离 fixed=TRUE/FALSE, 是否允许按照IUPAC代码进行简并碱基的匹配 algorithm="auto", 可设"auto", "naive-exact", "naive-inexact", "boyer-moore", "shift-or" or "indels"等值。但不能随便设置，得看前面几个参数的设置情况。 下面我们看看ABRE motif（ACGTGKC）的分布情况：
## 1模式查询1条序列：非模糊匹配，不能使用IUPAC简并代码，所以结果找不到（views:# NONE）matchPattern("ACGTGKC", DNA.str)
## Views on a 1000-letter DNAString subject ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC ## views: NONE
# 模糊匹配，使用IUPAC简并代码，K可匹配G/TmatchPattern("ACGTGKC", DNA.str, fixed = FALSE)
## Views on a 1000-letter DNAString subject ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC ## views: ## start end width ## [1] 17 23 7 [ACGTGTC] ## [2] 65 71 7 [ACGTGTC]
# 仅计数countPattern("ACGTGKC", DNA.str, fixed = FALSE)
## [1] 2
fixed参数设置为FALSE时，如果目标序列上有非确定碱基如N，得到的结果会很让人失望：
# 1模式查询n条序列，仅计数：vc <- vcountPattern("ACGTGKC", upstream.1k, fixed = FALSE)# 看看最多匹配数，太多了max(vc)
## [1] 970
# 用subset函数提取具有最多匹配数的序列子集subs <- subset(upstream.1k, vc == max(vc))# 原来匹配的都是NNN...，怪不得那么多matchPattern("ACGTGKC", subs[[1]], fixed = FALSE)
## Views on a 1000-letter DNAString subject ## subject: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNGAATTCGTCGACCAGGACGGCGGA ## views: ## start end width ## [1] 1 7 7 [NNNNNNN] ## [2] 2 8 7 [NNNNNNN] ## [3] 3 9 7 [NNNNNNN] ## [4] 4 10 7 [NNNNNNN] ## [5] 5 11 7 [NNNNNNN] ## [6] 6 12 7 [NNNNNNN] ## [7] 7 13 7 [NNNNNNN] ## [8] 8 14 7 [NNNNNNN] ## [9] 9 15 7 [NNNNNNN] ## ... ... ... ... ... ## [962] 962 968 7 [NNNNNNN] ## [963] 963 969 7 [NNNNNNN] ## [964] 964 970 7 [NNNNNNN] ## [965] 965 971 7 [NNNNNNN] ## [966] 966 972 7 [NNNNNNN] ## [967] 967 973 7 [NNNNNNN] ## [968] 968 974 7 [NNNNNNN] ## [969] 969 975 7 [NNNNNNN] ## [970] 970 976 7 [NNNNNNN]
从上面的结果可以看出，如果fixed=FALSE，模式和目标序列的碱基都要应用IUPAC简并规则。但实际应用是我们往往要忽略目标序列上的不确定碱基，这时应设置fixed="subject"，固定目标序列上的碱基：
vc <- vcountPattern("ACGTGKC", upstream.1k, fixed = "subject")
subs <- subset(upstream.1k, vc == max(vc))# 查看结果，这才是我们想要的：matchPattern("ACGTGKC", subs[[1]], fixed = "subject")
## Views on a 1000-letter DNAString subject ## subject: TTACCCGAGTTGGGTAATCGACTCGGTTGGTT...CCTTTTCATTTTCTTCTCCCTCTCCTGGGTT ## views: ## start end width ## [1] 171 177 7 [ACGTGGC] ## [2] 350 356 7 [ACGTGTC] ## [3] 421 427 7 [ACGTGTC] ## [4] 449 455 7 [ACGTGGC]
上面语句如果往清楚里写应是这样的：
matchPattern("ACGTGKC", subs[[1]], fixed = c(subject = TRUE, pattern = FALSE))
## Views on a 1000-letter DNAString subject ## subject: TTACCCGAGTTGGGTAATCGACTCGGTTGGTT...CCTTTTCATTTTCTTCTCCCTCTCCTGGGTT ## views: ## start end width ## [1] 171 177 7 [ACGTGGC] ## [2] 350 356 7 [ACGTGTC] ## [3] 421 427 7 [ACGTGTC] ## [4] 449 455 7 [ACGTGGC]
模式匹配的结果可进一步使用，比如要获取翻译起始点上游有2个以上ABRE motif的序列名称，这些基因可能对ABA有较强的响应：
names(subset(upstream.1k, vc > 2))
## [1] "AT1G24580 | chr1:8709180-8710179 FORWARD LENGTH=1000" ## [2] "AT1G09210 | chr1:2976751-2977750 REVERSE LENGTH=1000" ## [3] "AT1G19490 | chr1:6754000-6754999 REVERSE LENGTH=1000" ## [4] "AT1G07870 | chr1:2432058-2433057 REVERSE LENGTH=1000" ## [5] "AT1G51140 | chr1:18945754-18946753 REVERSE LENGTH=1000" ## [6] "AT1G20180 | chr1:6995210-6996209 FORWARD LENGTH=1000" ## [7] "AT1G77450 | chr1:29098954-29099953 FORWARD LENGTH=1000" ## [8] "AT1G72510 | chr1:27302366-27303365 FORWARD LENGTH=1000" ## [9] "AT1G60190 | chr1:22197403-22198402 FORWARD LENGTH=1000" ## [10] "AT1G79040 | chr1:29735016-29736015 FORWARD LENGTH=1000" ## [11] "AT1G55520 | chr1:20728092-20729091 REVERSE LENGTH=1000" ## [12] "AT2G04350 | chr2:1514805-1515804 FORWARD LENGTH=1000" ## [13] "AT2G22240 | chr2:9454095-9455094 REVERSE LENGTH=1000" ## [14] "AT3G63350 | chr3:23398468-23399467 FORWARD LENGTH=1000" ## [15] "AT3G33154 | chr3:13964599-13965598 FORWARD LENGTH=1000" ## [16] "AT3G42065 | chr3:14254659-14255658 FORWARD LENGTH=1000" ## [17] "AT3G63070 | chr3:23301365-23302364 FORWARD LENGTH=1000" ## [18] "AT3G49440 | chr3:18336318-18337317 REVERSE LENGTH=1000" ## [19] "AT3G10985 | chr3:3441528-3442527 FORWARD LENGTH=1000" ## [20] "AT3G15280 | chr3:5144694-5145693 REVERSE LENGTH=1000" ## [21] "AT3G03680 | chr3:906500-907499 FORWARD LENGTH=1000" ## [22] "AT3G02140 | chr3:386539-387538 REVERSE LENGTH=1000" ## [23] "AT3G53040 | chr3:19666488-19667487 REVERSE LENGTH=1000" ## [24] "AT4G21940 | chr4:11639802-11640801 FORWARD LENGTH=1000" ## [25] "AT4G25570 | chr4:13055624-13056623 REVERSE LENGTH=1000" ## [26] "AT4G27410 | chr4:13709150-13710149 REVERSE LENGTH=1000" ## [27] "AT5G50360 | chr5:20506339-20507338 REVERSE LENGTH=1000" ## [28] "AT5G10490 | chr5:3304762-3305761 REVERSE LENGTH=1000" ## [29] "AT5G15950 | chr5:5204869-5205868 FORWARD LENGTH=1000" ## [30] "AT5G28667 | chr5:10689062-10690061 FORWARD LENGTH=1000" ## [31] "AT5G06980 | chr5:2166468-2167467 FORWARD LENGTH=1000" ## [32] "AT5G15948 | chr5:5204871-5205870 FORWARD LENGTH=1000"
掩膜序列生成函数maskMotif调用的是matchPattern函数，可以在参数中设置matchPattern的参数。下面两条语句的结果等价：
as(maskMotif(DNA.str, "ACGTGKC", fixed = "subject"), "Views")
## Views on a 1000-letter DNAString subject ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC ## views: ## start end width ## [1] 1 16 16 [AGAGAAGGCTCCAAAC] ## [2] 24 64 41 [TCAGCCATATTATATCATTCGTCCGAGAGAAGGCTCCACAC] ## [3] 72 1000 929 [TCAGCCATATGTCTATCCGACGTAA...CTTTGCTCATTATATGTCCTCAGC]
gaps(matchPattern("ACGTGKC", DNA.str, fixed = "subject"))
## Views on a 1000-letter DNAString subject ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC ## views: ## start end width ## [1] 1 16 16 [AGAGAAGGCTCCAAAC] ## [2] 24 64 41 [TCAGCCATATTATATCATTCGTCCGAGAGAAGGCTCCACAC] ## [3] 72 1000 929 [TCAGCCATATGTCTATCCGACGTAA...CTTTGCTCATTATATGTCCTCAGC]
2、多模式匹配与PDict类
如果使用matchPattern或vmatchPattern进行多个模式的序列查询，那就必需使用循环。但是在R语言中进行循环运算是非常耗时的，如果要把高通量测序获得的数以千万计的reads匹配（map）到基因或基因组上，运行时间将无法让人忍受。这样的工作最好让C语言程序来做。Biostrings中的很多工作事实上也是调用C语言程序来完成的。和前面一样，多模式匹配的函数也可以获得匹配信息或仅计数：
matchPDict()：n个查询模式1条序列countPDict()：n个查询模式1条序列，仅计数vmatchPDict()：n个查询模式n条序列vcountPDict()：n个查询模式n条序列，仅计数
这些函数的参数包括：pdict, subject, max.mismatch, min.mismatch, with.indels, fixed, algorithm和verbose
pdict是TB_PDict类对象，它通过继承虚拟类PDict而来：
getClass("TB_PDict")
## Class "TB_PDict" [package "Biostrings"] ## ## Slots: ## ## Name: threeparts dict0 constant_width dups0 ## Class: PDict3Parts DNAStringSet logical Dups ## ## Name: elementType elementMetadata metadata ## Class: character DataTableORNULL list ## ## Extends: ## Class "PDict", directly ## Class "List", by class "PDict", distance 2 ## Class "Vector", by class "PDict", distance 3 ## Class "Annotated", by class "PDict", distance 4 ## ## Known Subclasses: "Expanded_TB_PDict"
TB_PDict的类定义看起来比较简单，但理解起来比较费事，尤其是threeparts这个slot。我们从pdict的构造函数PDict来做简单了解，它的使用方式为：
PDict(x, max.mismatch = NA, tb.start = NA, tb.end = NA, tb.width = NA, algorithm = "ACtree2", skip.invalid.patterns = FALSE)
x是为DNAStringSet对象或可以转成DNAStringSet对象的字符串或XStringViews对象。 如果PDict全部使用默认参数，那么x中只能包含确定的碱基，即A/T/G/C；如果要使用模糊匹配，那就得设定tb.start/tb.end/tb.width。 tb即Trusted band，直译为“信任品牌”，按其语言环境可理解为保真区。tb.start/tb.end/tb.width的作用是把查询模式/motif分成三个部分，即TB_PDict对象的threeparts，分别为head（头部），tb（保守区）和tail（尾部），其中tb不能含有简并碱基，而头部和尾部可以，否则出错：
# 这是错误的，因为tb.start和tb.end的设置导致tb含有非DNA碱基xx <- PDict("ATCNGGC", tb.start = 1, tb.end = 4)
## Error: non base DNA letter found in Trusted Band for pattern 1
# 下面是正确的：xx <- PDict("ATCNGGC", tb.start = 2, tb.end = 3)# head部分为Ahead(xx)
## A DNAStringSet instance of length 1 ## width seq ## [1] 1 A
# tb部分为TCtb(xx)
## A DNAStringSet instance of length 1 ## width seq ## [1] 2 TC
# tail部分为NGGCtail(xx)
## A DNAStringSet instance of length 1 ## width seq ## [1] 4 NGGC
pdict可以没有head或tail，但tb不可或缺，tb的长度也不能为0。下面获得的xx对象就没有head：
xx <- PDict(c("ATCNGGC", "ANTCGGCG"), tb.start = 1, tb.width = 1)class(xx)
## [1] "TB_PDict" ## attr(,"package") ## [1] "Biostrings"
# xx的head为NULLhead(xx)
## NULL
tb(xx)
## A DNAStringSet instance of length 2 ## width seq ## [1] 1 A ## [2] 1 A
tail(xx)
## A DNAStringSet instance of length 2 ## width seq ## [1] 6 TCNGGC ## [2] 7 NTCGGCG
PDict用于产生多个匹配模式，但tb.start/tb.end/tb.width这三个参数却分别只能设置一个整数数字，而tb长度又不允许为0，这很怪异。嫌麻烦的话只使用它的全默认设置就可以了。
PDict是预编译的匹配模式，而matchPDict调用C程序进行模式匹配，速度比使用matchPattern循环快得多：
# 产生随机序列的函数rndSeq <- function(dict, n) { paste(sample(dict, n, replace = T), collapse = "")
}## matchPattern和matchPDict运行时比较函数：参数len.p为查询模式的长度，n.p为查询模式的数量comp.runtime <- function(len.p, n.p) {
dna <- DNAString(rndSeq(DNA_BASES, 1e+07))
pat <- mapply(rndSeq, list(DNA_BASES), rep(len.p, n.p))
pd <- PDict(pat)
t.matchPattern <- system.time(for (i in 1:length(pat)) matchPattern(pat[i],
dna))
t.matchPDict <- system.time(matchPDict(pd, dna)) rbind(t.matchPattern, t.matchPDict)[, 1:3]
}# 1个长度为20bp的motif，matchPattern反而比matchPDict要快comp.runtime(20, 1)
## user.self sys.self elapsed ## t.matchPattern 0.02 0 0.01 ## t.matchPDict 0.09 0 0.10
# 10个长度为20bp的motifcomp.runtime(20, 10)
## user.self sys.self elapsed ## t.matchPattern 0.27 0 0.28 ## t.matchPDict 0.09 0 0.09
# 100个长度为20bp的motifcomp.runtime(20, 100)
## user.self sys.self elapsed ## t.matchPattern 2.41 0 2.45 ## t.matchPDict 0.07 0 0.09
# 1000个长度为20bp的motifcomp.runtime(20, 1000)
## user.self sys.self elapsed ## t.matchPattern 25.05 0.04 25.39 ## t.matchPDict 0.17 0.00 0.17
除上面4个函数外，matchPDict系列函数还有whichPDict和vwhichPDict函数，用于获取具有目标序列匹配的模式。但Biostrings的matchPDict系列函数和很多功能还在开发者的TODO list中，尚不完善，要做高通量测序reads的map工作，建议暂时还是先选用其他软件吧。
3、位置权重匹配
转录因子结合位点的序列一般用位置权重矩阵（PWM）、位置频率矩阵（PFM）或位置计数矩阵（PCM）表示。在Biostrings中，PWM由PWM函数使用DNAStringSet构建。PWM矩阵的行名称为'A', 'C', 'G', 'T'， 列数为motif的长度：
motif <- DNAStringSet(c(rep("ACGTGGC", 10), rep("ACGTGTC", 3)))
(pwm <- PWM(motif))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] ## A 0.1442 0.0000 0.0000 0.0000 0.0000 0.00000 0.0000 ## C 0.0000 0.1442 0.0000 0.0000 0.0000 0.00000 0.1442 ## G 0.0000 0.0000 0.1442 0.0000 0.1442 0.13487 0.0000 ## T 0.0000 0.0000 0.0000 0.1442 0.0000 0.09315 0.0000
使用PWM进行序列匹配的函数有：
matchPWM(pwm, subject, min.score = "80%", ...)countPWM(pwm, subject, min.score = "80%", ...)PWMscoreStartingAt(pwm, subject, starting.at = 1)
分别进行PWM匹配、计数和计算特定位置的PWM得分。min.score参数的设置直接影响匹配结果：
(xx <- matchPWM(pwm, DNA.str, min.score = 0.9))
## Views on a 1000-letter DNAString subject ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC ## views: ## start end width ## [1] 17 23 7 [ACGTGTC] ## [2] 65 71 7 [ACGTGTC]
(xx <- matchPWM(pwm, DNA.str, min.score = 0.7))
## Views on a 1000-letter DNAString subject ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC ## views: ## start end width ## [1] 17 23 7 [ACGTGTC] ## [2] 65 71 7 [ACGTGTC] ## [3] 105 111 7 [TCGTTGC] ## [4] 135 141 7 [AAGTGAC] ## [5] 266 272 7 [ATGTGGG] ## [6] 328 334 7 [ACGGTGC] ## [7] 372 378 7 [AATTGGC] ## [8] 641 647 7 [ATGTGGG] ## [9] 711 717 7 [TCGTGCC]
countPWM(pwm, DNA.str, min.score = 0.7)
## [1] 9
PWMscoreStartingAt(pwm, DNA.str, starting.at = start(xx))
## [1] 0.9583 0.9583 0.7116 0.7209 0.7116 0.7116 0.7116 0.7116 0.7209
Bioconductor的另外一个数据包MotifDb收集了两千多个DNA结合的motifs，并有详细的注释。但MotifDb提供的是PFM，要转成PWM还需先使用基因组背景序列计算出各碱基组成。不难，但稍嫌麻烦。
二、序列比对
Biostrings的pairwise比对函数为pairwiseAlignment：
pwa <- pairwiseAlignment("ATCGCAC", "ATCGAAAC", gapOpening = -3, gapExtension = -1)
pwa
## Global PairwiseAlignmentsSingleSubject (1 of 1) ## pattern: [1] ATCG--CAC ## subject: [1] ATCGAA-AC ## score: 2.891
score(pwa)
## [1] 2.891
pattern(pwa)
## [1] ATCG--CAC
subject(pwa)
## [1] ATCGAA-AC
aligned(pwa)
## A BStringSet instance of length 1 ## width seq ## [1] 8 ATCG--AC
另外还有一个很诡异的函数PairwiseAlignments，乍一看长得和pairwiseAlignment还真像。从它对序列长度的要求来看可能是用于score计算：
PairwiseAlignments("-PA--W-HEAE", "HEAGAWGHE-E")
## Global PairwiseAlignments (1 of 1) ## pattern: [1] PA--W-HEAE ## subject: [2] EAGAWGHE-E ## score: -6
PairwiseAlignments("-PA--W-HEAE", "HEAGAWGHE-E", gapOpening = -3, gapExtension = -1)
## Global PairwiseAlignments (1 of 1) ## pattern: [1] PA--W-HEAE ## subject: [2] EAGAWGHE-E ## score: -18
PairwiseAlignments("PAWHEAE", "HEAGAWGHEE")
## Error: 'pattern' and 'subject' must have the same number of characters
多序列比对的函数分DNA/RNA/AA设计了三个函数:
DNAMultipleAlignment(x = character(), start = NA, end = NA, width = NA, use.names = TRUE, rowmask = NULL, colmask = NULL)RNAMultipleAlignment(x = character(), start = NA, end = NA, width = NA, use.names = TRUE, rowmask = NULL, colmask = NULL)AAMultipleAlignment(x = character(), start = NA, end = NA, width = NA, use.names = TRUE, rowmask = NULL, colmask = NULL)
