R/BioC序列处理之二：Biostrings序列的基本操作

评论4,975

还是先获取随机DNA序列和其他序列对象：

  1. library(Biostrings)
  2. rndSeq <- function(dict, n) {
  3. paste(sample(dict, n, replace = T), collapse = "")
  4. }
  5. set.seed(0)
  6. # 用mapply和rndSeq函数获取5条序列（字符串）：
  7. DNA.raw <- mapply(rndSeq, list(DNA_BASES), rep(20, 5))
  8. names(DNA.raw) <- paste("SEQ", 1:5, sep = "-")
  9. # DNAString对象，1条序列
  10. DNA.str <- DNAString(DNA.raw[1])
  11. # DNAStringSet对象，含5条序列
  12. DNA.set <- DNAStringSet(DNA.raw)
  13. # Views对象
  14. DNA.vws <- successiveViews(DNA.str, width = rep(4, 5))

一、获取序列基本信息

包括获取名称（names）、长度（length）、字符个数（nchar）和对象头/尾（head/tail）等信息的函数。

函数的用法简单，但需注意XString类对象的返回结果和其他类型有些差别：

  1. # length函数----------------------
  2. length(DNA.raw)
  3. ## [1] 5
  4. # 结果为序列的长度
  5. length(DNA.str)
  6. ## [1] 20
  7. length(DNA.set)
  8. ## [1] 5
  9. length(DNA.vws)
  10. ## [1] 5
  11. # nchar函数-----------------------
  12. nchar(DNA.raw)
  13. ## SEQ-1 SEQ-2 SEQ-3 SEQ-4 SEQ-5
  14. ## 20 20 20 20 20
  15. nchar(DNA.str)
  16. ## [1] 20
  17. nchar(DNA.set)
  18. ## [1] 20 20 20 20 20
  19. nchar(DNA.vws)
  20. ## [1] 4 4 4 4 4
  21. # head/tail函数-------------------
  22. head(DNA.raw, 2)
  23. ## SEQ-1 SEQ-2
  24. ## "TCCGTATTGGAAAGCTCGTC" "TTAGACCACTCCGCATGTAG"
  25. # 结果为序列的前几个碱基
  26. head(DNA.str, 2)
  27. ## 2-letter "DNAString" instance
  28. ## seq: TC
  29. head(DNA.set, 2)
  30. ## A DNAStringSet instance of length 2
  31. ## width seq names
  32. ## [1] 20 TCCGTATTGGAAAGCTCGTC SEQ-1
  33. ## [2] 20 TTAGACCACTCCGCATGTAG SEQ-2
  34. head(DNA.vws, 2)
  35. ## Views on a 20-letter DNAString subject
  36. ## subject: TCCGTATTGGAAAGCTCGTC
  37. ## views:
  38. ## start end width
  39. ## [1] 1 4 4 [TCCG]
  40. ## [2] 5 8 4 [TATT]

二、序列转换

1、获取反向、互补、反向互补序列：

reverse(), complement(), reverseComplement()可以使用string, XString, XXXSet和Views类对象进行操作。下面看看reverse函数的结果：

  1. DNA.raw
  2. ## SEQ-1 SEQ-2 SEQ-3
  3. ## "TCCGTATTGGAAAGCTCGTC" "TTAGACCACTCCGCATGTAG" "CTGTGGTACGGCTCAAACGG"
  4. ## SEQ-4 SEQ-5
  5. ## "CTCCCGCCTATCTCCCTTCT" "TCGCCTAGAAAAAGTTTCCT"
  6. reverse(DNA.raw)
  7. ## SEQ-1 SEQ-2 SEQ-3
  8. ## "CTGCTCGAAAGGTTATGCCT" "GATGTACGCCTCACCAGATT" "GGCAAACTCGGCATGGTGTC"
  9. ## SEQ-4 SEQ-5
  10. ## "TCTTCCCTCTATCCGCCCTC" "TCCTTTGAAAAAGATCCGCT"
  11. reverse(DNA.str)
  12. ## 20-letter "DNAString" instance
  13. ## seq: CTGCTCGAAAGGTTATGCCT
  14. reverse(DNA.set)
  15. ## A DNAStringSet instance of length 5
  16. ## width seq names
  17. ## [1] 20 CTGCTCGAAAGGTTATGCCT SEQ-1
  18. ## [2] 20 GATGTACGCCTCACCAGATT SEQ-2
  19. ## [3] 20 GGCAAACTCGGCATGGTGTC SEQ-3
  20. ## [4] 20 TCTTCCCTCTATCCGCCCTC SEQ-4
  21. ## [5] 20 TCCTTTGAAAAAGATCCGCT SEQ-5
  22. reverse(DNA.vws)
  23. ## Views on a 20-letter DNAString subject
  24. ## subject: CTGCTCGAAAGGTTATGCCT
  25. ## views:
  26. ## start end width
  27. ## [1] 17 20 4 [GCCT]
  28. ## [2] 13 16 4 [TTAT]
  29. ## [3] 9 12 4 [AAGG]
  30. ## [4] 5 8 4 [TCGA]
  31. ## [5] 1 4 4 [CTGC]

2、序列翻译：

翻译函数translate()只能使用XString和XXXSet类对象（Biostrings version 2.29.3）：

  1. # 错误
  2. translate(DNA.raw)
  3. ## Error: unable to find an inherited method for function 'translate' for
  4. ## signature '"character"'
  5. translate(DNA.str)
  6. ## 6-letter "AAString" instance
  7. ## seq: SVLESS
  8. # 错误
  9. translate(DNA.set)
  10. ## A AAStringSet instance of length 5
  11. ## width seq
  12. ## [1] 6 SVLESS
  13. ## [2] 6 LDHSAC
  14. ## [3] 6 LWYGSN
  15. ## [4] 6 LPPISL
  16. ## [5] 6 SPRKSF
  17. # 错误
  18. translate(DNA.vws)
  19. ## Error: unable to find an inherited method for function 'translate' for
  20. ## signature '"XStringViews"'

3、DNA/RNA互转：

使用dna2rna, rna2dna或cDNA函数。这些函数对于数据对象有严格的要求：

  1. DNA.str
  2. ## 20-letter "DNAString" instance
  3. ## seq: TCCGTATTGGAAAGCTCGTC
  4. dna2rna(DNA.str)
  5. ## 20-letter "RNAString" instance
  6. ## seq: UCCGUAUUGGAAAGCUCGUC
  7. # 错误
  8. dna2rna(DNA.raw)
  9. ## Error: dna2rna() only works on DNA input
  10. # 错误
  11. cDNA(DNA.str)
  12. ## Error: cDNA() only works on RNA input
  13. cDNA(dna2rna(DNA.str))
  14. ## 20-letter "DNAString" instance
  15. ## seq: AGGCATAACCTTTCGAGCAG

4、其他：

R base包的chartr函数已经重载，可以应用于序列对象，但最好避免使用，因为会出现符号检查问题：

  1. chartr("T", "A", DNA.set)
  2. ## A DNAStringSet instance of length 5
  3. ## width seq names
  4. ## [1] 20 ACCGAAAAGGAAAGCACGAC SEQ-1
  5. ## [2] 20 AAAGACCACACCGCAAGAAG SEQ-2
  6. ## [3] 20 CAGAGGAACGGCACAAACGG SEQ-3
  7. ## [4] 20 CACCCGCCAAACACCCAACA SEQ-4
  8. ## [5] 20 ACGCCAAGAAAAAGAAACCA SEQ-5
  9. # 错误
  10. chartr("T", "U", DNA.set)
  11. ## Error: key 85 (char 'U') not in lookup table

其他一些R base包的字符串操作函数也可以使用序列对象，但注意返回结果的类型。下面两行代码的返回值都是字符串向量，而不是DNAStringSet对象：

  1. tolower(DNA.set)
  2. ## SEQ-1 SEQ-2 SEQ-3
  3. ## "tccgtattggaaagctcgtc" "ttagaccactccgcatgtag" "ctgtggtacggctcaaacgg"
  4. ## SEQ-4 SEQ-5
  5. ## "ctcccgcctatctcccttct" "tcgcctagaaaaagtttcct"
  6. gsub("T", "U", DNA.set)
  7. ## SEQ-1 SEQ-2 SEQ-3
  8. ## "UCCGUAUUGGAAAGCUCGUC" "UUAGACCACUCCGCAUGUAG" "CUGUGGUACGGCUCAAACGG"
  9. ## SEQ-4 SEQ-5
  10. ## "CUCCCGCCUAUCUCCCUUCU" "UCGCCUAGAAAAAGUUUCCU"

三、序列截取（sequence subset）

1、使用序列构造函数

Biostrings的XXXSet类和Views类序列构造函数本身就具备序列subset的功能，Views类对象可以通过类型转换获得XXXSet类对象。具体使用方法请参看上一篇文章。

2、subseq函数

subseq函数可以用于序列截取，也可以对选定序列进行修改：

  1. DNAstr <- DNAString(rndSeq(DNA_BASES, 50))
  2. (subseq(DNAstr, start = 20, end = 30))
  3. ## 11-letter "DNAString" instance
  4. ## seq: CGTCCATCGAA
  5. # 上面语句的subseq函数仅仅是截取了序列，没有改变原序列
  6. DNAstr
  7. ## 50-letter "DNAString" instance
  8. ## seq: GGCCTGAACTGTGCCAAGACGTCCATCGAAGGAAGTGGGTGGGACGCAGA
  9. # 下面语句的subseq函数改变了原序列
  10. subseq(DNAstr, start = 20, end = 30) <- DNAString("NNN")
  11. DNAstr
  12. ## 42-letter "DNAString" instance
  13. ## seq: GGCCTGAACTGTGCCAAGANNNGGAAGTGGGTGGGACGCAGA

3、特殊序列查找和截取

回文（palindrome）序列相关的函数，有：

  1. findPalindromes(subject, min.armlength = 4, max.looplength = 1, min.looplength = 0,
  2. max.mismatch = 0)
  3. palindromeArmLength(x, max.mismatch = 0, ...)
  4. palindromeLeftArm(x, max.mismatch = 0, ...)
  5. palindromeRightArm(x, max.mismatch = 0, ...)
  6. findComplementedPalindromes(subject, min.armlength = 4, max.looplength = 1,
  7. min.looplength = 0, max.mismatch = 0)
  8. complementedPalindromeArmLength(x, max.mismatch = 0, ...)
  9. complementedPalindromeLeftArm(x, max.mismatch = 0, ...)
  10. complementedPalindromeRightArm(x, max.mismatch = 0, ...)

有应用需求的自己去看看函数说明吧。

四、字符或寡核苷酸组合的统计

1、letterFrequency函数

这个函数用于统计指定字符的频率或比例：

  1. letterFrequency(DNA.set, DNA_BASES)
  2. ## A C G T
  3. ## [1,] 4 5 5 6
  4. ## [2,] 5 6 4 5
  5. ## [3,] 4 5 7 4
  6. ## [4,] 1 11 1 7
  7. ## [5,] 6 5 3 6
  8. letterFrequency(DNA.set, DNA_ALPHABET)
  9. ## A C G T M R W S Y K V H D B N - +
  10. ## [1,] 4 5 5 6 0 0 0 0 0 0 0 0 0 0 0 0 0
  11. ## [2,] 5 6 4 5 0 0 0 0 0 0 0 0 0 0 0 0 0
  12. ## [3,] 4 5 7 4 0 0 0 0 0 0 0 0 0 0 0 0 0
  13. ## [4,] 1 11 1 7 0 0 0 0 0 0 0 0 0 0 0 0 0
  14. ## [5,] 6 5 3 6 0 0 0 0 0 0 0 0 0 0 0 0 0
  15. letterFrequency(DNA.set, DNA_BASES, as.prob = TRUE)
  16. ## A C G T
  17. ## [1,] 0.20 0.25 0.25 0.30
  18. ## [2,] 0.25 0.30 0.20 0.25
  19. ## [3,] 0.20 0.25 0.35 0.20
  20. ## [4,] 0.05 0.55 0.05 0.35
  21. ## [5,] 0.30 0.25 0.15 0.30
  22. # 注意下面这种用法：
  23. letterFrequency(DNA.set, "GC", as.prob = TRUE)
  24. ## G|C
  25. ## [1,] 0.5
  26. ## [2,] 0.5
  27. ## [3,] 0.6
  28. ## [4,] 0.6
  29. ## [5,] 0.4

2、letterFrequencyInSlidingView函数

该函数按设置的窗口长度（view.width）一个个碱基滑动并统计字符频率或比例：

  1. letterFrequencyInSlidingView(DNA.set[[1]], view.width = 16, letter = DNA_BASES)
  2. ## A C G T
  3. ## [1,] 4 3 4 5
  4. ## [2,] 4 4 4 4
  5. ## [3,] 4 3 5 4
  6. ## [4,] 4 2 5 5
  7. ## [5,] 4 3 4 5
  8. letterFrequencyInSlidingView(DNA.set[[1]], 16, "GC", as.prob = TRUE)
  9. ## G|C
  10. ## [1,] 0.4375
  11. ## [2,] 0.5000
  12. ## [3,] 0.5000
  13. ## [4,] 0.4375
  14. ## [5,] 0.4375

3、alphabetFrequency函数

作用与letterFrequency函数类似，但按ALPHABET中的所有因子进行统计。baseOnly设置为TRUE可以对ALPHABET进行限制：

  1. alphabetFrequency(DNA.set)
  2. ## A C G T M R W S Y K V H D B N - +
  3. ## [1,] 4 5 5 6 0 0 0 0 0 0 0 0 0 0 0 0 0
  4. ## [2,] 5 6 4 5 0 0 0 0 0 0 0 0 0 0 0 0 0
  5. ## [3,] 4 5 7 4 0 0 0 0 0 0 0 0 0 0 0 0 0
  6. ## [4,] 1 11 1 7 0 0 0 0 0 0 0 0 0 0 0 0 0
  7. ## [5,] 6 5 3 6 0 0 0 0 0 0 0 0 0 0 0 0 0
  8. alphabetFrequency(DNA.set, as.prob = TRUE, baseOnly = TRUE)
  9. ## A C G T other
  10. ## [1,] 0.20 0.25 0.25 0.30 0
  11. ## [2,] 0.25 0.30 0.20 0.25 0
  12. ## [3,] 0.20 0.25 0.35 0.20 0
  13. ## [4,] 0.05 0.55 0.05 0.35 0
  14. ## [5,] 0.30 0.25 0.15 0.30 0

4、寡核苷酸统计：

有以下函数，分别统计2核苷酸组合、3核苷酸组合和寡核苷酸组合：

dinucleotideFrequency

trinucleotideFrequency

oligonucleotideFrequency

五、杂项函数

包括使用序列对象进行运算的一些函数（和符号）：

  1. # 反转对象元素的顺序，不是反转序列！
  2. rev(DNA.set)
  3. ## A DNAStringSet instance of length 5
  4. ## width seq names
  5. ## [1] 20 TCGCCTAGAAAAAGTTTCCT SEQ-5
  6. ## [2] 20 CTCCCGCCTATCTCCCTTCT SEQ-4
  7. ## [3] 20 CTGTGGTACGGCTCAAACGG SEQ-3
  8. ## [4] 20 TTAGACCACTCCGCATGTAG SEQ-2
  9. ## [5] 20 TCCGTATTGGAAAGCTCGTC SEQ-1
  10. # 判断两个对象是否相同
  11. DNA.set[[1]] == DNA.set[[2]]
  12. ## [1] FALSE
  13. # 判断对象是否包含在另外一个对象的元素中
  14. DNA.str %in% DNA.set
  15. ## [1] TRUE

Biostrings在2.0版后还添加了append函数，可以把几个序列集合合并，还可以使用c函数进行合并：

  1. append(DNA.set, DNAStringSet(DNA.str))
  2. ## A DNAStringSet instance of length 6
  3. ## width seq names
  4. ## [1] 20 TCCGTATTGGAAAGCTCGTC SEQ-1
  5. ## [2] 20 TTAGACCACTCCGCATGTAG SEQ-2
  6. ## [3] 20 CTGTGGTACGGCTCAAACGG SEQ-3
  7. ## [4] 20 CTCCCGCCTATCTCCCTTCT SEQ-4
  8. ## [5] 20 TCGCCTAGAAAAAGTTTCCT SEQ-5
  9. ## [6] 20 TCCGTATTGGAAAGCTCGTC
  10. c(DNA.set, DNAStringSet(DNA.str))
  11. ## A DNAStringSet instance of length 6
  12. ## width seq names
  13. ## [1] 20 TCCGTATTGGAAAGCTCGTC SEQ-1
  14. ## [2] 20 TTAGACCACTCCGCATGTAG SEQ-2
  15. ## [3] 20 CTGTGGTACGGCTCAAACGG SEQ-3
  16. ## [4] 20 CTCCCGCCTATCTCCCTTCT SEQ-4
  17. ## [5] 20 TCGCCTAGAAAAAGTTTCCT SEQ-5
  18. ## [6] 20 TCCGTATTGGAAAGCTCGTC

除以上列出的例子外还有duplicated, unique, sort, order, split, relist等。随着Biostrings的完善，可能会添加更多的函数，使序列的运算更符合R语言的运算习惯。

注：本文R语言的代码加亮显示由RStudio/knitr产生，使用R脚本对其产生的HTML代码进行解析后发布到本博客。

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

发表评论

匿名网友

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