R语言实现基因序列的匹配和比对
我们对字符串都很熟悉，那么面对大量的测序序列字符串，我们如何对其进行处理分析，获得最终的结果。在R语言中有学者专门针对字符串的处理开发了对应的包，命名为Biostrings。
接下来我们介绍下其安装以及主要的功能。
安装还是通过bioconductor进行安装，具体代码如下：
source("https://bioconductor.org/biocLite.R")
biocLite("Biostrings")
接下来生成我们需要的基础数据（摘自CSDN全文地址请点击：https://blog.csdn.net/u014801157/article/details/24372455?utm_source=copy
）：
rndSeq <- function(dict, n) {
paste(sample(dict, n, replace = T), collapse = "")
}
set.seed(0)
# 用mapply和rndSeq函数获取5条序列（字符串）：
DNA.raw <- mapply(rndSeq,list(DNA_BASES), rep(20, 5))
names(DNA.raw) <- paste("SEQ",1:5, sep = "-")
# DNAString对象，1条序列
DNA.str <- DNAString(DNA.raw[1])
# DNAStringSet对象，含5条序列
DNA.set <- DNAStringSet(DNA.raw)
# Views对象
DNA.vws <- successiveViews(DNA.str,width = rep(4, 5))
函数的介绍从此处开始：
1. reverse() 获取反向序列
2. complement() 获取互补的序列
3. reverseComplement() 获取反向互补的序列
4. translate() 翻译函数，他只能针对XString和XXXSet类对象。
XString 类允许我们创建、存储和使用不同类型的字符串。不过我们只被允许使用XString的一些子类: BString, DNAString, RNAString,和AAString.。那么我们首先看下Xstrings类的构建：
bstring = BString("I am a BStringobject")#换成指定的类名称的函数即可。
另外也可以利用另一个函数subseq（），其可以更加方便的构造任意长度的Xstrings。
subseq(bstring, start=1, end=3)
当然我们也可以将Xstrings进行字符串的转化，那么涉及到的函数是toString()。
5. letterFrequency() 获取序列中某些字符的频率。
其中主要的参数as.prob如果为TRUE那么所得的值就是频率，如果FALSE那么为个数。
示例如下：
另一种特殊的用法可能会更有用：
6. letterFrequencyInSlidingView() 函数主要是获取在指定长度序列中各字符的频率，并且将此指定长度作为窗口进行下移一个碱基，直至计算整个序列。示例如下：
letterFrequencyInSlidingView(DNA.set[[1]],view.width = 13, letter = DNA_BASES)
7. alphabetFrequency() 主要是对矩阵中所有的因子进行统计，并列出指定的频率：
接下来我们看下Biostrings中更高级的函数，那就是模式匹配和序列比对。
1. 单模式匹配主要包含以下函数：
matchPattern()：1个查询模式1条序列
countPattern()：1个查询模式1条序列，仅计数
vmatchPattern()：1个查询模式n条序列
vcountPattern()：1个查询模式n条序列，仅计数
2. 多模式的匹配函数如下：
matchPDict()：n个查询模式1条序列
countPDict()：n个查询模式1条序列，仅计数
vmatchPDict()：n个查询模式n条序列
vcountPDict()：n个查询模式n条序列，仅计数
首先我们导入我们需要的数据包：
source("https://bioconductor.org/biocLite.R")
biocLite("drosophila2probe")
biocLite("BSgenome.Dmelanogaster.UCSC.dm3")
利用以上数据构建我们的字典以及匹配数据：
接下来看我们的实例：
mi0 <- matchPDict(pdict0, chr3R)
注：我们上面所提到的所谓模式也就是指的序列的reads。
3. PWM() 位置频率矩阵计算。具体的实例如下：
motif <-DNAStringSet(c(rep("ACGTGGC", 10), rep("ACGTGTC", 3)))
pwm <- PWM(motif)
4. pairwiseAlignment() / PairwiseAlignments()序列的对比，实例如下：
pwa <-pairwiseAlignment("-PA--W-HEAE", "HEAGAWGHE-E", gapOpening= -3, gapExtension = -1)
PairwiseAlignments("-PA--W-HEAE","HEAGAWGHE-E", gapOpening = -3, gapExtension = -1)
