R/BioC序列处理之一：Biostrings常量与序列容器

评论2,648

序列说到底就是文本/字符串类型的数据，你完全可以用纯纯的R base函数来处理，只是太麻烦，而且效率很低。BioC的IRanges包从数据结构和运算规则等角度对生物序列做了很细致的定义，是使用R高效处理生物序列数据的基础。但是IRanges包定义的类和方法（函数）实在太多了，学习起来很乏味，有空自己啃吧。我们从Biostrings包开始。

一、用R base的函数来处理序列

先写一个用于产生随机DNA/RNA/AA序列的函数：

  1. # 按长度设置产生随机序列的函数
  2. rndSeq <- function(dict, n) {
  3. paste(sample(dict, n, replace = T), collapse = "")
  4. }

这个函数的参数dict用于指定字符字典，n 是要产生的序列长度，后面有几个地方要用到它。有关R处理文本/字符串的函数我已经在另一篇文章里面介绍过了，这里只演示一下如何用这些函数来做简单的序列处理（mapply函数的用法可以参考数据汇总/透视一文）：

  1. # 产生3条长度均为50bp的序列
  2. set.seed(1000)
  3. (seqs <- mapply(rndSeq, list(c("A", "T", "G", "C")), rep(50, 3)))
  4. ## [1] "TCAGGAGGATTCTCCATGAGACTAGGTCAAGAGCGGACGTCGAGCGGGTG"
  5. ## [2] "GGAGGGGGCTTACCACTCACGGTTCAAAAATACAGATTTTCCGAGGACAT"
  6. ## [3] "GTAGGCAGGTCGAGGACAAAGTTGCCCGTATGTACTACGAAGGGGTCGTC"
  7. # 用regexpr函数查找第一次匹配ATG的位置和序列长度
  8. (pos <- regexpr("ATG", seqs))
  9. ## [1] 16 -1 30
  10. ## attr(,"match.length")
  11. ## [1] 3 -1 3
  12. ## attr(,"useBytes")
  13. ## [1] TRUE
  14. as.numeric(pos)
  15. ## [1] 16 -1 30
  16. # 用gregexpr函数可以查找匹配ATG的全部起始位置，但后续处理不是很方便
  17. pos <- gregexpr("ATG", seqs)
  18. # 获取互补序列很简单
  19. (seqs_comp <- chartr("ATGC", "TACG", seqs))
  20. ## [1] "AGTCCTCCTAAGAGGTACTCTGATCCAGTTCTCGCCTGCAGCTCGCCCAC"
  21. ## [2] "CCTCCCCCGAATGGTGAGTGCCAAGTTTTTATGTCTAAAAGGCTCCTGTA"
  22. ## [3] "CATCCGTCCAGCTCCTGTTTCAACGGGCATACATGATGCTTCCCCAGCAG"
  23. # 获取反向序列不大容易，要经过几个步骤
  24. vseq <- function(x) {
  25. substring(x, 1:nchar(x), 1:nchar(x))
  26. } #把字符串向量化的函数
  27. seqs_rc <- lapply(seqs_comp, vseq)
  28. seqs_rc <- lapply(seqs_rc, rev)
  29. (seqs_rc <- sapply(seqs_rc, paste, collapse = "")) #这一步才得到反向互补序列
  30. ## [1] "CACCCGCTCGACGTCCGCTCTTGACCTAGTCTCATGGAGAATCCTCCTGA"
  31. ## [2] "ATGTCCTCGGAAAATCTGTATTTTTGAACCGTGAGTGGTAAGCCCCCTCC"
  32. ## [3] "GACGACCCCTTCGTAGTACATACGGGCAACTTTGTCCTCGACCTGCCTAC"

二、Biostrings定义的常量

包括DNA、RNA、AA、密码子的符号，载入Biostrings包后这些常量就可以直接使用，虽然看起来简单，但很方便：

  1. library(Biostrings)
  2. DNA_BASES
  3. ## [1] "A" "C" "G" "T"
  4. DNA_ALPHABET
  5. ## [1] "A" "C" "G" "T" "M" "R" "W" "S" "Y" "K" "V" "H" "D" "B" "N" "-" "+"
  6. RNA_BASES
  7. ## [1] "A" "C" "G" "U"
  8. RNA_ALPHABET
  9. ## [1] "A" "C" "G" "U" "M" "R" "W" "S" "Y" "K" "V" "H" "D" "B" "N" "-" "+"
  10. RNA_GENETIC_CODE
  11. ## UUU UUC UUA UUG UCU UCC UCA UCG UAU UAC UAA UAG UGU UGC UGA UGG CUU CUC
  12. ## "F" "F" "L" "L" "S" "S" "S" "S" "Y" "Y" "*" "*" "C" "C" "*" "W" "L" "L"
  13. ## CUA CUG CCU CCC CCA CCG CAU CAC CAA CAG CGU CGC CGA CGG AUU AUC AUA AUG
  14. ## "L" "L" "P" "P" "P" "P" "H" "H" "Q" "Q" "R" "R" "R" "R" "I" "I" "I" "M"
  15. ## ACU ACC ACA ACG AAU AAC AAA AAG AGU AGC AGA AGG GUU GUC GUA GUG GCU GCC
  16. ## "T" "T" "T" "T" "N" "N" "K" "K" "S" "S" "R" "R" "V" "V" "V" "V" "A" "A"
  17. ## GCA GCG GAU GAC GAA GAG GGU GGC GGA GGG
  18. ## "A" "A" "D" "D" "E" "E" "G" "G" "G" "G"
  19. AA_ALPHABET
  20. ## [1] "A" "R" "N" "D" "C" "Q" "E" "G" "H" "I" "L" "K" "M" "F" "P" "S" "T"
  21. ## [18] "W" "Y" "V" "U" "B" "Z" "X" "*" "-" "+"
  22. AMINO_ACID_CODE
  23. ## A R N D C Q E G H I L K
  24. ## "Ala" "Arg" "Asn" "Asp" "Cys" "Gln" "Glu" "Gly" "His" "Ile" "Leu" "Lys"
  25. ## M F P S T W Y V U B Z X
  26. ## "Met" "Phe" "Pro" "Ser" "Thr" "Trp" "Tyr" "Val" "Sec" "Asx" "Glx" "Xaa"
  27. GENETIC_CODE
  28. ## TTT TTC TTA TTG TCT TCC TCA TCG TAT TAC TAA TAG TGT TGC TGA TGG CTT CTC
  29. ## "F" "F" "L" "L" "S" "S" "S" "S" "Y" "Y" "*" "*" "C" "C" "*" "W" "L" "L"
  30. ## CTA CTG CCT CCC CCA CCG CAT CAC CAA CAG CGT CGC CGA CGG ATT ATC ATA ATG
  31. ## "L" "L" "P" "P" "P" "P" "H" "H" "Q" "Q" "R" "R" "R" "R" "I" "I" "I" "M"
  32. ## ACT ACC ACA ACG AAT AAC AAA AAG AGT AGC AGA AGG GTT GTC GTA GTG GCT GCC
  33. ## "T" "T" "T" "T" "N" "N" "K" "K" "S" "S" "R" "R" "V" "V" "V" "V" "A" "A"
  34. ## GCA GCG GAT GAC GAA GAG GGT GGC GGA GGG
  35. ## "A" "A" "D" "D" "E" "E" "G" "G" "G" "G"
  36. IUPAC_CODE_MAP
  37. ## A C G T M R W S Y K
  38. ## "A" "C" "G" "T" "AC" "AG" "AT" "CG" "CT" "GT"
  39. ## V H D B N
  40. ## "ACG" "ACT" "AGT" "CGT" "ACGT"

还有就是用于核酸和蛋白比对的取代矩阵，这些数据需用户自行载入（最后三个是函数，但可以不设参数）：

  1. data(BLOSUM45)
  2. data(BLOSUM50)
  3. data(BLOSUM62)
  4. data(BLOSUM80)
  5. data(BLOSUM100)
  6. data(PAM30)
  7. data(PAM40)
  8. data(PAM70)
  9. data(PAM120)
  10. data(PAM250)
  11. nucleotideSubstitutionMatrix()
  12. qualitySubstitutionMatrices()
  13. errorSubstitutionMatrices()

三、Biostrings定义的类（序列容器）

1、容纳单一序列的类

有XString、BString，DNAString，RNAString和AAString共5类，其中XString是虚拟类，不能产生XString对象，其他4类都是XString的子类：

  1. getClass("XString")
  2. ## Virtual Class "XString" [package "Biostrings"]
  3. ##
  4. ## Slots:
  5. ##
  6. ## Name: shared offset length elementMetadata
  7. ## Class: SharedRaw integer integer DataTableORNULL
  8. ##
  9. ## Name: metadata
  10. ## Class: list
  11. ##
  12. ## Extends:
  13. ## Class "XRaw", directly
  14. ## Class "XVector", by class "XRaw", distance 2
  15. ## Class "Vector", by class "XRaw", distance 3
  16. ## Class "Annotated", by class "XRaw", distance 4
  17. ##
  18. ## Known Subclasses: "BString", "DNAString", "RNAString", "AAString"

再追究XString的来源就复杂了，你可以用getClass函数一个个回溯着看。

除XString外，用与类名称同名的函数就可以由字符串向量（向量长度为1）产生相应的序列对象。BString可以用任意字符（当然要符合R的字符串规则），而DNAString，RNAString和AAString对字符有相应限制。但测试结果却是AAString对字符也不检查，可能函数有问题：

  1. xxx <- "qwertyuiopasdfghjklzxcvbnm,./;'@#$%^&*()_+~"
  2. (BString(xxx))
  3. ## 43-letter "BString" instance
  4. ## seq: qwertyuiopasdfghjklzxcvbnm,./;'@#$%^&*()_+~
  5. (DNAString(rndSeq(DNA_BASES, 20)))
  6. ## 20-letter "DNAString" instance
  7. ## seq: GGGTATGGTTCGTCTTTGCC
  8. # 用法错误
  9. (DNAString(xxx))
  10. ## Error: key 113 (char 'q') not in lookup table
  11. (RNAString(rndSeq(RNA_BASES, 20)))
  12. ## 20-letter "RNAString" instance
  13. ## seq: UCUAGAUUUUCACGACGCUG
  14. # 用法错误
  15. (RNAString(xxx))
  16. ## Error: key 113 (char 'q') not in lookup table
  17. (AAString(rndSeq(AA_ALPHABET, 20)))
  18. ## 20-letter "AAString" instance
  19. ## seq: LYPBXIC+UNHSFMIBPPVS
  20. # 下面语句并不出错
  21. (AAString(xxx))
  22. ## 43-letter "AAString" instance
  23. ## seq: qwertyuiopasdfghjklzxcvbnm,./;'@#$%^&*()_+~

2、容纳序列集的类

与容纳单一序列的类对应，容纳序列集（多个序列）的类有XStringSet、BStringSet，DNAStringSet，RNAStringSet和AAStringSet（下面统称为XXXset）。XStringSet是虚拟类，不能实例化，其他4类由它派生而来。此外，由XStringSet直接派生出来还有QualityScaledXStringSet：

  1. getClass("XStringSet")
  2. ## Virtual Class "XStringSet" [package "Biostrings"]
  3. ##
  4. ## Slots:
  5. ##
  6. ## Name: pool ranges elementType elementMetadata
  7. ## Class: SharedRaw_Pool GroupedIRanges character DataTableORNULL
  8. ##
  9. ## Name: metadata
  10. ## Class: list
  11. ##
  12. ## Extends:
  13. ## Class "XRawList", directly
  14. ## Class "XVectorList", by class "XRawList", distance 2
  15. ## Class "List", by class "XRawList", distance 3
  16. ## Class "Vector", by class "XRawList", distance 4
  17. ## Class "Annotated", by class "XRawList", distance 5
  18. ##
  19. ## Known Subclasses:
  20. ## Class "BStringSet", directly
  21. ## Class "DNAStringSet", directly
  22. ## Class "RNAStringSet", directly
  23. ## Class "AAStringSet", directly
  24. ## Class "QualityScaledXStringSet", directly
  25. ## Class "XStringQuality", by class "BStringSet", distance 2
  26. ## Class "PhredQuality", by class "BStringSet", distance 3
  27. ## Class "SolexaQuality", by class "BStringSet", distance 3
  28. ## Class "IlluminaQuality", by class "BStringSet", distance 3
  29. ## Class "QualityScaledBStringSet", by class "BStringSet", distance 2
  30. ## Class "QualityScaledDNAStringSet", by class "DNAStringSet", distance 2
  31. ## Class "QualityScaledRNAStringSet", by class "RNAStringSet", distance 2
  32. ## Class "QualityScaledAAStringSet", by class "AAStringSet", distance 2

序列集类型的对象和单条序列一样用与类名称同名的函数由字符串向量产生。下面的代码用 rndSeq 和 mapply 函数产生了一个DNA序列向量，包含有4条长短不一的序列，然后用它生成了一个 DNAStringSet 对象：

  1. set.seed(1000)
  2. DNAs <- mapply(rndSeq, list(DNA_BASES), c(10000, 20000, 50000, 25000))
  3. names(DNAs) <- paste("SEQ", 1:4, sep = "-")
  4. (DNAset <- DNAStringSet(DNAs))
  5. ## A DNAStringSet instance of length 4
  6. ## width seq names
  7. ## [1] 10000 CTAGGAGGACCTCTTACGAGA...TTGCCCTAACCTCTCCCATG SEQ-1
  8. ## [2] 20000 TACATGATGGGGCCATTGGTG...TGGCTTGAAATGAGGGAATA SEQ-2
  9. ## [3] 50000 CGCGCGGTGACTTAACGTGGA...GATGACTCTACGACTATACC SEQ-3
  10. ## [4] 25000 CAATAATTCGGCCTAGGTGCG...GGCAAACCCACTCGCACCGG SEQ-4

但更多的时候是用这些容器来放置由fasta格式文件读取得到序列，这个过程由相应的readXXXSet函数完成。readXXXSet 函数在老版本的Biostrings中是read.XXXSet，现在旧版函数已被弃用。下面例子用file.choose函数选择了包含拟南芥基因组序列FASTA格式文件，有7条序列，读入后得到一个DNAStringSet对象：

  1. (chrs <- readDNAStringSet(file.choose(), "fasta"))
  2. ## A DNAStringSet instance of length 7
  3. ## width seq names
  4. ## [1] 30427671 CCCTAAACCCTAAACCCTA...TAGGGTTTAGGGTTTAGGG Chr1 CHROMOSOME
  5. ## [2] 19698289 NNNNNNNNNNNNNNNNNNN...TAGGGTTTAGGGTTTAGGG Chr2 CHROMOSOME
  6. ## [3] 23459830 NNNNNNNNNNNNNNNNNNN...AACCCTAAACCCTAAACCC Chr3 CHROMOSOME
  7. ## [4] 18585056 NNNNNNNNNNNNNNNNNNN...TTAGGGTTTAGGGTTTAGG Chr4 CHROMOSOME
  8. ## [5] 26975502 TATACCATGTACCCTCAAC...GATTTAGGGTTTTTAGATC Chr5 CHROMOSOME
  9. ## [6] 154478 ATGGGCGAACGACGGGAAT...TAACTTGGTCCCGGGCATC chloroplast CHROM...
  10. ## [7] 366924 GGATCCGTTCGAAACAGGT...AATGGAAACAAACCGGATT mitochondria CHRO...
  11. rm(chrs)

当然也有readBStringSet、readRNAStringSet 和 readAAStringSet 函数。还可以直接用构造函数由 DNAStringSet 对象直接产生 RNAStringSet 对象，但 AAStringSet 对象不行：

  1. RNAStringSet(DNAset)
  2. ## A RNAStringSet instance of length 4
  3. ## width seq names
  4. ## [1] 10000 CUAGGAGGACCUCUUACGAGA...UUGCCCUAACCUCUCCCAUG SEQ-1
  5. ## [2] 20000 UACAUGAUGGGGCCAUUGGUG...UGGCUUGAAAUGAGGGAAUA SEQ-2
  6. ## [3] 50000 CGCGCGGUGACUUAACGUGGA...GAUGACUCUACGACUAUACC SEQ-3
  7. ## [4] 25000 CAAUAAUUCGGCCUAGGUGCG...GGCAAACCCACUCGCACCGG SEQ-4

XXXSet对象构造函数同时也是XXXSet取子序列的函数，它还有几个参数，分别是start, end, width 和 use.names。前三个参数很有用，比如从染色体上截取两个位置间的序列或者截取某些位点上游或下游N个碱基的序列：

  1. # 使用start和end 参数获取子序列
  2. DNAStringSet(DNAset, start = seq(1000, length = 4, by = 1000), end = seq(1500,
  3. length = 4, by = 1000))
  4. ## A DNAStringSet instance of length 4
  5. ## width seq names
  6. ## [1] 501 GTTGCCGAACAGAGCACGGCT...AGGGAATCGTTAGCGATGAC SEQ-1
  7. ## [2] 501 AGCGGTGATGTCTGACATTGA...CTATCGCGTGGGACTAGCAC SEQ-2
  8. ## [3] 501 CCAGCGTAGTTGGGAGAATTG...ACTGCACCCGCTTCGTTGTA SEQ-3
  9. ## [4] 501 AGCACTGAGGAACAGCTGTAG...GGGAGCCGACACTAAAATTC SEQ-4
  10. # 使用 end（或start）和width参数获取子序列
  11. DNAStringSet(DNAset, end = c(234, 3000, 1029), width = 100)
  12. ## A DNAStringSet instance of length 4
  13. ## width seq names
  14. ## [1] 100 TCATGAAGGGGCTGCTCGGGT...ACTATTTCCCCGCTTGCAGG SEQ-1
  15. ## [2] 100 GGGATCATGATCGTACGCTAT...CACGCTCCAGGCCTATGAGG SEQ-2
  16. ## [3] 100 CCCATGTGGCGTTCTTATAGG...CATCAGTCAATAATCATACG SEQ-3
  17. ## [4] 100 CCACTATGGCGCGACGTAGAC...TTGTTTATCTTACTCCTCAT SEQ-4

start, end, width的组合运算使用向量循环规则，即短向量循环使用。

如果原序列是XXXSet类对象，start/end/width的向量长度不能超过原序列：

  1. class(DNAset)
  2. ## [1] "DNAStringSet"
  3. ## attr(,"package")
  4. ## [1] "Biostrings"
  5. class(DNAset[1])
  6. ## [1] "DNAStringSet"
  7. ## attr(,"package")
  8. ## [1] "Biostrings"
  9. DNAStringSet(DNAset, start = seq(1000, length = 2, by = 1000), width = c(5,
  10. 10, 15))
  11. ## A DNAStringSet instance of length 4
  12. ## width seq names
  13. ## [1] 5 GTTGC SEQ-1
  14. ## [2] 10 AGCGGTGATG SEQ-2
  15. ## [3] 15 GCTTAAGTCCCATCA SEQ-3
  16. ## [4] 5 TGTCC SEQ-4
  17. DNAStringSet(DNAset[1], start = seq(1000, length = 2, by = 1000), width = c(5,
  18. 10, 15))
  19. ## Error: 'start', 'end' or 'width' is longer than 'refwidths'

如果原序列是XString类对象，对start/end/width的向量长度没有限制，但取子序列的位置（start/end）应在原序列的范围之内：

  1. class(DNAset[[1]])
  2. ## [1] "DNAString"
  3. ## attr(,"package")
  4. ## [1] "Biostrings"
  5. DNAStringSet(DNAset[[1]], start = seq(1000, length = 2, by = 1000), width = c(5,
  6. 10, 15))
  7. ## A DNAStringSet instance of length 3
  8. ## width seq
  9. ## [1] 5 GTTGC
  10. ## [2] 10 ATAAATGGGA
  11. ## [3] 15 GTTGCCGAACAGAGC
  12. DNAStringSet(DNAset[[1]], start = seq(1000, length = 2, by = 1e+05), width = c(5,
  13. 10, 15))
  14. ## Error: solving row 2: 'allow.nonnarrowing' is FALSE and the supplied start
  15. ## (101000) is > refwidth + 1
  16. DNAStringSet(DNAset[[1]], start = seq(1000, length = 2, by = 1000), width = 1e+05)
  17. ## Error: solving row 1: 'allow.nonnarrowing' is FALSE and the solved end
  18. ## (100999) is > refwidth

组合使用3个参数中的2个基本可以满足序列截取的需要。XXXSet 对象是列表，可取元素或子集，元素是对应的String（从上面的代码已经可以看出）。

3、XStringViews类

XStringViews类由通过继承Views类（IRanges包定义）而来，XStringViews类对象用于存储同一条序列（目标序列）上的一系列“视野”，即序列区域或子序列。每一个视野由起始点（start）和终止点（end）确定，也隐含了序列的长度信息。

  1. getClass("XStringViews")
  2. ## Class "XStringViews" [package "Biostrings"]
  3. ##
  4. ## Slots:
  5. ##
  6. ## Name: subject ranges elementType elementMetadata
  7. ## Class: XString IRanges character DataTableORNULL
  8. ##
  9. ## Name: metadata
  10. ## Class: list
  11. ##
  12. ## Extends:
  13. ## Class "Views", directly
  14. ## Class "List", by class "Views", distance 2
  15. ## Class "Vector", by class "Views", distance 3
  16. ## Class "Annotated", by class "Views", distance 4
  17. ##
  18. ## Known Subclasses: "XStringPartialMatches"

XStringViews和显微镜下观察物体类似，从一大片区域内挑选感兴趣的区域存储起来。用于挑选（即获取对象）的是Views()函数（IRanges包提供），用法和XXXSet的对象构造函数差不多，但只能使用XString和长度为1的字符串向量，XXXSet类对象即使只包含有一条序列也不行：

  1. Views(DNAset[[1]], start = seq(1000, length = 2, by = 1000), width = 100)
  2. ## Views on a 10000-letter DNAString subject
  3. ## subject: CTAGGAGGACCTCTTACGAGATCAGGCTAAGA...ACGTTCACTATTTGCCCTAACCTCTCCCATG
  4. ## views:
  5. ## start end width
  6. ## [1] 1000 1099 100 [GTTGCCGAACAGAGCACGGCTCGAT...ACGACCGTTTTGTATAGGAATGAC]
  7. ## [2] 2000 2099 100 [ATAAATGGGATCCGATTCAGGTTCC...CCAAATATGAGCACCGGCGCCACT]
  8. Views(DNAs[1], start = seq(1000, length = 2, by = 1000), width = 100)
  9. ## Views on a 10000-letter BString subject
  10. ## subject: CTAGGAGGACCTCTTACGAGATCAGGCTAAGA...ACGTTCACTATTTGCCCTAACCTCTCCCATG
  11. ## views:
  12. ## start end width
  13. ## [1] 1000 1099 100 [GTTGCCGAACAGAGCACGGCTCGAT...ACGACCGTTTTGTATAGGAATGAC]
  14. ## [2] 2000 2099 100 [ATAAATGGGATCCGATTCAGGTTCC...CCAAATATGAGCACCGGCGCCACT]
  15. Views(DNAs, start = seq(1000, length = 2, by = 1000), width = 100)
  16. ## Error: zero or more than one input sequence
  17. Views(DNAset[1], start = seq(1000, length = 2, by = 1000), width = 100)
  18. ## Error: unable to find an inherited method for function 'Views' for
  19. ## signature '"DNAStringSet"'

IRanges包提供successiveViews方法同样可以用来产生XStringViews类对象，获得的是目标序列上连续的区域，它的用法是：

  1. successiveViews(subject, width, gapwidth = 0, from = 1)

使用该函数可以非常方便地获得长度和间隔都很规律的序列，例如：

  1. (successiveViews(DNAset[[1]], width = rep(20, 4), gapwidth = 0, from = 1))
  2. ## Views on a 10000-letter DNAString subject
  3. ## subject: CTAGGAGGACCTCTTACGAGATCAGGCTAAGA...ACGTTCACTATTTGCCCTAACCTCTCCCATG
  4. ## views:
  5. ## start end width
  6. ## [1] 1 20 20 [CTAGGAGGACCTCTTACGAG]
  7. ## [2] 21 40 20 [ATCAGGCTAAGAGTGGATGC]
  8. ## [3] 41 60 20 [TGAGTGGGCGGGAGGGGGTC]
  9. ## [4] 61 80 20 [CATTATCTATGGCCTAAAAA]
  10. (successiveViews(DNAset[[1]], width = rep(20, 4), gapwidth = 80, from = 1))
  11. ## Views on a 10000-letter DNAString subject
  12. ## subject: CTAGGAGGACCTCTTACGAGATCAGGCTAAGA...ACGTTCACTATTTGCCCTAACCTCTCCCATG
  13. ## views:
  14. ## start end width
  15. ## [1] 1 20 20 [CTAGGAGGACCTCTTACGAG]
  16. ## [2] 101 120 20 [GCAGGTAGGCTGAGGATAAA]
  17. ## [3] 201 220 20 [ACGCCCTGGGGGGGACTATT]
  18. ## [4] 301 320 20 [ACTAGTTTCGAGACGAGCAA]
  19. (successiveViews(DNAset[[1]], width = rep(20, 4), gapwidth = -19, from = 1))
  20. ## Views on a 10000-letter DNAString subject
  21. ## subject: CTAGGAGGACCTCTTACGAGATCAGGCTAAGA...ACGTTCACTATTTGCCCTAACCTCTCCCATG
  22. ## views:
  23. ## start end width
  24. ## [1] 1 20 20 [CTAGGAGGACCTCTTACGAG]
  25. ## [2] 2 21 20 [TAGGAGGACCTCTTACGAGA]
  26. ## [3] 3 22 20 [AGGAGGACCTCTTACGAGAT]
  27. ## [4] 4 23 20 [GGAGGACCTCTTACGAGATC]

4、MaskedXString掩膜序列类

序列的掩膜方法也是由IRanges包定义的，Biostrings包只是将该方法应用到了XString类对象上。和XString一样，MaskedXString也是虚拟类，能够实例化（产生对象实体）的是它的子类："MaskedBString", "MaskedDNAString", "MaskedRNAString", "MaskedAAString"类。

  1. getClass("MaskedXString")
  2. ## Virtual Class "MaskedXString" [package "Biostrings"]
  3. ##
  4. ## Slots:
  5. ##
  6. ## Name: unmasked masks
  7. ## Class: XString MaskCollection
  8. ##
  9. ## Known Subclasses: "MaskedBString", "MaskedDNAString", "MaskedRNAString", "MaskedAAString"

用于产生掩膜序列对象的函数是maskMotif()，如：

  1. DNAset[[1]]
  2. ## 10000-letter "DNAString" instance
  3. ## seq: CTAGGAGGACCTCTTACGAGATCAGGCTAAGAGT...GTACGTTCACTATTTGCCCTAACCTCTCCCATG
  4. (mDNA <- maskMotif(DNAset[[2]], "CATTAG"))
  5. ## 20000-letter "MaskedDNAString" instance (# for masking)
  6. ## seq: TACATGATGGGGCCATTGGTGGACAGCGTTTTTA...CTCATGGTACTCGTGGCTTGAAATGAGGGAATA
  7. ## masks:
  8. ## maskedwidth maskedratio active desc
  9. ## 1 12 6e-04 TRUE CATTAG-blocks

通过类型转换可以把MaskedXString类对象转成Views类对象，但是BStrings有更直接的方式：mask()函数。两种方式是完全等价的：

  1. as(mDNA, "Views")
  2. ## Views on a 20000-letter DNAString subject
  3. ## subject: TACATGATGGGGCCATTGGTGGACAGCGTTTT...CATGGTACTCGTGGCTTGAAATGAGGGAATA
  4. ## views:
  5. ## start end width
  6. ## [1] 1 3129 3129 [TACATGATGGGGCCATTGGTGGAC...GACCTACCACGCATCGTAATCCAT]
  7. ## [2] 3136 4979 1844 [GTTCGTCCGTTCAGAAGATATCCG...AAAGCGCGTGCCACTCTGGTGTGC]
  8. ## [3] 4986 20000 15015 [AATCGATAGAGTGTAGGGCCCAGG...CTCGTGGCTTGAAATGAGGGAATA]
  9. mask(DNAset[[2]], "CATTAG")
  10. ## Views on a 20000-letter DNAString subject
  11. ## subject: TACATGATGGGGCCATTGGTGGACAGCGTTTT...CATGGTACTCGTGGCTTGAAATGAGGGAATA
  12. ## views:
  13. ## start end width
  14. ## [1] 1 3129 3129 [TACATGATGGGGCCATTGGTGGAC...GACCTACCACGCATCGTAATCCAT]
  15. ## [2] 3136 4979 1844 [GTTCGTCCGTTCAGAAGATATCCG...AAAGCGCGTGCCACTCTGGTGTGC]
  16. ## [3] 4986 20000 15015 [AATCGATAGAGTGTAGGGCCCAGG...CTCGTGGCTTGAAATGAGGGAATA]

有时候被掩膜掉的区域才是我们敢兴趣的，比如寻找酶切位点，这时候我们要对掩膜区和非掩膜区进行反转，使用gaps()函数可以完成。而去除掩膜可以用unmasked()函数，去除掩膜对象的序列（仅保留掩膜信息）用masks()函数：

  1. as(gaps(mDNA), "Views")
  2. ## Views on a 20000-letter DNAString subject
  3. ## subject: TACATGATGGGGCCATTGGTGGACAGCGTTTT...CATGGTACTCGTGGCTTGAAATGAGGGAATA
  4. ## views:
  5. ## start end width
  6. ## [1] 3130 3135 6 [CATTAG]
  7. ## [2] 4980 4985 6 [CATTAG]
  8. unmasked(mDNA)
  9. ## 20000-letter "DNAString" instance
  10. ## seq: TACATGATGGGGCCATTGGTGGACAGCGTTTTTA...CTCATGGTACTCGTGGCTTGAAATGAGGGAATA
  11. masks(mDNA)
  12. ## MaskCollection of length 1 and width 20000
  13. ## masks:
  14. ## maskedwidth maskedratio active desc
  15. ## 1 12 6e-04 TRUE CATTAG-blocks

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

发表评论

匿名网友

拖动滑块以完成验证
加载中...