R/BioC序列处理之三：Biostrings模式匹配和序列比对

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/T
matchPattern("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部分为A
head(xx)
## A DNAStringSet instance of length 1
## width seq
## [1] 1 A
# tb部分为TC
tb(xx)
## A DNAStringSet instance of length 1
## width seq
## [1] 2 TC
# tail部分为NGGC
tail(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为NULL
head(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的motif
comp.runtime(20, 10)
## user.self sys.self elapsed
## t.matchPattern 0.27 0 0.28
## t.matchPDict 0.09 0 0.09
# 100个长度为20bp的motif
comp.runtime(20, 100)
## user.self sys.self elapsed
## t.matchPattern 2.41 0 2.45
## t.matchPDict 0.07 0 0.09
# 1000个长度为20bp的motif
comp.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)

序列比对在其他软件做得太多了，应该都会，不再啰嗦。

原文来自：http://blog.csdn.net/u014801157/article/details/24372465