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

评论3,712

Biostrings最后一节，介绍模式匹配和序列比对的相关函数和操作。

下面我们使用拟南芥基因转录起始点上游1kb的序列进行分析。序列文件可以从TAIR网站（http://www.arabidopsis.org）下载。先用readDNAStringSet函数从FASTA文件中读取序列并查看头2个序列的信息：

  1. library(Biostrings)
  2. upstream.1k <- readDNAStringSet(file.choose(), "fasta")
  3. head(upstream.1k, 2)
  4. ## A DNAStringSet instance of length 2
  5. ## width seq names
  6. ## [1] 1000 ACAAGCATGCTTAGCATACGT...ATCTTAATGTGGATAGTGCT AT1G08520 | chr1:...
  7. ## [2] 1000 CGTAACAGTAACAACTATATT...CTCCTGGAGAAGAGAAGACT AT1G08530 | chr1:...
  8. # 获取第29条序列
  9. 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. #
  2. # 1模式查询1条序列：非模糊匹配，不能使用IUPAC简并代码，所以结果找不到（views:
  3. # NONE）
  4. matchPattern("ACGTGKC", DNA.str)
  5. ## Views on a 1000-letter DNAString subject
  6. ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC
  7. ## views: NONE
  8. # 模糊匹配，使用IUPAC简并代码，K可匹配G/T
  9. matchPattern("ACGTGKC", DNA.str, fixed = FALSE)
  10. ## Views on a 1000-letter DNAString subject
  11. ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC
  12. ## views:
  13. ## start end width
  14. ## [1] 17 23 7 [ACGTGTC]
  15. ## [2] 65 71 7 [ACGTGTC]
  16. # 仅计数
  17. countPattern("ACGTGKC", DNA.str, fixed = FALSE)
  18. ## [1] 2

fixed参数设置为FALSE时，如果目标序列上有非确定碱基如N，得到的结果会很让人失望：

  1. # 1模式查询n条序列，仅计数：
  2. vc <- vcountPattern("ACGTGKC", upstream.1k, fixed = FALSE)
  3. # 看看最多匹配数，太多了
  4. max(vc)
  5. ## [1] 970
  6. # 用subset函数提取具有最多匹配数的序列子集
  7. subs <- subset(upstream.1k, vc == max(vc))
  8. # 原来匹配的都是NNN...，怪不得那么多
  9. matchPattern("ACGTGKC", subs[[1]], fixed = FALSE)
  10. ## Views on a 1000-letter DNAString subject
  11. ## subject: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNGAATTCGTCGACCAGGACGGCGGA
  12. ## views:
  13. ## start end width
  14. ## [1] 1 7 7 [NNNNNNN]
  15. ## [2] 2 8 7 [NNNNNNN]
  16. ## [3] 3 9 7 [NNNNNNN]
  17. ## [4] 4 10 7 [NNNNNNN]
  18. ## [5] 5 11 7 [NNNNNNN]
  19. ## [6] 6 12 7 [NNNNNNN]
  20. ## [7] 7 13 7 [NNNNNNN]
  21. ## [8] 8 14 7 [NNNNNNN]
  22. ## [9] 9 15 7 [NNNNNNN]
  23. ## ... ... ... ... ...
  24. ## [962] 962 968 7 [NNNNNNN]
  25. ## [963] 963 969 7 [NNNNNNN]
  26. ## [964] 964 970 7 [NNNNNNN]
  27. ## [965] 965 971 7 [NNNNNNN]
  28. ## [966] 966 972 7 [NNNNNNN]
  29. ## [967] 967 973 7 [NNNNNNN]
  30. ## [968] 968 974 7 [NNNNNNN]
  31. ## [969] 969 975 7 [NNNNNNN]
  32. ## [970] 970 976 7 [NNNNNNN]

从上面的结果可以看出，如果fixed=FALSE，模式和目标序列的碱基都要应用IUPAC简并规则。但实际应用是我们往往要忽略目标序列上的不确定碱基，这时应设置fixed="subject"，固定目标序列上的碱基：

  1. vc <- vcountPattern("ACGTGKC", upstream.1k, fixed = "subject")
  2. subs <- subset(upstream.1k, vc == max(vc))
  3. # 查看结果，这才是我们想要的：
  4. matchPattern("ACGTGKC", subs[[1]], fixed = "subject")
  5. ## Views on a 1000-letter DNAString subject
  6. ## subject: TTACCCGAGTTGGGTAATCGACTCGGTTGGTT...CCTTTTCATTTTCTTCTCCCTCTCCTGGGTT
  7. ## views:
  8. ## start end width
  9. ## [1] 171 177 7 [ACGTGGC]
  10. ## [2] 350 356 7 [ACGTGTC]
  11. ## [3] 421 427 7 [ACGTGTC]
  12. ## [4] 449 455 7 [ACGTGGC]

上面语句如果往清楚里写应是这样的：

  1. matchPattern("ACGTGKC", subs[[1]], fixed = c(subject = TRUE, pattern = FALSE))
  2. ## Views on a 1000-letter DNAString subject
  3. ## subject: TTACCCGAGTTGGGTAATCGACTCGGTTGGTT...CCTTTTCATTTTCTTCTCCCTCTCCTGGGTT
  4. ## views:
  5. ## start end width
  6. ## [1] 171 177 7 [ACGTGGC]
  7. ## [2] 350 356 7 [ACGTGTC]
  8. ## [3] 421 427 7 [ACGTGTC]
  9. ## [4] 449 455 7 [ACGTGGC]

模式匹配的结果可进一步使用，比如要获取翻译起始点上游有2个以上ABRE motif的序列名称，这些基因可能对ABA有较强的响应：

  1. names(subset(upstream.1k, vc > 2))
  2. ## [1] "AT1G24580 | chr1:8709180-8710179 FORWARD LENGTH=1000"
  3. ## [2] "AT1G09210 | chr1:2976751-2977750 REVERSE LENGTH=1000"
  4. ## [3] "AT1G19490 | chr1:6754000-6754999 REVERSE LENGTH=1000"
  5. ## [4] "AT1G07870 | chr1:2432058-2433057 REVERSE LENGTH=1000"
  6. ## [5] "AT1G51140 | chr1:18945754-18946753 REVERSE LENGTH=1000"
  7. ## [6] "AT1G20180 | chr1:6995210-6996209 FORWARD LENGTH=1000"
  8. ## [7] "AT1G77450 | chr1:29098954-29099953 FORWARD LENGTH=1000"
  9. ## [8] "AT1G72510 | chr1:27302366-27303365 FORWARD LENGTH=1000"
  10. ## [9] "AT1G60190 | chr1:22197403-22198402 FORWARD LENGTH=1000"
  11. ## [10] "AT1G79040 | chr1:29735016-29736015 FORWARD LENGTH=1000"
  12. ## [11] "AT1G55520 | chr1:20728092-20729091 REVERSE LENGTH=1000"
  13. ## [12] "AT2G04350 | chr2:1514805-1515804 FORWARD LENGTH=1000"
  14. ## [13] "AT2G22240 | chr2:9454095-9455094 REVERSE LENGTH=1000"
  15. ## [14] "AT3G63350 | chr3:23398468-23399467 FORWARD LENGTH=1000"
  16. ## [15] "AT3G33154 | chr3:13964599-13965598 FORWARD LENGTH=1000"
  17. ## [16] "AT3G42065 | chr3:14254659-14255658 FORWARD LENGTH=1000"
  18. ## [17] "AT3G63070 | chr3:23301365-23302364 FORWARD LENGTH=1000"
  19. ## [18] "AT3G49440 | chr3:18336318-18337317 REVERSE LENGTH=1000"
  20. ## [19] "AT3G10985 | chr3:3441528-3442527 FORWARD LENGTH=1000"
  21. ## [20] "AT3G15280 | chr3:5144694-5145693 REVERSE LENGTH=1000"
  22. ## [21] "AT3G03680 | chr3:906500-907499 FORWARD LENGTH=1000"
  23. ## [22] "AT3G02140 | chr3:386539-387538 REVERSE LENGTH=1000"
  24. ## [23] "AT3G53040 | chr3:19666488-19667487 REVERSE LENGTH=1000"
  25. ## [24] "AT4G21940 | chr4:11639802-11640801 FORWARD LENGTH=1000"
  26. ## [25] "AT4G25570 | chr4:13055624-13056623 REVERSE LENGTH=1000"
  27. ## [26] "AT4G27410 | chr4:13709150-13710149 REVERSE LENGTH=1000"
  28. ## [27] "AT5G50360 | chr5:20506339-20507338 REVERSE LENGTH=1000"
  29. ## [28] "AT5G10490 | chr5:3304762-3305761 REVERSE LENGTH=1000"
  30. ## [29] "AT5G15950 | chr5:5204869-5205868 FORWARD LENGTH=1000"
  31. ## [30] "AT5G28667 | chr5:10689062-10690061 FORWARD LENGTH=1000"
  32. ## [31] "AT5G06980 | chr5:2166468-2167467 FORWARD LENGTH=1000"
  33. ## [32] "AT5G15948 | chr5:5204871-5205870 FORWARD LENGTH=1000"

掩膜序列生成函数maskMotif调用的是matchPattern函数，可以在参数中设置matchPattern的参数。下面两条语句的结果等价：

  1. as(maskMotif(DNA.str, "ACGTGKC", fixed = "subject"), "Views")
  2. ## Views on a 1000-letter DNAString subject
  3. ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC
  4. ## views:
  5. ## start end width
  6. ## [1] 1 16 16 [AGAGAAGGCTCCAAAC]
  7. ## [2] 24 64 41 [TCAGCCATATTATATCATTCGTCCGAGAGAAGGCTCCACAC]
  8. ## [3] 72 1000 929 [TCAGCCATATGTCTATCCGACGTAA...CTTTGCTCATTATATGTCCTCAGC]
  9. gaps(matchPattern("ACGTGKC", DNA.str, fixed = "subject"))
  10. ## Views on a 1000-letter DNAString subject
  11. ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC
  12. ## views:
  13. ## start end width
  14. ## [1] 1 16 16 [AGAGAAGGCTCCAAAC]
  15. ## [2] 24 64 41 [TCAGCCATATTATATCATTCGTCCGAGAGAAGGCTCCACAC]
  16. ## [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而来：

  1. getClass("TB_PDict")
  2. ## Class "TB_PDict" [package "Biostrings"]
  3. ##
  4. ## Slots:
  5. ##
  6. ## Name: threeparts dict0 constant_width dups0
  7. ## Class: PDict3Parts DNAStringSet logical Dups
  8. ##
  9. ## Name: elementType elementMetadata metadata
  10. ## Class: character DataTableORNULL list
  11. ##
  12. ## Extends:
  13. ## Class "PDict", directly
  14. ## Class "List", by class "PDict", distance 2
  15. ## Class "Vector", by class "PDict", distance 3
  16. ## Class "Annotated", by class "PDict", distance 4
  17. ##
  18. ## Known Subclasses: "Expanded_TB_PDict"

TB_PDict的类定义看起来比较简单，但理解起来比较费事，尤其是threeparts这个slot。我们从pdict的构造函数PDict来做简单了解，它的使用方式为：

  1. PDict(x, max.mismatch = NA, tb.start = NA, tb.end = NA, tb.width = NA, algorithm = "ACtree2",
  2. 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不能含有简并碱基，而头部和尾部可以，否则出错：

  1. # 这是错误的，因为tb.start和tb.end的设置导致tb含有非DNA碱基
  2. xx <- PDict("ATCNGGC", tb.start = 1, tb.end = 4)
  3. ## Error: non base DNA letter found in Trusted Band for pattern 1
  4. # 下面是正确的：
  5. xx <- PDict("ATCNGGC", tb.start = 2, tb.end = 3)
  6. # head部分为A
  7. head(xx)
  8. ## A DNAStringSet instance of length 1
  9. ## width seq
  10. ## [1] 1 A
  11. # tb部分为TC
  12. tb(xx)
  13. ## A DNAStringSet instance of length 1
  14. ## width seq
  15. ## [1] 2 TC
  16. # tail部分为NGGC
  17. tail(xx)
  18. ## A DNAStringSet instance of length 1
  19. ## width seq
  20. ## [1] 4 NGGC
  21. pdict可以没有head或tail，但tb不可或缺，tb的长度也不能为0。下面获得的xx对象就没有head：
  22. xx <- PDict(c("ATCNGGC", "ANTCGGCG"), tb.start = 1, tb.width = 1)
  23. class(xx)
  24. ## [1] "TB_PDict"
  25. ## attr(,"package")
  26. ## [1] "Biostrings"
  27. # xx的head为NULL
  28. head(xx)
  29. ## NULL
  30. tb(xx)
  31. ## A DNAStringSet instance of length 2
  32. ## width seq
  33. ## [1] 1 A
  34. ## [2] 1 A
  35. tail(xx)
  36. ## A DNAStringSet instance of length 2
  37. ## width seq
  38. ## [1] 6 TCNGGC
  39. ## [2] 7 NTCGGCG

PDict用于产生多个匹配模式，但tb.start/tb.end/tb.width这三个参数却分别只能设置一个整数数字，而tb长度又不允许为0，这很怪异。嫌麻烦的话只使用它的全默认设置就可以了。

PDict是预编译的匹配模式，而matchPDict调用C程序进行模式匹配，速度比使用matchPattern循环快得多：

  1. # 产生随机序列的函数
  2. rndSeq <- function(dict, n) {
  3. paste(sample(dict, n, replace = T), collapse = "")
  4. }
  5. #
  6. # matchPattern和matchPDict运行时比较函数：参数len.p为查询模式的长度，n.p为查询模式的数量
  7. comp.runtime <- function(len.p, n.p) {
  8. dna <- DNAString(rndSeq(DNA_BASES, 1e+07))
  9. pat <- mapply(rndSeq, list(DNA_BASES), rep(len.p, n.p))
  10. pd <- PDict(pat)
  11. t.matchPattern <- system.time(for (i in 1:length(pat)) matchPattern(pat[i],
  12. dna))
  13. t.matchPDict <- system.time(matchPDict(pd, dna))
  14. rbind(t.matchPattern, t.matchPDict)[, 1:3]
  15. }
  16. # 1个长度为20bp的motif，matchPattern反而比matchPDict要快
  17. comp.runtime(20, 1)
  18. ## user.self sys.self elapsed
  19. ## t.matchPattern 0.02 0 0.01
  20. ## t.matchPDict 0.09 0 0.10
  21. # 10个长度为20bp的motif
  22. comp.runtime(20, 10)
  23. ## user.self sys.self elapsed
  24. ## t.matchPattern 0.27 0 0.28
  25. ## t.matchPDict 0.09 0 0.09
  26. # 100个长度为20bp的motif
  27. comp.runtime(20, 100)
  28. ## user.self sys.self elapsed
  29. ## t.matchPattern 2.41 0 2.45
  30. ## t.matchPDict 0.07 0 0.09
  31. # 1000个长度为20bp的motif
  32. comp.runtime(20, 1000)
  33. ## user.self sys.self elapsed
  34. ## t.matchPattern 25.05 0.04 25.39
  35. ## 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的长度：

  1. motif <- DNAStringSet(c(rep("ACGTGGC", 10), rep("ACGTGTC", 3)))
  2. (pwm <- PWM(motif))
  3. ## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
  4. ## A 0.1442 0.0000 0.0000 0.0000 0.0000 0.00000 0.0000
  5. ## C 0.0000 0.1442 0.0000 0.0000 0.0000 0.00000 0.1442
  6. ## G 0.0000 0.0000 0.1442 0.0000 0.1442 0.13487 0.0000
  7. ## T 0.0000 0.0000 0.0000 0.1442 0.0000 0.09315 0.0000

使用PWM进行序列匹配的函数有：

  1. matchPWM(pwm, subject, min.score = "80%", ...)
  2. countPWM(pwm, subject, min.score = "80%", ...)
  3. PWMscoreStartingAt(pwm, subject, starting.at = 1)

分别进行PWM匹配、计数和计算特定位置的PWM得分。min.score参数的设置直接影响匹配结果：

  1. (xx <- matchPWM(pwm, DNA.str, min.score = 0.9))
  2. ## Views on a 1000-letter DNAString subject
  3. ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC
  4. ## views:
  5. ## start end width
  6. ## [1] 17 23 7 [ACGTGTC]
  7. ## [2] 65 71 7 [ACGTGTC]
  8. (xx <- matchPWM(pwm, DNA.str, min.score = 0.7))
  9. ## Views on a 1000-letter DNAString subject
  10. ## subject: AGAGAAGGCTCCAAACACGTGTCTCAGCCATA...AATAAATCTTTGCTCATTATATGTCCTCAGC
  11. ## views:
  12. ## start end width
  13. ## [1] 17 23 7 [ACGTGTC]
  14. ## [2] 65 71 7 [ACGTGTC]
  15. ## [3] 105 111 7 [TCGTTGC]
  16. ## [4] 135 141 7 [AAGTGAC]
  17. ## [5] 266 272 7 [ATGTGGG]
  18. ## [6] 328 334 7 [ACGGTGC]
  19. ## [7] 372 378 7 [AATTGGC]
  20. ## [8] 641 647 7 [ATGTGGG]
  21. ## [9] 711 717 7 [TCGTGCC]
  22. countPWM(pwm, DNA.str, min.score = 0.7)
  23. ## [1] 9
  24. PWMscoreStartingAt(pwm, DNA.str, starting.at = start(xx))
  25. ## [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：

  1. pwa <- pairwiseAlignment("ATCGCAC", "ATCGAAAC", gapOpening = -3, gapExtension = -1)
  2. pwa
  3. ## Global PairwiseAlignmentsSingleSubject (1 of 1)
  4. ## pattern: [1] ATCG--CAC
  5. ## subject: [1] ATCGAA-AC
  6. ## score: 2.891
  7. score(pwa)
  8. ## [1] 2.891
  9. pattern(pwa)
  10. ## [1] ATCG--CAC
  11. subject(pwa)
  12. ## [1] ATCGAA-AC
  13. aligned(pwa)
  14. ## A BStringSet instance of length 1
  15. ## width seq
  16. ## [1] 8 ATCG--AC

另外还有一个很诡异的函数PairwiseAlignments，乍一看长得和pairwiseAlignment还真像。从它对序列长度的要求来看可能是用于score计算：

  1. PairwiseAlignments("-PA--W-HEAE", "HEAGAWGHE-E")
  2. ## Global PairwiseAlignments (1 of 1)
  3. ## pattern: [1] PA--W-HEAE
  4. ## subject: [2] EAGAWGHE-E
  5. ## score: -6
  6. PairwiseAlignments("-PA--W-HEAE", "HEAGAWGHE-E", gapOpening = -3, gapExtension = -1)
  7. ## Global PairwiseAlignments (1 of 1)
  8. ## pattern: [1] PA--W-HEAE
  9. ## subject: [2] EAGAWGHE-E
  10. ## score: -18
  11. PairwiseAlignments("PAWHEAE", "HEAGAWGHEE")
  12. ## Error: 'pattern' and 'subject' must have the same number of characters

多序列比对的函数分DNA/RNA/AA设计了三个函数:

  1. DNAMultipleAlignment(x = character(), start = NA, end = NA, width = NA, use.names = TRUE,
  2. rowmask = NULL, colmask = NULL)
  3. RNAMultipleAlignment(x = character(), start = NA, end = NA, width = NA, use.names = TRUE,
  4. rowmask = NULL, colmask = NULL)
  5. AAMultipleAlignment(x = character(), start = NA, end = NA, width = NA, use.names = TRUE,
  6. rowmask = NULL, colmask = NULL)

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

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

发表评论

匿名网友

拖动滑块以完成验证
加载失败