R语言的字符串操作

R语言主要擅长于数值向量和矩阵操作,但是让他去做字符串操作也可以。

字符串的基本操作类型:

  • 查找和替换
  • 大小写转换
  • 字符数统计
  • 字符串连接和拆分

就我所知,有两套处理函数,一套是Hadley大神的stringr,一套是R自带的。

stringr使用指南

stringr函数主要分为四类:

1. 字符操作:操作字符向量中的单个符 str_length, str_sub, str_dup

2. 添加,移除和操作空白符 str_pad, str_trim, str_wrap

3. 大小写转换处理 str_to_lower, str_to_upper, str_to_title

4. 模式匹配函数 str_detect, str_subset, str_count, str_locate, str_locate_all, str_match, str_match_all, str_replace, str_replace_all, str_split_fix, str_split, str_extract, str_extract_all 

单个字符的处理

字符长度str_length,等价于nchar

1
str_length("abc")

根据位置信息提取或替换字符, 类似于substr()

1
2
3
4
5
6
7
x <- c("abcdef","ghijk")
# 第三个
str_sub(x,3,3)
# 第二个到倒数第二个
str_sub(x,2, -2)
# 替换
str_sub(x,3,3) <- X

重复字符串,不同于rep

1
str_dup(x,c(2,3))

空白符

前后增加空白字符, str_pad()

1
2
3
x <- c("abc", "defghi")
str_pad(x, 10)
str_pad(x, 10, "both")

移除空白字符, str_trim()

1
2
3
x <- c("   b   ", "c   ", "   d")
str_trim(x)
str_trim(x,left)

更好的排版,让每一行的看起来一样长, str_wrap()

1
2
3
4
nature <- c("Nature Methods' Points of Significance column ",
"on statistics explains many key statistical and ","experimental design concepts. Other resources include an online",
" plotting tool and links to statistics guides from other publishers.")
cat(str_wrap(nature, width=40))

大小写转换

  • 全部大写 str_to_upper, 类似于基础R的toupper()
  • 全部小写 str_to_lower, 类似于基础R的tolower()
  • title形式 str_to_title

模式匹配

功能: decect, locate, extract, match, replace, split

测试数据:

1
2
3
4
5
6
7
strings <- c(
    "apple",
    "219 733 8965",
    "329-293-8753",
    "Work: 579-499-7527; Home: 543.355.3679"
)

号码的正则形式:

1
phone <- "([2-9][0-9]{2})[- .]([0-9]{3})[- .]([0-9]{4})"

检测字符串是否符合特定模式, str_detect

1
str_detect(strings, phone)

提取字符串(全部内容), str_subset

1
str_subset(strings, phone)

统计匹配次数,str_count

1
str_count(strings, phone)

定位首个匹配的位置,str_locate(返回matrix), 或所有匹配的位置, str_locate_all(返回list)

1
2
str_locate(strings, phone)
str_locate_all(strings, phone)

str_extract提取首个匹配,返回字符串向量, str_extract_all()提取所有匹配,返回list

1
2
3
str_extract(strings, phone)
str_extract_all(strings,phone)
str_extract_all(strings, phone, simplify=TRUE)

str_match提取首个匹配中()分组内的内容, str_match_all()则是全部,因此是list

1
2
str_match(strings, phone)
str_match_all(strings, phone)

str_replace()替代第一个匹配, str_replace_all替代所有匹配

1
2
str_replace(strings, phone, "XXX-XXX-XXXX")
str_replace_all(strings, phone, "XXX-XXX-XXXX")

str_split_fixed()返回固定数目,全部拆分,str_split()可变拆分

1
2
str_split_fixed("a-b-c", "-")
str_spilt("a-b-c","-", n=2)

4类字符串描述引擎

1. 默认是正则表达式`vignette(“regular-expression”)

2. 逐byte固定匹配, fixed()

3. locale-sensive 字符匹配, coll()

4. 字符边界分析, boundary()

正则表达式练习

基本匹配

最简单的模式就是匹配某个完整的字符

1
2
x <- c("apple", "banana", "pear")
str_extract(x, "an")

如果需要忽略大小写, ignore_case =TRUE

1
2
bananas <- c("banana", "Banana", "BANANA")
str_dectect(bananas, regrex("banana", ignore_case=TRUE))

可以用点.匹配任意字符,但是默认不包括\n,需要用dotall=TURE开启

1
2
str_detect("\nX\n", ".X.")
str_detect("\nX\n", regex(".X.", dotall=TRUE))

转义(R中一坑)

正则表达式中有一些是特殊字符,比如说刚才的顿号,因此为了匹配这些特殊字符,我们需要对其转义。在Linux命令行里,转义用的是\, 所以可以直接用\..

但是R里面的坑就出现了,我们用字符表示正则表达式,\ 在字符里被用作转义符号。然后我们需要先把\转义成字符,然后才能进一步转义,\\.

如果要匹配\ ,就需要\\\\,不可以思议,难以释怀,不知道被坑了多少次。

或者你用\Q...\E 类似于Python的 r’….’, 原意匹配

特殊字符

一些比较常用的字符匹配

  • \d: 任意数字, \D:任意非数字.
  • \s: 任意空白字符,\S:任意非空白字符
  • \w: 匹配单词
  • \b: 匹配字符边界, \B:非字符边界
  • [abc], [a-z], [^abc], [^-]:匹配字母,和不匹配

在R里面,需要对”"进行转义,所以上面的\在R里都要写成,\
下面是一些预编译好的字符集,顾名思义

  • [:punct:]
  • [:alpha:]
  • [:lower:]
  • [:upper:]
  • [:digit:]
  • [:xdigit:]
  • [:alnum:]
  • [:cntrl:]
  • [:graph:]
  • [:print:]
  • [:space:]
  • [:blank:]

匹配abc或def

1
str_detect(c("abc","def","ghi"), "abc|def")

分组

匹配grey或gray

1
str_extract(c("grey","gray"), "gr(e|a)y")

分组可以用\1, \2进行提取,

定位

  • ^: 字符串开始, 如^a
  • $: 字符串结束, 如a$

如果字符串有多行,那么就需要regex(multiline=TRUE)。此时,

  • \A: 输入开头
  • \z: 输入结尾
  • \Z: 头尾

重复

  • ?: 0或1
  • +: 大于等于1
  • *: 大于等于0
  • {n}: n次
  • {n,m}: n到m次
  • {n,}: 大于那次

默认是贪婪模式, 在上述字符后添加”?” 则为非贪婪模式。

PS: 下面是R语言自带的字符处理函数,我已经放弃他们了。

基础R包函数

nchar(): 函数返回字符串长度
paste()paste0(): 连接若干个字符串
sprintf():格式化输出,下面举例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
sprintf("%f", pi)
sprintf("%.3f", pi)
sprintf("%1.0f", pi)
sprintf("%5.1f", pi)
sprintf("%05.1f", pi)
sprintf("%+f", pi)
sprintf("% f", pi)
sprintf("%-10f", pi) # left justified
sprintf("%e", pi)
sprintf("%E", pi)
sprintf("%g", pi)
sprintf("%g",   1e6 * pi) # -> exponential
sprintf("%.9g", 1e6 * pi) # -> "fixed"
sprintf("%G", 1e-6 * pi)

toupper(): 大写转换
tolower(): 小写转换
substr(): 提取或替换一个字符串向量的子串

1
2
3
4
5
x <- "abcde"
substr(x,1,2)
# ab
substr(x,1,2) <- 2333
# 233cde

上面都是一些普通的函数,很好理解,下面都是一些和正则表达式相关的函数,如grep, grepl, regexpr, gregexpr, sub, gsub, strsplit
因此必须介绍一下R语言的正则表达式写法了。

  • R语言是用的扩展正则表达式(Extended Regular Expressions)
  • 元字符:\ | ( ) [ { ^ $ * + ?
  • 非元字符转义后:\a as BEL, \e as ESC, \f as FF, \n as LF, \r as CR and \t as TAB
  • 一些定义字符集合[:alnum:], [:alpha:], [:blank:], [:cntrl:], [:digit:], [:graph:], [:lower:], [:print:], [:punct:], [:space:], [:upper:],[:xdigit:]
  • 找出“组”字符串 
  • 默认是贪婪模式,可以通过用?改变为非贪婪模式

这些是基本知识,可以百度到每个字符的具体解释,或者看文档?regexp
不说基础知识了,看下应用吧。我常用的操作一般是找到某个字符串,或者对字符串进行替换

比如说,我想找到所有以P开头,且不是P结尾的字符,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
test <- c("Python", "Perl", "PHP", "JAVA", "C", "C++")
grep("^P.*?[^P]$", test)
[1] 1 2
grep("^P.*?[^P]$", test,value=TRUE)
[1] "Python" "Perl"
grepl("^P.*?[^P]$", test)
[1]  TRUE  TRUE FALSE FALSE FALSE FALSE
regexpr("^P.*?[^P]$", test)
[1]  1  1 -1 -1 -1 -1
attr(,"match.length")
[1]  6  4 -1 -1 -1 -1
attr(,"useBytes")
> gregexpr("^P.*?[^P]$", test)
[[1]]
[1] 1
attr(,"match.length")
[1] 6
attr(,"useBytes")
[1] TRUE

[[2]]
[1] 1
attr(,"match.length")
[1] 4
attr(,"useBytes")
[1] TRUE

其中grep()默认是返回下标,如果设置value=TRUE,则返回字符串,grepl()返回是否配对的逻辑判断, regexpr则是返回匹配范围,如果不匹配结果是-1,gregexpr和前者功能一致,只不过返回的是列表形式。
:忽略大小写ignore.case = TRUE

现在我想把C++替换成C--。我先试着找到C++

1
2
> grep("\+\+",test)
错误: 由""\+"开头的字符串中存在'\+',但没有这种逸出号

什么情况,为什么\+不能把+这个元字符转义?难不成+在R里面不是元字符?我测试下

1
2
3
grep("++",test,value=TRUE)
Error in grep("++", test, value = TRUE) : 
  正规表现’++'不对,原因是'Invalid use of repetition operators'

啊!看来+还是元字符,难道是\ 叛变革命了,我试试看。

1
2
3
4
5
> grep("\23","test\23",value=TRUE)
[1] "test\023"
> grep("\\23","test\23",value=TRUE)
Error in grep("\\23", "test\023", value = TRUE) : 
  正规表现’\23'不对,原因是'Invalid back reference'

看来\ 是主要任务是把非元字符转义,如果想把元字符转义成普通字符,只能是\\元字符

1
2
grep("\\+\\+",test,value=TRUE)
[1] "C++"

回到我们之前的替换任务sub只对第一个匹配进行替换,gsub对所有匹配替换。

1
2
 sub("\\+\\+","--",test)
[1] "Python" "Perl"   "PHP"    "JAVA"   "C"      "C--" 

最后还可以用strsplit对字符串进行分割,返回的是一个列表

1
2
x <- c(as = "asfef", qu = "qwerty", "yuiop[", "b", "stuff.blah.yech")
strsplit(x, "e")

从零开始学Circos绘制圈图(一)

一般基因组文章都会有下面这种酷炫图,用来描述基因组的基因密度分布,转座子的密度分布,和其他物种或者多倍体的多套染色体间的共线性关系,以及其他各种你只要测序就能加上的信息,比如说你要是测了ATAC-seq,加上全基因组开放状态,要是测了多个组织,多个时期的RNA-seq,那就加上热图展现这种变化关系。

circos绘制基因组

当然除了基因组文章,其他类型的文章也可以考虑这种图。接下来我将会写一些列教程(可能有视频),通过教别人学Circos的方式来自学Circos。

环境配置

建议在Linux环境下配置Circos,之后只要用conda就能配置好分析环境

1
2
3
# 安装
## circos
conda create -c bioconda -n circos circos

测试软件安装结果

1
2
3
4
5
6
# 测试circos
conda activate circos
# 确认安装
circos -V
# 显示如下
# circos | v 0.69-8 | 15 Jun 2019 | Perl 5.026002

可以从http://circos.ca/software/download/下载官方的教程文件,分别是

处理过程

Circos依赖于一些列的配置文件,用来定义复杂图形的各个部分,最终加工成图形。

因此,用Circos画图是一个不断增添内容的过程,你要不断根据输出结果来调整输入参数。

并且整个分析中,你还要拥有过关的数据预处理的能力,这是因为Circos不是数据处理工具,它只是展示你已有的数据。

处理过程

快速开始

不管怎么样,先快速绘制出一个Circos图再说。

果子老师说过,我们不是先成为了老司机才开车,而是开车多了才成为了老司机。

第一步,先新建一个文件夹,用于存放本次分析的所有数据和配置文件

1
mkdir -p my_first_circos && cd my_first_circos

然后用vim karyotype.tair10.txt编辑文本,新增如下内容

1
2
3
4
5
chr - chr1 chr1 0 30427617 black
chr - chr2 chr2 0 19698289 black
chr - chr3 chr3 0 23459830 black
chr - chr4 chr4 0 18585056 black
chr - chr5 chr5 0 26975502 black

之后创建一个circos.conf文件,用于增加各类配置参数

1
touch circos.conf

vim circos.conf,增加我们的第一条记录,染色体信息

1
karyotype = karyotype.tair10.txt

然而要想真正的出图,还需要增加至少以下配置语句才行

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
<ideogram>
<spacing>
default = 0.005r
</spacing>
radius = 0.90r
thickness = 20p
fill = yes
stroke_color = dgrey
stroke_thickness = 2p
</ideogram>

<image>
<<include etc/image.conf>>
</image>

<<include etc/colors_fonts_patterns.conf>>
<<include etc/housekeeping.conf>>

在当前路径下运行circos -conf circos.conf, 最终效果图如下

第一张图

虽然图比较丑,但是至少我们成功运行了人生第一次的circos, 这就相当于买了一套毛坯房,后面要做的事情就是不断装修。

比如说,我们至少可以让不同染色体拥有不同的颜色,修改之前的karyotype.tair10.txt中的最后一列

1
2
3
4
5
chr - chr1 chr1 0 30427617 chr1
chr - chr2 chr2 0 19698289 chr2
chr - chr3 chr3 0 23459830 chr3
chr - chr4 chr4 0 18585056 chr4
chr - chr5 chr5 0 26975502 chr5

在当前路径下运行circos -conf circos.conf, 最终效果图如下

1564117790558

这就引出了第一个知识点,配色

为了实现配色,需要circos.conf文件了有一个和配色有关的语句

1
<<include etc/colors_fonts_patterns.conf>>

这里<<>>表示通过相对路径的方式加载另外一个配置文件,它的实际路径是和circos所在目录同级的etc,可用下面语句看到colors_fonts_patterns.conf的内容

1
2
circos_path=$(dirname `which circos`)
less ${circos_path%bin}/etc/colors_fonts_patterns.conf

你会发现,这个文件里还嵌套其他的配置文件。最终通过层层排查,你才知道etc/colors.ucsc.conf才是实际定义我们填写的颜色名的文件,而颜色的定义如下:

1
2
3
4
5
chr1  = 153,102,0
chr2 = 102,102,0
chr3 = 153,153,30
chr4 = 204,0,0
chr5 = 255,0,0

还有一个问题,为什么这里用的是两个尖括号<<,而不是一个尖括号<呢?这是因为<已经被用于分隔不同的语句块,如下语句就表示etc/image.conf里的配置信息是用来调整和image有关的配置,而不是去调整ideogram的配置。

1
2
3
<image>
<<include etc/image.conf>>
</image>

以上是快速开始部分,后续将会在此基础上,做出发表级别的图。

如何用软件模拟NGS数据

如何用软件模拟NGS数据

为了评价一个工具的性能,通常我们都需要先模拟一批数据。这样相当于有了参考答案,才能检查工具的实际表现情况。因此对于我们而言,面对一个新的功能,可以先用模拟的数据测试下不同工具的优缺点。有如下几个工具值得推荐一下:

  • ‘wgsim/dwgsim’: 从全基因组中获取测序reads
  • ‘msbar’: EMBOSS其中一个工具,能够从单个序列中模拟随机突变
  • ‘biosed’: EMBOSS的一个工具,可以按照我们给定突变位点模拟
  • ‘ReadSim’: 专门用于模拟PacBio/Nanopore这类仪器产生的long read
  • ‘Art’: 目前最复杂的模拟工具,能够模拟测序仪测序引入的错误位点
  • ‘Metasim’: 用于模拟宏基因组得到的reads
  • ‘Polyester’: 用于模拟RNA-seq

值得注意的是,这些工具模拟效果是有限,比如建库操作中超声破碎会出现的误差就很难模拟。但是最好的用途就是看看不同生物学事件在数据的情况,比如说发生了“大规模倒置”的基因组得到的数据比对到参考基因组上会是什么情况。

使用dwgsim进行模拟

wgismdwgsim能够根据参考基因组模拟出测序reads,主要是二倍体基因组的SNPs和插入缺失(INDEL)多态位点。wgism容易安装,但是参考答案是以简单的文本格式保存,不容易可视化。dwgsimwgism启发,虽然安装稍微麻烦了点,但是参考答案是以VCF格式保存,很方便可视化。

1
2
3
4
5
6
7
8
9
10
# 请先安装好ncurse
# 安装dwgsim
mkdir -p ~/scr
mkdir -p ~/.local/bin
cd ~/src
git clone --recursive https://github.com/nh13/DWGSIM.git
cd DWGSIM
make
ln -s ~/src/DWGSIM/dwgsim ~/.local/bin/dwgsim
ln -s ~/src/DWGSIM/dwgsim_eval ~/.local/bin/dwgsim_eval

简单地模拟一批数据

1
2
3
4
5
6
7
# efetch 需要用到conda安装启动
# conda create -n entrez entrez-direct
# conda activate entrez
# 获取参考基因组
efetch -db=nuccore -format=fasta -id=AF086833 > genome.fa
# 模拟数据
~/.local/bin/dwgsim genome.fa data

会得到如下数据

1
2
3
4
5
|-- data.bfast.fastq.gz # 用于bfast
|-- data.bwa.read1.fastq.gz # 用于BWA的R1
|-- data.bwa.read2.fastq.gz # 用于BWA的R2
|-- data.mutations.txt
|-- data.mutations.vcf # VCF形式擦

随后将这批数据用BWA比对,以bcftools检测变异和参考答案比较一下。

1
2
3
4
5
# conda install bwa samtools bcftools
bwa index genome.fa
bwa mem genome.fa data.bwa.read1.fastq.gz data.bwa.read2.fastq.gz | samtools sort -o data.bwa.bam
samtools mpileup -uf genome.fa data.bwa.bam | bcftools call -mv -o data.bwa.vcf
samtools index data.bwa.bam

利用使用IGV可视化,检查分析结果和真集的一致性

IGV检查

说明samtools+bcftools找变异这个组合其实还是靠谱的,至少在动植物领域研究里应该够用。

通过grep来学习正则表达式

正则表达式(regular expression, regex)是一个重要且实用的概念,我时常提起却从未细谈。一怕能力不够说不清楚反而会误导人,二是已经有无数前人撰文介绍。考虑日常用到的grep,sed,awk里经常需要用到正则表达式,于是开一个小系列,介绍如何在grepsed,awk中适用正则。

什么是正则表达式

正则表达式(regular expression)的概念,最初来自于20世纪40年代的两位神经学家(Warren McCulloch, Walter Pitts)研究神经元时提出的想法。后来数学家Stephen Kleene在代数学中正式描述了这种被他称之为“正则集合”的模型。并且,他还发明了一套简洁的方法表示正则集合,也就是正则表达式。

目前最快速的文本搜索工具grep就内置了正则表达式。grep起源于Unix中的ed编辑器的一条命令g/Regular Expression/p, 读作“Global Reular Expression Print”,即运用正则表达式的全局输出。由于这个功能太过实用,于是就从ed中独立出来,演变成了grep以及扩展版本的egrep。都知道grep因为有正则表达式所以很强大,但是正则表达式为何如此强大呢?

正则表达式的强大之处在于它是一套语法,分为两个部分,元字符(metacharacters)普通文本字符(normal text characters)

以语言类比,“我爱正则表达式”这句话可以抽象成“主谓宾”结构,主语是”我”,谓语是”爱”,宾语是“正则表达式”。这种语法还适用于其他语言,比如说英语就是”I love regular expression”. 这种语法结构就是元字符,而构成句子的语言就是普通文字字符。

元字符一览

正则表达式有很多流派,不同流派之间的差异在于对元字符的支持程度。以下的元字符适用于GNU版本的grep, sed, awk. mac自带的是BSD版本。

匹配单个字符的元字符

元字符 匹配对象
. 匹配单个任意字符
[…] 匹配单个列出的字符
[^…] 匹配单个未列出的字符
\char 转义元字符成普通字符

提供技术功能的元字符

元字符 匹配对象
匹配0或1次
* 匹配0到n次
+ 至少一次,最多不限
{min,max} 至少min次, 最多max次

匹配位置的元字符

元字符 匹配对象
^ 匹配一行的开头
$ 匹配一行的结尾

其他元字符

元字符 匹配对象
竖线 匹配任意分割的表达式
(…) 限定多选结构的范围,标注量词的作用范围,为反向引用捕获元素
\1, \2 反向引用元素

在grep中使用正则表达式

grep的强大之处它所做的事情就只有在文本搜索”正则表达式“定义的**模式(pattern)**,如果找到就打印出来。可以使用man egrep查看所支持的参数。

1
2
3
4
5
6
7
8
9
10
11
12
egrep [options] pattern [file]
egrep [options] [-e pattern]... [-f FILE]... [FILE...]
# 参数参数
-e PATTERN: 定义多个模式
-f FILE: 从文本中读取模式
-w: 匹配整个单词
-v: 反向匹配
-i: 忽略大小写
-x: 仅仅选择整行匹配结果
-c: 计数
-n: 输出表明行号
-A/-B NUM: 同时输出后/前几行

: grep有基础和扩展两个模式,基础模式支持的元字符较少,而egrep表示扩展的grep,支持的元字符较多。

检查测序结果的接头

在分析的质控阶段,需要检查给出的结果是否含有接头(adapter)。除了fastqc, 还能用grep检测是否有接头序列。

你可以去illumina的官网根据公司测序的protocol查找对应的接头:

https://support.illumina.com/downloads/illumina-customer-sequence-letter.html

或者是直接看FASTQC的配置文件中检测的接头

adapter

1
2
# 下载测试数据
fastq-dump -X 10000 SRR1553606 --split-files

然后用部分的接头对fastq文件进行搜索

1
egrep -B1 'TCGGAA' SRR1553606_1.fastq

pattern search

被找到的序列并非出现最开头,你可能希望是开头有几个其他碱基,然后跟着接头序列.

1
egrep -B1 '^[ATGC]{0,5}TCGGAA' SRR1553606_1.fastq

enhance pattern

基因搜索

一般而言,大家都回去TAIR上查找基因的区间。但是实际上我们可以先下载好拟南芥参考基因组及其注释文件,然后直接用正则表达式确定区间,用bedtools提取序列进行了。

1
2
3
4
5
6
# 下载序列
wget -c -4 -q http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas
# 下载GFF文件
wget -c -4 -q http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff
# 基因的信息
wget http://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/gene_description_20131231.txt.gz

在注释文件中查找某一个基因如SPL9, 并获取他的位置信息.

1
grep -i -w 'spl9' gene_description_20131231.txt

image.png

1
grep  -e 'AT2G42200.1' TAIR10_GFF3_genes.gff | grep 'mRNA'

image.png

当然以上命令可以连成管道, 并且和bedtools合用,就可用直接获取DNA序列。

1
bedtools getfasta -fi TAIR10.fa -bed  <(grep -i -w 'spl9' gene_description_20131231.txt | cut -f 1 | xargs -i grep -i {} TAIR10_GFF3_genes.gff | grep mRNA | cut -f '1,4,5')

你可以将其定义成一个函数,存放在配置文件中, 随后就能进行调用该函数。

bash function

use function