在我没有启动多少应用的时候，macOS已经显示它使用了22.09GB内存。其中App内存是15.81GB, 我并没有打开那么多App.

这估计跟configd有关，因为configd占用了20.55G内存。那么configd是真的占用了内存，还是就是声明自己会用到那么多内存呢？

我尝试着调用了比较多的内存，直接用了29.3Gb内存

1
2
3
4
5
6
7
8
9

cols <- 8189
rows <- 320127
mat1 <- matrix(data = 0, nrow=320127, ncol = 8189)
print(object.size(mat1), unit="GB")
# 19.5 Gb
mat2 <- matrix(data = 0L, nrow=320127, ncol = 8189)
print(object.size(mat2), unit="GB")
# 9.8 Gb

此时看内存，发现rsession-arm64内存到了29.55GB, 而configd还是20.55GB，已使用内存仅仅比之前多了不到4GB。

在我操作的前后，内存压力始终是绿色，说明系统没啥压力，那看来，有压力的显然是我。

我显然是之前用Windows系统，被某些软件的绿色气泡给吓到了，当气泡变红时，我就压力变大，就得点这个球释放内存。

而Unix系统的思路是调度所有内存，所有内存都要尽可能去干活，不要空着。

通过查阅文档，我们也可以知道configd是macOS特有的系统配置的后台驻留程序(System Configuration Daemon).configd管理本地系统的各类配置，维持数据映射到目标位置和当前系统状态，当数据发生改变时给应用程序发送信息，同时以 loadable bundles 形式托管一组配置代理。

最后，既然内存没压力，就不要自己给自己找压力，想着释放内存了。

目前(截止到2022年2月份), Arm64架构的R能够直接安装CRAN绝大部分的R包，但是Bioconductor的R包都需要编译（这也就为什么我推荐用Intel架构的R，可以避免遇到要编译R包里的坑）。

编译R包最怕遇到的问题，就是遇到R以外的底层库的依赖，比如说Rhtslib, 就要求编译htslib, 而htslib依赖lzma是MacOS默认没有安装，因此我们会因为缺少xz库的支持，导致htslib编译失败，从而导致Rhtslib安装失败，间接导致所有依赖于Rhtslib的R包无法使用。还有一些R包编译需要gsl库，就会因为找不到 gsl-config 而报错。

经过我不断踩坑后，我整理了解决方案。后续，如果你也打算在Arm64的Mac上安装原生的R，就按照我的思路来配置

首先，你需要安装homebrew, 国内用户推荐按照清华镜像里的方案来配置，https://mirrors.tuna.tsinghua.edu.cn/help/homebrew/ , 能够避免网络问题导致的坑。在海外的朋友，就按照官方站点, brew.sh 的命令来配置。

homebrew会安装在 /opt/homebrew下，该目录里面包括后续安装的软件(bin), 头文件(include)和动态链接库(lib).

之后，我们需要为R配置相关环节，使得R能够调用到homebrew安装的环境。

我们需要打开R，然后运行如下命令

1

file.edit(file.path(Sys.getenv("R_HOME"), "etc", "Makeconf"))

如此就能打开R的Makeconf文件，R就通过该文件里面的参数来进行R包和底层依赖的编译。我们搜索搜索 CFLAG, 添加 -I/opt/homebrew/inlucde， 搜索LDFLAGS, 添加 -L/opt/homebrew/lib, 效果如下

1
2

CFLAGS = -I/opt/homebrew/include -falign-functions=64 -Wall -g -O2 $(LTO)
LDFLAGS = -L/opt/R/arm64/lib -L/opt/homebrew/lib

重启R/RStudio之后，你就能编译安装Rhtslib（当然你需要先运行 brew install xz，把lzma.h 装好`)

1

BiocManager::install("Rhtslib")

此外你可能还会遇到gsl-config not found报错，尽管你明明用brew install gsl安装了软件，也能够在命令行中成功调用 gsl-config，但是R就是视而不见。 这是因为R读取的变量，默认就来自于MacOS的环境变量来自于 /etc/path文件和 /ect/path.d目录里。 因此，你即便在.zshrc, 或者.bashrc里面加入了 homebrew 的路径，R也是不认的。

解决方案，打开R，然后运行如下命令,

1

file.edit("~/.Rprofile")

添加如下内容，修改R启动时的环境变量

1
2

Sys.setenv("PATH"=paste0( "/opt/homebrew/bin:/opt/homebrew/sbin:",  Sys.getenv("PATH")))

重启R/Rstudio之后，R就能够用到homebrew安装的程序了。

或者，我们也可以下载到本地，用 R CMD INSTALL进行编译安装，以DirichletMultinomial为例

1
2

wget https://mirrors.tuna.tsinghua.edu.cn/bioconductor/packages/3.14/bioc/src/contrib/DirichletMultinomial_1.36.0.tar.gz
R CMD INSTALL DirichletMultinomial_1.36.0.tar.gz

总结一下:

1. 通过homebrew来安装依赖环境
2. 编译R的Makeconf使得R能够调用homebrew里的include和lib
3. 编译R的.Rprofile修改PATH，使得R能够用到hombrew的bin的程序

当使用biopython的SeqIO模块加载fasta之后，每条fasta都会记录成 Bio.SeqRecord.SeqRecord, 会有id, name, description, dbxrefs 等属性.

举个例子，输入数据是 test.fa

1
2
3
4
5
6

# test.fa
#>a|b c
#ATCG
seq = SeqIO.read("test.fa", "fasta")
seq
SeqRecord(seq=Seq('ATCG'), id='a|b', name='a|b', description='a|b c', dbxrefs=[])

我们可以将其保存成另一个文件

1
2

with open("test2.fa", "w") as f:
    SeqIO.write(seq, f, "fasta")

输出的header和输入的header是完全一样的。

假如，你想修改他的名字, 比如说把 “a|b” 改成 “a”, 你会发现一个非常诡异的事情，就是它完全不会按照你想的来.

1
2
3
4
5
6

seq.id = "a"
# 输出 a a|b c
seq.name = "a"
# 输出a|b c
seq.description = "a c"
# 输出 a|b a c

如果你改了他的ID，他会输出 ID + description, 如果你改了name，相当于没改，如果你改了description，就真的只有description（即header中空格后的信息）会发生变化。

假如我想输出 “>a c”, 那我就必须下面这个形式

1
2

seq.id = "a"
seq.description = "c"

我实在是没想懂他背后的逻辑, 修改输出的header居然如此不符合预期

代码和数据都在GitHub上，见 https://github.com/xuzhougeng/rice_annotation

水稻有两个主要的注释项目组

• RAP-DB: 编号格式形如Os01g0100100
• RGAP: 编号格式形如LOC_Os0101010

虽然这两个注释项目组使用的是同一个水稻基因组，但是注释结果却存在一些差异。你可以使用任意一套注释进行分析，通过ID的对应关系，将其中一套编号转换成另外一套编号，我用的是RAP-DB版本，这是因为KEGG数据库使用的就是RAP-DB的编号形式。

RAP-DB的网站上提供了很多数据方便我们使用，比如说光基因序列文件就至少提供了为UTR+CDS, UTR+CDS+Intron, CDS三类，非常贴心。美中不足的一点是，他们提供的GFF文件并不规整。

他们提供了2个下载链接，分别对应GFF和GTF

• Gene structure and function information in GFF format.
• Gene structure (only exon) information in GTF format.

GFF下载之后是一个压缩包，解压能得到locus.gff, transcripts_exon.gff, transcripts.gff. GTF下载之后只有一个gtf文件。

然而，无论是GFF还是GTF文件，他们都没有提供完整信息。GTF里面只有transcript和exon. GFF要么只有gene(locus.gff), 要么只有mRNA和exon(transcripts_exon.gff), 要么是mRNA,UTR和CDS(transcripts.gff). 我就非常纳闷了，为什么没有一个文件，同时有gene, mRNA, UTR, CDS, exon这些完整信息呢？

算了，既然他们没有给，我就自己动手，于是我写了一些代码，基于他们提供loucs.gff和transcripts.gff进行整理，代码我上传到GitHub。

我在代码中，将chr01-12改成了Chr1-Chr12, 所以对应的基因组文件，我也手动调整了。整理好的FASTA, GTF, GFF文件，我都一并上传到GitHub上，方便后续使用，地址为 https://github.com/xuzhougeng/rice_annotation

既然整理好了FASTA和GFF文件，为了方便后续在R语言里面使用，我还生成了对应BSgenome和TxDb，安装方法如下（生成代码和输出R包都在GitHub上）

安装方式(以BSgenome为例):

1. 在R里安装：install.pakcages("/路径/到/BSgenome.Osativa.MSU.xzg_1.0.tar.gz", repos=NULL, type="source)
2. 在命令行里安装: R CMD INSTALL BSgenome.Osativa.MSU.xzg_1.0.tar.gz

为了避免和Biocondutor其他MSU的包冲突，所以我用自己名字的缩写xzg作为后缀。

安装之后，就可以通过library进行加载

1
2

library(BSgenome.Osativa.MSU.xzg)
library(TxDb.Osativa.MSU.xzg)

在我用mac的时候，我一直有一个问题，那就是mac键盘上一些按键为什么会用哪个符号表示？比如说，最常见的command键用的⌘表示。通过一番检索，我找到了一个视频The Secret History of Mac Keyboard Keys ，这里简单的整理成文字，进行介绍下。

command(⌘),最初用的是苹果的logo, 称为apple key. 但是乔布斯觉得随处使用苹果的logo不利于用户区分, 就让设计师重新找logo, 最终设计师看到这个最早用在Scandinavia的高速公路旁的露营地上使用的图标, 替换了原来的apple logo。

option(⌥), 在windows上对应alt键，两者都有切换的含义, 所以这个图标可以想象成铁道上两条轨道。

control(^)用的符号叫做脱字符(caret). windows和mac都有control键，但是windows的control通常对应mac上command键, 这就让我们在不同系统间切换变得头疼了. 通常而言在终端上, windows和mac的control键时相同含义，但是对于其他跨平台的软件, windows用的control键到mac上都要改成command键。之所以用脱字符, 是因为这个符号指向上方, 让你去脱离当下掌控(control)全局。

脱字符在我们小时候写作文时比较常用,我们写完一句话之后,发现漏掉一个或几个字,于是就用这个符号加我们需要补充的内容, 把脱落掉的文字补充回去。

Caps Lock(⇪) 和shift(⇧) 键, 在英文输入模式下, 前者用于锁定大写输入, 后者用于临时输入大写字母, 在中文输入模式下, 前者默认功能是切换中英文, 后者还是临时输入大写字母. 基于英文输入语境, 这两个功能的符号就可以理解成保持大写和临时大写了.

delete(⌫) 对应windows上backspace键, 功能和符号我觉得是比较对应的, 就是删除(叉掉)内容. 有一点要注意到是 windows系统可以用键盘上的Del删除文件, 但是mac删除文件的快捷键默认是command + del, 对应windows键盘的win + backspace, 而不能用win + Del. 这是因为Del对应的是⌦, 表示向后删除.

最近切换到了MacOS平台进行办公，就不能用Windows下好用的XShell，用上了传统在命令行输入 ssh -p port user@address的方式进行登录了。

作为一个‘懒惰’的人，我肯定是要避免重复的运行登录命令了。回溯用过的命令进行复用是一种方式，但还是需要输入密码，所以我的操作方式如下

第一步: 通过编辑 ~/.ssh/config文件， 为指定服务器增加别名

1
2
3
4

Host 别名
    HostName 服务器地址
    User 用户名
    Port 端口

这样子就能用 ssh 目标服务器的别名的方式登陆指定服务器，不必写后面的内容

第二步: 将本地的ssh公钥上传到目标服务器，实现免密登陆。

1
2
3
4

# 生成密钥
ssh-keygen
# 上传公钥
ssh-copy-id 目标服务器的别名

通过上述两步的配置，后续登陆服务器就就能节约不少时间，同时基于SSH的scp在服务器之间传送也不需要输入密码了。

视频地址: https://www.bilibili.com/video/BV1F3411E7td/

各位好，我是洲更，欢迎收看本期教程。

在这一期视频中，我将会根据我的个人经验介绍如何在M1芯片的Mac上配置R语言的分析环境。

自从苹果在2021年10月份发布搭载M1 Pro或者是M1 Max芯片的14寸和16寸的MacBook Pro之后，目前在苹果官方网站就只能够购买M1芯片的MacBook Pro了。也就意味着，往后我们购买到的Mac都是ARM加购，而不是原来intel的x86架构。

我在付完钱，等待2个星期之后，终于也在2021年的11月的最后一天，到我手上这台搭载M1 Max芯片的16寸的MacBook Pro 。

于是，我终于有设备可以让我介绍，如何在M1芯片的电脑上配置R语言的分析环境了。

插句题外话：不得不说，苹果的供应链系统真的强，只要你愿意等，你就能用它标的价格购买，而不像我之前想买的PS5和英伟达的显卡，要么就是几万人抢几百台，要么就得要多花个几百几千从其他店铺买。

在这次教程中，我会介绍如何在M1芯片的Mac电脑上搭建R的分析环境，主要包括如下几个部分

• R语言的安装
• RStudio的安装
• R包安装所需的环境，
• R包运行环境的配置
• R包安装

尽管，我主要是以M1芯片的Mac介绍如何配置R语言环境，但由于目前许多软件都没有跟上，所以下载软件大部分还是为intel芯片开发的，需要通过rosetta2转译后才能在M1芯片的mac使用，所以这篇教程也适用于搭载intel芯片的mac电脑。

R语言的安装

首先，我们需要安装最核心的组件，也就是R语言。

我们先打开浏览器，搜索R，找到R语言的官网,http://r-project.org, 选择Download R，国外用户选择主站，国内用户推荐选择清华镜像源

然后点击 Download R for macOS， 此处，我们就会有两个选择， 一个是intel版本的R，一个是arm64版本，也就是专门给M1芯片开发的R。

那么到底要选择哪个呢？

因为我这个视频是面向初学者的，比较基础，所以我无脑推荐使用Intel 64-bit的R语言，牺牲一点性能去避免不必要的折腾，是非常值得的。

后续，等到M1芯片相关的R语言环境变好了，我再出一期视频介绍。

下载完之后，我们只需要不断的下一步，就可以了。

RStudio安装

说完了R语言的安装，RStudio的安装就非常的容易了，因为他没有选择，目前就只有一个intel的版本。

同样通过搜索找到Rstudio的官方网站，然后找到他的下载方式，然后点击下载Mac版本的Rstudio。下载地址为https://www.rstudio.com/products/rstudio/download/#download

下载完成后，我们双击打开，将其拖动到Applications中。

然后打开启动台，等待Rstudio出现，然后点击它打开它。

这里会出现一个弹窗，问我们是不是需要安装一个命令行开发工具，我们选择取消。后续我们会专门去安装这个工具。

我们简单的打一个 print(“hello world”)测试即可。

在安装完R和Rstudio后，我还想给大家推荐一个非常好用的工具， Rswitch

MacOS版本的RStudio是不支持切换不同的版本的R的，安装了这个工具后，我们就能够实现类似于Windows那种不同版本R语言切换的功能，

然而，这个软件比较小众，百度上我们暂时搜索不到（不过微软的必应可以），不过，这里我们选择直接输入网址， (https://rud.is/rswitch/)。

我们点击下载，等待下载完成，最后将其拖动到应用程序中就算安装完成了。

我们在启动台里将其打开，就可以在屏幕的右上角找到它。点击它，就能够用它来切换不同版本的R，以及用它来启动Rstudio。

编译环境

在编译环境上，我们需要安装2个工具， xcode和gfortran。

xcode是苹果提供的开发工具，里面包括编译C和C++代码所需的工具。安装过程很简单，

• 首先，我们打开启动台，
• 然后，在其他里面，找到终端，并打开
• 使用快捷键cmd + + 就可以放大终端
• 在终端里输入 sudo xcode-select --install，敲回车键
• 此时提示我们需要输入密码，注意，为了保证安全，所以在输入密码的时候，是看不到我们输入内容，我们输入结束后，敲回车就可以确定了。
• 接着就有一个窗口弹出来，我们点击同意。

安装时间取决于网络环境，所以你的时间可能和我的不一样。我们等待它完成。

安装成功之后，我们在终端里面输入如下命令，如果没有提示找不到，就说明安装成功了。

1
2

which clang
which clang++

除了xcode外，R包编译还会用到gfortran。跟xcode不同，gfortran的安装需要分两种情况，一种是为intel版本的R准备的版本，一种是为M1芯片版本的R准备的版本。

那么gfortran从哪里下载呢？我们先在浏览器中输入这个地址，https://mac.r-project.org/tools/index.html，打开之后，就可以看到下载方式了。

因为我们安装的是intel 64位的R，所以我们需要选择intel的版本下载就行。

如果你发现自己无法打开，可以尝试https://download.csdn.net/download/u012110870/78708981

intel的版本安装就是不断的下一步。

运行环境

接下来，我们来准备R语言的运行环境。

在运行环境上，我们需要安装两个工具，一个是Java,一个是，XQuartz,

有些包在安装的时候看起来没啥问题，但是在运行的时候却会报错，比如说rJava，究其原因，就是因为我们没有安装java环境。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

library(rJava)

The operation couldn’t be completed. Unable to locate a Java Runtime.
Please visit http://www.java.com for information on installing Java.

错误: package or namespace load failed for ‘rJava’:
 loadNamespace()里算'rJava'时.onLoad失败了，详细内容：
  调用: fun(libname, pkgname)
  错误: JVM could not be found
此外: Warning messages:
1: In system("/usr/libexec/java_home", intern = TRUE) :
  运行命令'/usr/libexec/java_home'的状态是1
2: In fun(libname, pkgname) :
  Cannot find JVM library 'NA/lib/server/libjvm.dylib'
Install Java and/or check JAVA_HOME (if in doubt, do NOT set it, it will be detected)

那如何安装Java呢？这时候我们又需要借助于清华镜像了。

我们打开浏览器，搜索清华镜像，然后打开，选择 AdoptOpenJDK，我们选择最新的17，选择JDK，我们装的是intel版本的R，因此点击x64，点击mac，点击这个tar.gz结尾的文件进行下载。

下载结束之后，我们需要安装我们之前下载的安装包。这一步需要用到终端，请跟着视频输入命令。

这里，我选择将文件拖动到终端中，来获取文件路径

1
2
3

sudo mkdir -p /opt/jdk_x64
sudo tar xf OpenJDK17U-jdk_x64_mac_hotspot_17.0.1_12.tar -C /opt/jdk_x64

按照上述步骤装完Java只完成一半，另外一半，我们需要在RStudio里操作。

我们打开Rstudio，输入如下命令 file.edit("~/.Rprofile") 打开R的配置文件，然后属于如下内容。由于你下载的版本可能和我的不一样，所以在输入过程中，要不断的使用tab 补全功能，确保路径正确。

1

Sys.setenv("JAVA_HOME"="/opt/jdk_x64/jdk-17.0.1+12/Contents/Home")

重启Rstudio后，R语言就会根据我们的配置来确定Java所在目录下。

除了java外，XQuartz，是后续安装的R包rgl不可或缺的一部分。如果没有安装XQuartz，即便安装了rgl，加载也会出现如下报错。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

> library(rgl)
Error in dyn.load(dynlib <- getDynlib(dir)) : 
  无法载入共享目标对象‘/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so’：:
  dlopen(/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so, 0x0006):
   Library not loaded: /opt/X11/lib/libGLU.1.dylib
  Referenced from: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/rgl/libs/rgl.so
  Reason: tried: '/opt/X11/lib/libGLU.1.dylib' (no such file), '/usr/lib/libGLU.1.dylib' (no such file)
错误: package or namespace load failed for ‘export’:
 loadNamespace()里算'rgl'时.onLoad失败了，详细内容：
  调用: rgl.init(initValue, onlyNULL)
  错误: OpenGL is not available in this build
此外: Warning messages:
1:   Loading rgl's DLL failed. 
  This build of rgl depends on XQuartz, which failed to load.
 See the discussion in https://stackoverflow.com/a/66127391/2554330 
2: Trying without OpenGL...

我们打开浏览器，输入它的官方站点: https://www.xquartz.org/

XQuartz目前是2.8.1版本，我们选择下载，双击打开，开始安装。安装结束之后，mac会注销一下。之后，我们打开启动台，打开其他，就能看到它了。

R包安装

在完成环境配置之后，我们终于可以安装R包了。

我们打开Rstudio，为了提供安装速度，我建议国内的用户修改一下镜像，提高安装速度，推荐国内的用户，按照如下方式修改修改。

首先，我们打开R的配置文件，然后我们打开清华镜像，点击CRAN旁边的问号，然后复制这一行到 R的配置文件里；接着我们找到Bioconductor，复制这一行到R的配置文件里。注意用cmd + s进行保存。

1
2
3
4
5
6

file.edit("~/.Rprofile")

options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")

之后，我们关掉配置文件，并重启Rstudio。

我们新建一个R脚本。

我们首先安装一个tidyverse，来测试CRAN的包的安装。

输入install.packages(“tidyverse”), 然后用快捷键 cmd + 回车运行。

首先蹦跶出一堆红字，这是它依赖的其他R包。

然后提示 有一个R包，它虽然有二进制版本，但源代码是最新的，问我们是不是要编译最新的。因为我们已经配置好了编译环境，所以这里输入yes。

后续，又出现了大量的红字，虽然很红，但是不是报错，都是正常信息。此时出现了黑字，这些事编译信息。

等他安装结束之后，我们输入 library(tidyverse) 来确认是否安装成功。 虽然有很多信息，但是没有看到error就说明成功了

接着，我们安装xlsx 来测试java环境。 没有报错说明，java环境正常。

然后是 rgl；我这里提示R出现fatal error，这并不是我们的环境配置出问题了，而是因为我的操作是在虚拟机里面执行的，是虚拟机性能不行，在加载时出错了。

接着，我们安装， BiocManager 用来安装 Bioconductor上的包。安装成功后，我们安装DESeq2，来测试Bioconductor，同样它也有很多依赖。

安装结束之后，提醒我们是否需要升级包，我这里输入n 表示不想升级，你们可以试试输入y。

最后我们试试安装github上的包，我们以ComplexHeat为例。我们输入如下命令， 运行的时候提示devtools没有找到，这说明我们没有安装这个包。

我们先去安装这个包； 这里有提示没有 这个函数， 原因是我输入出错了。改正就行了。

安装devtools后，我们再次运行之前的命令，这次就没有问题了。 我们加载也是成功的。

参考资料

• https://mac.r-project.org/
• https://mac.r-project.org/tools/index.html
• https://kaopubear.top/blog/2021-07-06-macos-arm64-R-tips/
• https://code2care.org/q/install-native-java-jdk-jre-on-apple-silicon-m1-mac

GMAP最早用于讲EST/cDNA序列比对到参考基因组上，可以用于基因组结构注释。后来高通量测序时代，又开发了GSNAP支持高通量数据比对，这篇文章主要介绍GMAP，毕竟高通量转录组数据比对大家更喜欢用STAR, HISTA2等软件。

软件安装

下面是我源码安装的代码

1
2
3
4
5

wget http://research-pub.gene.com/gmap/src/gmap-gsnap-2018-07-04.tar.gz
tar xf gmap-gsnap-2018-07-04.tar.gz
cd gmap-2018-07-04/
./configure --prefix=$HOME/opt/biosoft/gmap
make -j 20

软件使用

如下步骤假设你有一个物种的基因组序列和对应的CDS序列，分别命名为”reference.fa”和”cds.fa”

第一步:构建GMAP/GSNAP索引数据库

GMAP/GSNAP对FASTA文件中每个记录下的序列的长度有一定限制, 每一条不能超过4G, 能应付的了大部分物种了。

构建索引分为两种情况考虑，第一种是一个fasta文件包含所有的序列

1

~/opt/biosoft/gmap/bin/gmap_build -d reference reference.fa

第二种则是每个染色体的序列都单独存放在一个文件夹里，比如说你下载人类参考基因组序列解压后发现有N多个fasta文件, 然后你就想用其中几条染色体构建索引

1

~/opt/biosoft/gmap/bin/gmap_build -d reference Chr1.fa Chr2.fa Chr3.fa ...

注: 这里的-d表示数据K库的名字，默认把索引存放在gmap安装路径下的share里，可以用-D更改.此外还有一个参数-k用于设置K-mer的长度, 默认是15, 理论上只有大于4GB基因组才会有两条一摸一样的15bp序列(当然是完全随机情况下)。

第二步:正式比对

建立完索引之后就可以将已有的CDS或者EST序列和参考基因组序列进行比较。

1

~/opt/biosoft/gmap/bin/gmap -t 10 -d reference -f gff3_gene cds.fa > cds_gene.gff3

其中-t设置线程数, -d表示参考基因组数据库的名字, 都是常规参数。我比较感兴趣的参数是如何将序列输出成GFF格式. GMAP允许多种格式的输出，比如说-S只看联配的总体情况，而-A会显示每个比对上序列的联配情况, 还可以输出蛋白序列(-P)或者是genomic序列(-E). 但是做结构注释要的gff文件，参数就是-f gff3_gene, -f gff3_match_cdna, -f gff3_match_est。

参考文献

要想对一个软件有更好的认识，最好还是看看他们文章是怎么说的。

• GMAP: a genomic mapping and alignment program for mRNA and EST sequences
  Bioinformatics 2005 21:1859-1875 Abstract Full Text, Thomas D. Wu and Colin K. Watanabe
• Fast and SNP-tolerant detection of complex variants and splicing in short reads
  Bioinformatics 2010 26:873-881 AbstractFull Text, Thomas D. Wu and Serban Nacu

本文基于biopython>=1.79

配对序列联配是一个非常基础的生信分析，biopython提供了Bio.Align.PairwiseAligner基于动态规划处理两个序列的联配问题。

如下是一个简单的案例：

1
2
3

from Bio import Align
aligner = Align.PairwiseAligner()
alignments = aligner.align("ACCGGTTT", "ACGGGTT")

它首先初始化了一个PairwiseAligner对象，然后用该对象对两条序列进行联配，联配结果记录在了alignments。

alignments里面记录得分最高的联配结果（不止一个），我们可以通过循环来遍历其联配结果。

1
2
3

for alignment in alignments:
    print("Score = %.1f:" % alignment.score)
    print(alignment)

PairwiseAligner的运行速度极快，这是因为他底层通过C语言实现，而非Python代码，如下是他的结构体。

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

typedef struct {
    PyObject_HEAD
    Mode mode;
    Algorithm algorithm;
    double match;
    double mismatch;
    double epsilon;
    double target_internal_open_gap_score;
    double target_internal_extend_gap_score;
    double target_left_open_gap_score;
    double target_left_extend_gap_score;
    double target_right_open_gap_score;
    double target_right_extend_gap_score;
    double query_internal_open_gap_score;
    double query_internal_extend_gap_score;
    double query_left_open_gap_score;
    double query_left_extend_gap_score;
    double query_right_open_gap_score;
    double query_right_extend_gap_score;
    PyObject* target_gap_function;
    PyObject* query_gap_function;
    Py_buffer substitution_matrix;
    PyObject* alphabet;
    int* mapping;
    int wildcard;
} Aligner;

在结构体中，有一个非常关键的algorithm, 决定了联配时所需的算法。主要分为三个大类 NeedlemanWunschSmithWaterman, WatermanSmithBeyer和Gotoh.算法无法手动指定，而是根据match, mismatch, gap等参数的得分来确定(共6种)，具体的决策过程如下

当gap得分与gap长度有关，也就是 gap_score 赋值一个函数时，在WatermanSmithBeyer基础上，当mode是global, 使用 “Waterman-Smith-Beyer global alignment algorithm”；当mode是local，使用 “Waterman-Smith-Beyer local alignment algorithm”

当gap打开的得分和gap延伸的得分相同时，在NeedlemanWunschSmithWaterman基础上，如果mode是global， 使用Needleman-Wunsch；如果mode是local， 使用Smith-Waterman

其他情况下，也就是得分不同，则是在Gotoh基础上，如果mode是global，选择 “Gotoh global alignment algorithm”，如果mode是local，则选择”Gotoh local alignment algorithm”

我们举几个简单的例子来说明

场景1: 通过函数定义为gap打分

1
2
3
4
5
6
7
8

from Bio import Align
aligner = Align.PairwiseAligner()
aligner.gap_score = lambda x: x*2
aligner.algorithm
# 'Waterman-Smith-Beyer global alignment algorithm'
aligner.mode = "local"
aligner.algorithm
# 'Waterman-Smith-Beyer local alignment algorithm'

场景2: gap打开和延伸得分相同

1
2
3
4
5
6

from Bio import Align
aligner = Align.PairwiseAligner()
aligner.open_gap_score = -1
aligner.extend_gap_score = -1
aligner.algorithm
# 'Needleman-Wunsch'

场景3: gap打开和延伸得分不同

1
2
3
4
5
6

from Bio import Align
aligner = Align.PairwiseAligner()
aligner.open_gap_score = -1
aligner.extend_gap_score = -2
aligner.algorithm
# 'Gotoh global alignment algorithm'

除了为match和mismatch设置固定值外，我们还可以考虑使用得分矩阵为不同碱基或者氨基酸设置不同的得分，例如BLOSUM62

1
2
3
4
5
6
7

from Bio import Align
aligner = Align.PairwiseAligner()
# add matrix 
from Bio.Align import substitution_matrices
aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
# align
alignments = aligner.align("ACCGGTTT", "ACGGGTT")

biopython的替换矩阵来自于ftp://ftp.ncbi.nih.gov/blast/matrices/, 保存在Bio/Align/substitution_matrices/data 目录下

在搞定联配之后，我们如何获取联配后的结果？实际上对于一个 PairwiseAlignment 对象，例如 alignment, 我们调用 print(alignment)实际等价于 print(format(alignment)), 因此我们可以通过如下代码，获取联配后的query和target序列

1
2
3

seqs = format(alignment).split("\n")
target = seqs[0]
query = seqs[2]

最后总结下本文内容：

1. biopython的配对序列联配算法是C编写，运行速度极快
2. 联配算法会根据 gap score和mode 自动确定
3. 基于format函数可以获取联配后的序列