当前位置:天才代写 > tutorial > R语言教程 > R语言教程之序列处理

R语言教程之序列处理

2018-05-24 08:00 星期四 所属: R语言教程 浏览:731

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)

 

    关键字:

天才代写-代写联系方式