R字符串处理

字符串处理是生物信息当中经常遇到的问题。常见的操作有:赋值,获取字符串长度,截取,替换,联接,分割,比较,排序,匹配,部分匹配,正则匹配,格式化,分行。对于生物字符串,基本操作有:互补,反向,反向互补,翻译,转录,逆转录,碱基频率统计,模板搜索,搜索回文结构以及间隔成对匹配,比对。在文件操作方面,有读取及写入fasta序列,fastq序列等等。下面我就逐一介绍在R中如何操作。

常规操作

 
> #赋值
> mychar <- c("ACTACCACTAACCACT","TCATCCATTCGTGGG","GTTGTTCCATAG")
> #获取字串长度
> nchar(mychar)
[1] 16 15 12
> length(mychar)
[1] 3
> #截取
> substr(mychar, 2, 4)
[1] "CTA" "CAT" "TTG"
> substr(mychar, 2, 4)<-"TTT"
> mychar
[1] "ATTTCCACTAACCACT" "TTTTCCATTCGTGGG"  "GTTTTTCCATAG"    
> #替换
> dna2rna<-function(inputStr){
+ 	if(!is.character(inputStr))
+ 		stop("need character input")
+ 	is = toupper(inputStr)
+ 	chartr("T", "U", is)
+ }
> dna2rna(mychar)
[1] "AUUUCCACUAACCACU" "UUUUCCAUUCGUGGG"  "GUUUUUCCAUAG"    
> #联接
> paste("prefix",mychar,"postfix",sep="-")
[1] "prefix-ATTTCCACTAACCACT-postfix" "prefix-TTTTCCATTCGTGGG-postfix"  "prefix-GTTTTTCCATAG-postfix"    
> paste(mychar,collapse="--")
[1] "ATTTCCACTAACCACT--TTTTCCATTCGTGGG--GTTTTTCCATAG"
> #分割
> strsplit(mychar,"C")
[[1]]
[1] "ATTT" ""     "A"    "TAA"  ""     "A"    "T"   
 
[[2]]
[1] "TTTT"  ""      "ATT"   "GTGGG"
 
[[3]]
[1] "GTTTTT" ""       "ATAG"  
 
> strsplit(mychar,"[CG]")
[[1]]
[1] "ATTT" ""     "A"    "TAA"  ""     "A"    "T"   
 
[[2]]
[1] "TTTT" ""     "ATT"  ""     "T"    ""     ""    
 
[[3]]
[1] ""      "TTTTT" ""      "ATA"  
 
> #比较
> mychar[1]>mychar[2]
[1] FALSE
> #排序
> sort(mychar)
[1] "ATTTCCACTAACCACT" "GTTTTTCCATAG"     "TTTTCCATTCGTGGG" 
> #匹配
> exT<-c("Intron", "Exon", "promoter","enhancer")
> match("Exon", exT)
[1] 2
> "promoter" %in% exT
[1] TRUE
> #部分匹配
> pmatch("E", exT)
[1] 2
> pmatch("x", exT)
[1] NA
> charmatch("E", exT)
[1] 2
> charmatch("x", exT)
[1] NA
> pmatch(c("Exo","enh"), exT)
[1] 2 4
> #格式化
> format(1:10)
 [1] " 1" " 2" " 3" " 4" " 5" " 6" " 7" " 8" " 9" "10"
> format(1:10, trim = TRUE)
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10"
> 
> zz <- data.frame("(row names)"= c("aaaaa", "b"), check.names=FALSE)
> format(zz)
  (row names)
1       aaaaa
2           b
> format(zz, justify = "left")
  (row names)
1       aaaaa
2       b    
> 
> ## use of nsmall
> format(13.7)
[1] "13.7"
> format(13.7, nsmall = 3)
[1] "13.700"
> format(c(6.0, 13.1), digits = 2)
[1] " 6" "13"
> format(c(6.0, 13.1), digits = 2, nsmall = 1)
[1] " 6.0" "13.1"
> 
> ## use of scientific
> format(2^31-1)
[1] "2147483647"
> format(2^31-1, scientific = TRUE)
[1] "2.147484e+09"
> 
> ## a list
> z <- list(a=letters[1:3], b=(-pi+0i)^((-2:2)/2), c=c(1,10,100,1000),
+           d=c("a", "longer", "character", "string"))
> format(z, digits = 2)
                                                             a 
                                                     "a, b, c" 
                                                             b 
"-0.32+0.00i, 0.00-0.56i, 1.00+0.00i, 0.00+1.77i, -3.14+0.00i" 
                                                             c 
                                            "1, 10, 100, 1000" 
                                                             d 
                  "a        , longer   , character, string   " 
> format(z, digits = 2, justify = "left", trim = FALSE)
                                                                a 
                                                        "a, b, c" 
                                                                b 
"-0.32+0.00i,  0.00-0.56i,  1.00+0.00i,  0.00+1.77i, -3.14+0.00i" 
                                                                c 
                                         "   1,   10,  100, 1000" 
                                                                d 
                     "a        , longer   , character, string   " 
> #还有两种C风格的格式化手段,formatC和sprintf。这里就不列举了。
> #分行
> x <-paste(readLines(file.path(R.home(), "COPYING")), collapse="\n")
> strwrap(x, 30, prefix = "QIUWORLD: ")[1:6]
[1] "QIUWORLD: GNU GENERAL PUBLIC"  "QIUWORLD: LICENSE Version 2," 
[3] "QIUWORLD: June 1991"           "QIUWORLD: "                   
[5] "QIUWORLD: Copyright (C) 1989," "QIUWORLD: 1991 Free Software" 
> writeLines(strwrap(x, 30, prefix="QIUWORLD: ")[1:6])
QIUWORLD: GNU GENERAL PUBLIC
QIUWORLD: LICENSE Version 2,
QIUWORLD: June 1991
QIUWORLD: 
QIUWORLD: Copyright (C) 1989,
QIUWORLD: 1991 Free Software

正则匹配(regular expression)
正则匹配是一个非常常用的字符搜索手段,在数据挖掘中起着非常重要的作用。所以虽然它是一种常规手段,但我还是另起一段来专门讲述这个概念。
在R当中,可以使用三种正则:

  1. 扩展正则
  2. 基本正则
  3. Perl风格正则

正则的使用主要涉汲以下7个函数:grep, grepl, sub, gsub, regexpr, gregrexpr, regexec。而象strsplit, apropos以及browseEnv都是基于这7个函数基础之上的。
我们先从正则讲起。
假设我们现在需要从一堆字符当中找到一个符合一定规则的字符串,比如说从一个表格中找到所有人的email地址,或者说找到一段文字中所有的URL地址,你会如何做呢?嗯,回答肯定是正则了。正则就是做这个用的。我们知道一个email地址通常都是这样的(最简单情行),xxxxxx@ppp.ddd,其中,xxxxxx可能是任意字母,数字,以及下划线,点等组成,而ppp.ddd就是一个域名地址。它们之间以@相隔。在正则下是这样表示,
^[A-Za-z0-9\\._%+-]+@[A-Za-z0-9\\.-]+\\.[A-Za-z]{2,4}$

> pattern="^[A-Za-z0-9\\._%+-]+@[A-Za-z0-9\\.-]+\\.[A-Za-z]{2,4}$"
> str<-c("abc","someone@qiuworld.com","efg","anotherone@gmail.com","thirdone@yahoo.cn")
> #grepl会返回一个逻辑值,l就代表logical, g就代表global
> grepl(pattern,str)
[1] FALSE  TRUE FALSE  TRUE  TRUE
> #grep会返回匹配的id
> grep(pattern,str)
[1] 2 4 5
> #regexpr会返回一个数字,1表示匹配,-1表示不匹配,还会返回两个属性,匹配的长度以及是否使用useBytes。useBytes一般很少会使用到false,因为我们不处理宽字符。
> regexpr(pattern,str)
[1] -1  1 -1  1  1
attr(,"match.length")
[1] -1 20 -1 20 17
attr(,"useBytes")
[1] TRUE
 
> #regexec会返回一个list,下面的内容是第一个匹配及其长度
> regexec("\\w+@\\w+\\.[a-zA-Z]{2,4}","abc@qiuworld.com,efd@qiuworld.com")
[[1]]
[1] 1
attr(,"match.length")
[1] 16
 
> #gregexpr也会返回一个list, 下面的内容是每一个匹配及其长度以及useBytes。g就代表global
> gregexpr("\\w+@\\w+\\.[a-zA-Z]{2,4}","abc@qiuworld.com,efd@qiuworld.com")
[[1]]
[1]  1 18
attr(,"match.length")
[1] 16 16
attr(,"useBytes")
[1] TRUE
 
> #sub和gsub都用来做正则替换,其区别只在于g所代表的global。sub只替换遇到的第一个匹配,而gsub会替换所有的匹配。
> #需要注意的是,这里的匹配都是对应的一个字符串而言的,如果是多个字符串,要区别为每一个来单独对待。
> sub("\\w+@\\w+\\.[a-zA-Z]{2,4}","sub_function","abc@qiuworld.com,efd@qiuworld.com")
[1] "sub_function,efd@qiuworld.com"
> gsub("\\w+@\\w+\\.[a-zA-Z]{2,4}","gsub_function","abc@qiuworld.com,efd@qiuworld.com")
[1] "gsub_function,gsub_function"
</c>
正则的使用我们已经看到了,但是让人看不明白就是字符. \ | ( ) [ { ^ $ * + - ? 这些符号都是什么意思啊?
下面我就先来仔细讲讲perl中的正则符号
<pre lang="c">
.	# 除了换行以外的任意字符
^	# 一行字符串的起始,它并不代表第一个字符,只代表这里开始新的一行字符串。
$	# 一行字符串的结束,它并不代表最后一个字符(因为换行符并不会被包含进匹配当中),只代表一行字符串到这里结束。
*	# 零个或者多个之前的字符
+	# 一个或者多个之前的字符
?	# 零个或者一个之前的字符
# 示例部分
t.e	# t后面跟任意一个非换行字符然后跟字符e
	# 它可以匹配        the
	#                 tre
	#                 tle
	# 但是不匹配 te
	#           tale
^f	# 一行字符以f起始
^ftp	# 一行字符以ftp起始
e$	# 一行字符以e结尾
tle$	# 一行字符以tle结尾
und*	# un跟着零个或者多个d字符
	# 它会匹配 un
	#         und
	#         undd
	#         unddd (etc)
.*	# 任意一个无换行的字符串,
	#  . 匹配任何一个非换行字符
	#  * 将匹配一个或者多个之前的字符.
^$	# 一个空行.
 
# 在正则中有方括号[],代表可以匹配其中任何一个字符。而^在[]中就有了新的意义,代表“非”, -代表了“之间”
[qjk]		# q,j,k中任意一个字符
[^qjk]		# 非q,j,k的任意其它字符
[a-z]		# a至z中任意一个小写字符
[^a-z]		# 非任意一个a至z小写字符的其它字符(可以是大写字符)
[a-zA-Z]	# 任意一个英文字母
[a-z]+		# 一个或者多个小写英文字母
 
# |代表或者 小括号(...)可以把或者的部分括起来。注意小括号可能还有别的用途,但是在R当中先不使用。
 
jelly|cream	# jelly或者cream
(eg|le)gs	# eggs或者legs
(da)+		# da或者dada或者dadada 或者更多个da的重复
 
# 大括号括住1至2个数字,代表重复次数。
*	# 零个或者多个之前的字符
+	# 一个或者多个之前的字符
?	# 零个或者一个之前的字符
{n}	# n个之前的字符
{n,}	# 大于等于n个之前的字符
{n,m}	# n至m个之前的字符
 
# 下面的是一些字符被转义符\转义会赋以了一些新的(有可能是约定俗成的)意义
\n		# 换行符
\t		# tab
\w		# 任意字母(包括下划线)或者数字
		# 等同于[a-zA-Z0-9_]
\W		# \w的反义.
		# 等同于[^a-zA-Z0-9_]
\d		# 任意一个数字,等同于[0-9]
\D		# \d的反义,等同于[^0-9]
\s		# 任意一个空格,比如,
		# space, tab, newline, 等
\S		# \s的反义,任意一个非空格
\b		# 单词的边界, 不能使用在[]内
\B		# \b的反义
 
# 很明显,对于保留字符$, |, [, ), \, / 等等都需要转义字符\来转义表示:
 
\|		# 竖线
\[		# \[左方括号 \]右方括号
\)		# \(左小括号 \)右小括号
\*		# 星号
\^		# 开号
\/		# 反斜杠
\\		# 斜杠

接下来再讲一下POSIX中定义的一些特殊意义的字符

[:alnum:]	# [:alpha:]和[:digit:]的组合
[:alpha:]	# [:lower:]和[:upper:]的组合
[:blank:]	# 空格(space, tab),不包括换行符
[:cntrl:]	# 控制符,在ASCII码中000~037或者177
[:digit:]	# 任意数字:0-9
[:graph:]	# [:alnum:]和[:punct:]的组合
[:lower:]	# 当前字符集的小写字母
[:print:]	# 可打印出来的字符,[:graph:]以及空格
[:punct:]	# 标点符号,包括:^ ! " # $ % & ' ( ) * + - . / : ; < = > ? @ [ ] \ _ { } ` ~
[:space:]	# 空格,包括tab, newline, vertical tab, form feed, carriage return, and space
[:upper:]	# 当前字符集的大写字母
[:xdigit:]	# 16进制数 0-9a-fA-F

还有两个锚点特殊字符

^	# 一行字符串的起始,它并不代表第一个字符,只代表这里开始新的一行字符串。
$	# 一行字符串的结束,它并不代表最后一个字符(因为换行符并不会被包含进匹配当中),只代表一行字符串到这里结束。
\<	# 单词左边界
\>	# 单词右边界

弄清楚了这些正则符号,我们再回过头来看一点之前的

pattern <- "^[A-Za-z0-9\\._%+-]+@[A-Za-z0-9\\.-]+\\.[A-Za-z]{2,4}$"
可以改写为
pattern <- "^[\\w\\._%+-]+@[\\w\\.-]+\\.[A-Za-z]{2,4}$"
或者
pattern <- "^[[:alnum:]\\._%+-]+@[[:alnum:]\\.-]+\\.[[:alpha:]]{2,4}$"

有人会问了,为什么转义字符都要写两次啊?因为R本身也把\当成转义字符,所以在写pattern的时候,就需要使用\\来表示转义字符。还有一种办法就是设置fixed为TRUE。
那么参数中perl是什么意思呢?其实就是指是否使用PCRE的算法,我们来看实例:

> regexpr("foo|foobar","myfoobar")
[1] 3
attr(,"match.length")
[1] 6
attr(,"useBytes")
[1] TRUE
> regexpr("foo|foobar","myfoobar", perl=TRUE)
[1] 3
attr(,"match.length")
[1] 3
attr(,"useBytes")
[1] TRUE

生物字符串

必须掌握的一个包是Biostrings。它有一个基础类,BString,下面有三个继承类:DNAString, RNAString以及AAString。比对的结果会放在BStringAlign类当中。
我们还是先来看看如何做基本操作:互补,反向,反向互补,翻译,转录,逆转录,碱基频率统计。

> dna<-DNAString("TCTCCCAACCCTTGTACCAGT")
> Biostrings::dna2rna(dna)
  21-letter "RNAString" instance
seq: UCUCCCAACCCUUGUACCAGU
> rna<-transcribe(dna)
> rna
  21-letter "RNAString" instance
seq: AGAGGGUUGGGAACAUGGUCA
> rna2dna(rna)
  21-letter "DNAString" instance
seq: AGAGGGTTGGGAACATGGTCA
> cD<-cDNA(rna)
> cD
  21-letter "DNAString" instance
seq: TCTCCCAACCCTTGTACCAGT
> codons(rna)
  Views on a 21-letter RNAString subject
subject: AGAGGGUUGGGAACAUGGUCA
views:
    start end width
[1]     1   3     3 [AGA]
[2]     4   6     3 [GGG]
[3]     7   9     3 [UUG]
[4]    10  12     3 [GGA]
[5]    13  15     3 [ACA]
[6]    16  18     3 [UGG]
[7]    19  21     3 [UCA]
> AA<-translate(rna)
> AA
  7-letter "AAString" instance
seq: RGLGTWS
> complement(dna)
  21-letter "DNAString" instance
seq: AGAGGGTTGGGAACATGGTCA
> reverseComplement(dna)
  21-letter "DNAString" instance
seq: ACTGGTACAAGGGTTGGGAGA
> reverse(dna)
  21-letter "DNAString" instance
seq: TGACCATGTTCCCAACCCTCT
 
> library("BSgenome.Hsapiens.UCSC.hg19")
Loading required package: BSgenome
Loading required package: GenomicRanges
> Hsapiens
Human genome
| 
| organism: Homo sapiens (Human)
| provider: UCSC
| provider version: hg19
| release date: Feb. 2009
| release name: Genome Reference Consortium GRCh37
| 
| single sequences (see '?seqnames'):
|   chr1                   chr2                   chr3                 
|   chr4                   chr5                   chr6                 
|   chr7                   chr8                   chr9                 
|   chr10                  chr11                  chr12                
|   chr13                  chr14                  chr15                
|   chr16                  chr17                  chr18                
|   chr19                  chr20                  chr21                
|   chr22                  chrX                   chrY                 
|   chrM                   chr1_gl000191_random   chr1_gl000192_random 
|   chr4_ctg9_hap1         chr4_gl000193_random   chr4_gl000194_random 
|   chr6_apd_hap1          chr6_cox_hap2          chr6_dbb_hap3        
|   chr6_mann_hap4         chr6_mcf_hap5          chr6_qbl_hap6        
|   chr6_ssto_hap7         chr7_gl000195_random   chr8_gl000196_random 
|   chr8_gl000197_random   chr9_gl000198_random   chr9_gl000199_random 
|   chr9_gl000200_random   chr9_gl000201_random   chr11_gl000202_random
|   chr17_ctg5_hap1        chr17_gl000203_random  chr17_gl000204_random
|   chr17_gl000205_random  chr17_gl000206_random  chr18_gl000207_random
|   chr19_gl000208_random  chr19_gl000209_random  chr21_gl000210_random
|   chrUn_gl000211         chrUn_gl000212         chrUn_gl000213       
|   chrUn_gl000214         chrUn_gl000215         chrUn_gl000216       
|   chrUn_gl000217         chrUn_gl000218         chrUn_gl000219       
|   chrUn_gl000220         chrUn_gl000221         chrUn_gl000222       
|   chrUn_gl000223         chrUn_gl000224         chrUn_gl000225       
|   chrUn_gl000226         chrUn_gl000227         chrUn_gl000228       
|   chrUn_gl000229         chrUn_gl000230         chrUn_gl000231       
|   chrUn_gl000232         chrUn_gl000233         chrUn_gl000234       
|   chrUn_gl000235         chrUn_gl000236         chrUn_gl000237       
|   chrUn_gl000238         chrUn_gl000239         chrUn_gl000240       
|   chrUn_gl000241         chrUn_gl000242         chrUn_gl000243       
|   chrUn_gl000244         chrUn_gl000245         chrUn_gl000246       
|   chrUn_gl000247         chrUn_gl000248         chrUn_gl000249       
| 
| multiple sequences (see '?mseqnames'):
|   upstream1000  upstream2000  upstream5000  
| 
| (use the '$' or '[[' operator to access a given sequence)
> chr22NoN<-mask(Hsapiens$chr22,"N")
> alphabetFrequency(Hsapiens$chr22, baseOnly=TRUE)
      A       C       G       T   other 
9094775 8375984 8369235 9054551       0 
> alphabetFrequency(Hsapiens$chr22)
      A       C       G       T       M       R       W       S 
9094775 8375984 8369235 9054551       0       0       0       0 
      Y       K       V       H       D       B       N       - 
      0       0       0       0       0       0       0       0 
      + 
      0 
> hasOnlyBaseLetters(Hsapiens$chr22)
[1] TRUE
> uniqueLetters(Hsapiens$chr22)
[1] "A" "C" "G" "T"
> GC_content <- letterFrequency(Hsapiens$chr22, letters="CG")
> GC_content
     C|G 
16745219 
> GC_content <- letterFrequency(Hsapiens$chr22, letters="CG")/letterFrequency(Hsapiens$chr22, letters="ACGT")
> GC_content
      C|G 
0.4798807

下面来看看如何操作模板搜索,搜索回文结构以及间隔成对匹配,比对。对于模板搜索,可以使用matchPattern

> TATA = "TATAAAA"
> mT = matchPattern(TATA, chr22NoN)
> countPattern(TATA, chr22NoN)
[1] 5276
> mmT = matchPattern(TATA, chr22NoN, max.mismatch = 1)
> length(mmT)
[1] 102178
> mismatch(TATA, mmT[1:5])
[[1]]
[1] 2
 
[[2]]
[1] 5
 
[[3]]
[1] 7
 
[[4]]
[1] 1
 
[[5]]
[1] 2
 
> #批量搜索。需要注意的是,现在在批量搜索中,做为模板的长度必须是一样的,不能不一样长。
> library(hgu133a2probe)
> dict<-hgu133a2probe$sequence
> length(dict)
[1] 247899
> unique(nchar(dict))
[1] 25
> dict[1:3]
[1] "CACCCAGCTGGTCCTGTGGATGGGA" "GCCCCACTGGACAACACTGATTCCT"
[3] "TGGACCCCACTGGCTGAGAATCTGG"
> pdict<-PDict(dict, max.mismatch=1)
> vindex <- matchPDict(pdict, Hsapiens$chr22)
WARNING: using 'max.mismatch=0' with an MTB_PDict object that was preprocessed for allowing up to 1 mismatching letter(s) per match is not optimal
> length(vindex)
[1] 247899
> count_index<-countIndex(vindex)
> sum(count_index)
[1] 51471
> table(count_index)
count_index
     0      1      2      3      4      5      6      7      8      9 
243903   3549    187    100     23     25      5      4      7      4 
    10     11     12     13     14     15     17     19     27     30 
     5      4      2      2      3      1      1      1      2      1 
    31     32     33     35     68     90    147    153    179    186 
     1      1      1      1      1      1      1      1      1      1 
   190    194    196    197    205    215    250    265    275    284 
     1      1      1      1      1      1      1      1      1      1 
   289    297    309    310    324    331    336    338    365    384 
     1      1      1      1      1      1      1      1      2      1 
   413    417    421    444    454    461    468    479    480    487 
     1      1      1      1      1      2      1      1      1      1 
   493    503    515    789    823    886    904    921    932    953 
     1      1      1      1      1      1      1      1      1      1 
   973   1146   1147   1206   1227   1270   1271   1284   1299   1301 
     1      1      1      1      1      1      1      1      1      1 
  1304   1306   1311   1317   1958   2130   2760   2773 
     1      1      1      1      1      1      1      1 
> dict[count_index == max(count_index)]
[1] "CTGTAATCCCAGCACTTTGGGAGGC"
> countPattern("CTGTAATCCCAGCACTTTGGGAGGC",Hsapiens$chr22)
[1] 2773
 
> #搜索回文结构,以及成对匹配。
> chr22_pals <- findPalindromes(chr22NoN, min.armlength=40, max.looplength=20)
> nchar(chr22_pals)
 [1]  83  96 107  94  81  90  88  91  88 136  91  88 106 100 100  88
[17]  97  81  82  81  85  89  93  97 101 105 109 111 107 103  99  95
[33]  91  87  83  81  86  83  85  91  83  89 113  83  96  98  97 127
[49]  95  80  85  88  83 100  97  94  87  83 105 125 104  81  83  93
> palindromeArmLength(chr22_pals)
 [1]  83  96 107  94  81  43  88  40  40  64  40  40 106  43  43  40
[17]  40  81  82  81  85  89  93  97 101 105 109 111 107 103  99  95
[33]  91  87  83  81  86  83  85  91  83  89  47  83  96  98  97 127
[49]  95  80  85  88  83  40  40  41  41  83 105  53 104  81  83  93
> head(palindromeLeftArm(chr22_pals))
  Views on a 51304566-letter DNAString subject
subject: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNNNNNNNNNNN
views:
       start      end width
[1] 16288595 16288677    83 [GCGCGTGCGGCGTGCGCGT...TGCGCGTGCGGCGTGCGCG]
[2] 16288595 16288690    96 [GCGCGTGCGGCGTGCGCGT...TGCGCGTGCGGCGTGCGCG]
[3] 16288596 16288702   107 [CGCGTGCGGCGTGCGCGTG...GTGCGCGTGCGGCGTGCGC]
[4] 16288609 16288702    94 [CGCGTGCGGCGTGCGCGTG...GTGCGCGTGCGGCGTGCGC]
[5] 16288622 16288702    81 [CGCGTGCGGCGTGCGCGTG...GTGCGCGTGCGGCGTGCGC]
[6] 18778213 18778255    43 [ATAGATAGATAGATAGATA...ATAGATAGATAGATAGATA]
> ans <- alphabetFrequency(chr22_pals, baseOnly=TRUE)
> head(ans)
      A  C  G  T other
[1,]  0 26 45 12     0
[2,]  0 30 52 14     0
[3,]  0 34 57 16     0
[4,]  0 30 50 14     0
[5,]  0 26 43 12     0
[6,] 46  0 21 23     0
> Lpattern<-"CTCCGAG"
> Rpattern<-"GTTCACA"
> LRans<-matchLRPatterns(Lpattern, Rpattern, 500, Hsapiens$chr22)
> length(LRans)
[1] 21
> head(LRans)
  Views on a 51304566-letter DNAString subject
subject: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNNNNNNNNNNN
views:
       start      end width
[1] 16602263 16602567   305 [CTCCGAGGGTTTGAATGAT...TGCTGCTGCTCGGTTCACA]
[2] 17620584 17620854   271 [CTCCGAGTAGCTGGGATTA...CCTTGCTTTTCAGTTCACA]
[3] 21741812 21742230   419 [CTCCGAGTGCCAAGGCCAG...AGATGCCCAAGGGTTCACA]
[4] 23793654 23793992   339 [CTCCGAGTTAACCGCCATT...TCTACAAGCGAGGTTCACA]
[5] 33218204 33218260    57 [CTCCGAGGCTGGCTCCTGT...GTCTCTGGGTATGTTCACA]
[6] 33692643 33693134   492 [CTCCGAGGTTCAAGCGATT...CTACCTGGCTTTGTTCACA]
 
> #比对
> aa1<-AAString("HXBLVYMGCHFDCXVBEHIKQZ")
> aa2<-AAString("QRNYMYCFQCISGNEYKQN")
> needwunsQS(aa1, aa2, "BLOSUM62", gappen=3)
Global PairwiseAlignments (1 of 1)
pattern: [1] HXBLVYMGCHFDCXV-BEHIKQZ 
subject: [1] QRN--YMYC-FQCISGNEY-KQN 
score: 39 
> needwunsQS(aa1, aa2, "BLOSUM62", gappen=8)
Global PairwiseAlignments (1 of 1)
pattern: [1] HXBLVYMGCHFDCXVBEHIKQZ 
subject: [1] QRN--YMYC-FQCISGNEYKQN 
score: 17 
> dna1<-DNAString("CTCCGAGGGTTTGAATGAT")
> dna2<-DNAString("CTCCGAGTAGCTGGGATTA")
> mat <- matrix(-5L, nrow = 4, ncol = 4)
> for (i in seq_len(4)) mat[i, i] <- 0L
> rownames(mat) <- colnames(mat) <- DNA_ALPHABET[1:4]
> needwunsQS(dna1, dna2, mat, gappen = 0)
Global PairwiseAlignments (1 of 1)
pattern:  [5] GAGGGTTTGAATGA 
subject: [15] ----GA-----T---T-A 
score: 0 
 
> needwunsQS(dna1, dna2, mat, gappen = 1)
Global PairwiseAlignments (1 of 1)
pattern: [1] CTCC--G-AG--GGTTTGA-ATGA 
subject: [1] CTCCGAGTAGCTGG---GAT-T-A 
score: -12 
 
> needwunsQS(dna1, dna2, mat, gappen = 5)
Global PairwiseAlignments (1 of 1)
pattern: [1] CTCCGAG-GGTTTGAATGA 
subject: [1] CTCCGAGTAGCTGGGATTA 
score: -35

文件操作

读写fastq, fasta文件应该是一项非常常规的操作了。可以使用readDNAStringSet, writeXStringSet

## ---------------------------------------------------------------------
## A. READ/WRITE FASTA FILES
## ---------------------------------------------------------------------
filepath <- system.file("extdata", "someORF.fa", package="Biostrings")
fasta.info(filepath, seqtype="DNA")
x <- readDNAStringSet(filepath)
x
out1 <- tempfile()
writeXStringSet(x, out1)
 
## ---------------------------------------------------------------------
## B. READ/WRITE FASTQ FILES
## ---------------------------------------------------------------------
filepath <- system.file("extdata", "s_1_sequence.txt",
                        package="Biostrings")
fastq.geometry(filepath)
readDNAStringSet(filepath, format="fastq")
 
library(BSgenome.Celegans.UCSC.ce2)
## Create a "sliding window" on chr I:
sw_start <- seq.int(1, length(Celegans$chrI)-50, by=50)
sw <- Views(Celegans$chrI, start=sw_start, width=10)
my_fake_shortreads <- as(sw, "XStringSet")
my_fake_ids <- sprintf("ID%06d",  seq_len(length(my_fake_shortreads)))
names(my_fake_shortreads) <- my_fake_ids
my_fake_shortreads
 
## Fake quality ';' will be assigned to each base in 'my_fake_shortreads':
out2 <- tempfile()
writeXStringSet(my_fake_shortreads, out2, format="fastq")
 
## Passing qualities thru the 'qualities' argument:
my_fake_quals <- rep.int(BStringSet("DCBA@?>=<;"),
                         length(my_fake_shortreads))
my_fake_quals
out3 <- tempfile()
writeXStringSet(my_fake_shortreads, out3, format="fastq",
                qualities=my_fake_quals)
 
## ---------------------------------------------------------------------
## C. SERIALIZATION
## ---------------------------------------------------------------------
saveXStringSet(my_fake_shortreads, "my_fake_shortreads", dirpath=tempdir())

6 thoughts on “R字符串处理

  1. 求教!read.table的时候特殊字符被转换了,比如数字开头的字符串前面被加个x,逗号被转换成点,等等,有那个参数可以设置下不转换这些特殊字符不?

    1. 30是指的字符的宽度。参考:
      x < - paste(letters, collapse=" ") strwrap(x, 10) strwrap(x, 20) strwrap(x, 46) nchar(strwrap(x, 46)[1])
      注意,必须要有空格,否则无法分词。

发表评论

电子邮件地址不会被公开。 必填项已用*标注