在使用R语言多年以后，我终于开始去学习Rcpp，利用C++来提高运行速度。其实当你能熟练的使用一门语言后，再去学一门新的语言，并没有想象中的那么难，更何况Rcpp把很多脏活累活都给包办了，在里面调用C++还是挺方便。

C++是一门静态编译面向对象的编程语言，R是动态解释性面向对象语言，那么有一个不同就在于，你需要先声明一个变量，才能调用该变量。而在声明变量的时候，你就会遇到一个R语言中不怎么思考的问题，我这个变量要存放什么样的数据呢？

数据类型

R的基本数据类型有六种 “logical”, “integer”, “numeric” (等价于”double”), “complex”, “character” 和 “raw”，但常用的就四种， “integer”, “numeric”,”logical”和”character”。

对于数值数据而言，C++的整型和浮点型可以定义多种精度，而R语言则简化成两种”integer”(32 bit)和”numeric”(53 bits)，并且默认情况下数值都是浮点型。

在数据结构上，R语言提供了向量(vector)，矩阵(matrix)，数组(array)，数据框(data.frame)和列表(list)，唯独没有标量。其中向量可以认为是R语言的基本单位，是其他数据结构的基础。

原本能够让C++处理R对象，以及将C++处理完结果转成R能识别的对象，你需要了解R的内部结构，但是Rcpp通过定义了一系列数据类型（如下），使得我们能够非常容易地在R和C++之间交换对象。

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
27
28
29

// Rcpp/include/Rcpp/vector/instantiation.h
namespace Rcpp{

    typedef Vector<CPLXSXP> ComplexVector ;
    typedef Vector<INTSXP> IntegerVector ;
    typedef Vector<LGLSXP> LogicalVector ;
    typedef Vector<REALSXP> NumericVector ;
    typedef Vector<REALSXP> DoubleVector ;
    typedef Vector<RAWSXP> RawVector ;

    typedef Vector<STRSXP> CharacterVector ;
    typedef Vector<STRSXP> StringVector ;
    typedef Vector<VECSXP> GenericVector ;
    typedef Vector<VECSXP> List ;
    typedef Vector<EXPRSXP> ExpressionVector ;

    typedef Matrix<CPLXSXP> ComplexMatrix ;
    typedef Matrix<INTSXP> IntegerMatrix ;
    typedef Matrix<LGLSXP> LogicalMatrix ;
    typedef Matrix<REALSXP> NumericMatrix ;
    typedef Matrix<RAWSXP> RawMatrix ;

    typedef Matrix<STRSXP> CharacterMatrix ;
    typedef Matrix<STRSXP> StringMatrix ;
    typedef Matrix<VECSXP> GenericMatrix ;
    typedef Matrix<VECSXP> ListMatrix ;
    typedef Matrix<EXPRSXP> ExpressionMatrix ;

}

以R语言的向量为例，Rcpp考虑到R语言所有可能的数据类型，提供了9种Vector。比如说”NumericVector”, “IntegerVector”, “LogicalVector”,” CharacterVector”就是对应着”integer”, “numeric”,”logical”和”character”这四种数据类型。因此，你可以根据你输入向量的数据类型来进行选择。

我们以一个简单的求和函数作为例子

1
2
3
4
5
6
7
8
9

double sumC(NumericVector v){

    double total = 0;
    int num = v.size();
    for (int i = 0; i< num; i++){
        total += v[i];
    }
    return total;
}

sumC这个函数读取一个数值向量，对其进行结合，返回一个浮点型标量。当然这个浮点型标量在R语言中就表现为一个长度为1的向量。

由于R语言数据类型和C++的不是一一对应，因此在处理过程中，有些时候需要对数据进行转换

1

std::string(x[i])

变量创建

我们可以通过下面这些方式创建向量

1
2
3
4
5
6
7
8
9
10
11

// v <- rep(0, 3)
NumericVector v (3);
// v <- rep(1, 3)
NumericVector v (3,1);
// v <- c(1,2,3) 
// C++11 Initializer list
NumericVector v = {1,2,3}; 
// v <- c(1,2,3)
NumericVector v = NumericVector::create(1,2,3);
// v <- c(x=1, y=2, z=3)
NumericVector v = NumericVector::create(Named("x",1), Named("y")=2 , _["z"]=3);

创建矩阵的方法如下

1
2
3
4
5
6

// m <- matrix(0, nrow=2, ncol=2)
NumericMatrix m1( 2 );
// m <- matrix(0, nrow=2, ncol=3)
NumericMatrix m2( 2 , 3 );
// m <- matrix(v, nrow=2, ncol=3)
NumericMatrix m3( 2 , 3 , v.begin() );

数据框和列表的创建依赖于已有向量，和R语言中创建形式相似。

1
2
3
4
5
6
7
8
9

// Creating DataFrame df from Vector v1, v2
DataFrame df = DataFrame::create(v1, v2);
// When giving names to columns
DataFrame df = DataFrame::create( Named("V1") = v1 , _["V2"] = v2 );

// Create list L from vector v1, v2
List L = List::create(v1, v2);
// When giving names to elements
List L = List::create(Named("name1") = v1 , _["name2"] = v2);a

元素获取和赋值

在选取元素之前一定要注意，C++是以0为基，而R是以1为基。

我们可以用[]和()进行数据选取向量中的元素并复制，支持利用数值向量或者逻辑向量来选取多个数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

// 新建向量
NumericVector v  {10,20,30,40,50};
// 设置向两名
v.names() = CharacterVector({"A","B","C","D","E"});
// 准备用于获取数据的向量
NumericVector   numeric = {1,3};
IntegerVector   integer = {1,3};
CharacterVector character = {"B","D"};
LogicalVector   logical = {false, true, false, true, false};
// 获取数据
double x1 = v[0];
double x2 = v["A"];
NumericVector res1 = v[numeric];
NumericVector res2 = v[integer];
NumericVector res3 = v[character];
NumericVector res4 = v[logical];
// 赋值
v[0]   = 100;
v["A"] = 100;
NumericVector v2 {100,200};
v[numeric]   = v2;
v[integer]   = v2;
v[character] = v2;
v[logical]   = v2;

对于矩阵而言，建议只用()进行数据选取，因为没有[row_index, col_index]的操作。

1
2
3
4
5
6
7
8
9
10

// 创建一个5x5的矩阵
NumericMatrix m( 5, 5 );
// 获取0,2的元素
double x = m( 0 , 2 );
// 选择第1行
NumericVector v = m( 0 , _ );
// 选择第3列
NumericVector v = m( _ , 2 );
// 选择第1-2行，3-4列
NumericMatrix m2 = m( Range(0,1) , Range(2,3) );

对于DataFrame和List，两则都只支持[]选择其中向量。

成员函数

成员函数(member function)是C++面向对象编程中的一个概念，通过调用成员函数可以对类进行操作。

举个例子，对于R语言而言，向量是可以有名字的，例如

1

x <- c(a = 1, b = 2, c = 3)	

R语言是一个面向对象的编程语言，所以这种命名向量其实认为是一个向量拥有了一个名为names的属性而已。

在R语言中attr函数能够增加对象的属性，因此上面代码等价于

1
2

y <- c(1,2,3)
attr(y, 'names') <- c('a','b','c')

而C++中的向量和矩阵拥有一个成员函数.name()就可以获取/修改/增加变量命名。

1
2
3
4

NumericVector addname2(NumericVector x, CharacterVector y){
    x.names() = y;
    return x;
}

更通用的是.attr函数,它能够修改和增加任意属性

不同数据类型拥有不同的成员函数，参考

• 向量: https://teuder.github.io/rcpp4everyone_en/080_vector.html#member-functions
• 矩阵: https://teuder.github.io/rcpp4everyone_en/100_matrix.html#member-functions-1
• 数据框: https://teuder.github.io/rcpp4everyone_en/140_dataframe.html#member-functions-2
• 列表: https://teuder.github.io/rcpp4everyone_en/150_list.html#member-functions-3
• 字符串: https://teuder.github.io/rcpp4everyone_en/170_string.html#member-functions-4
• 日期: https://teuder.github.io/rcpp4everyone_en/180_date.html#member-functions-5
• 时间: https://teuder.github.io/rcpp4everyone_en/190_datetime.html#member-functions-6

缺失值

同R一样，缺失值之间以及缺失值和正常值之间的比较没有意义，因此在写代码过程时要提前考虑到缺失值的可能，或者是用R语言处理完缺失值，把不含缺失值的数据作为C++函数的输入。

R语言采用比特模式对每一种数据类型进行标注，也就是针对每一种数据类型（Integer，Numeric，Character，Logical）保留一个比特作为缺失数据的标签值。

为了在C++中处理R的缺失值，Rcpp提供了四个变量对应R的四种数据类型

• NA_INTEGER: 整型，Integer
• NA_STRING: 字符串，Character
• NA_LOGICAL：逻辑性，Logical
• NA_REAL：双精度浮点型， Numeric

写一个函数进行讲解下缺失值的使用。

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

double sumC(NumericVector v, bool na_rm){

    double total = 0;
    for (int i = 0; i< v.size(); i++){
        if( ! NumericVector::is_na(v[i])){
            total += double(v[i]);
        } else if (NumericVector::is_na(v[i]) & na_rm){
            continue ;
        } else {
            return NA_REAL;
        }
    }
    return total;
}

NumericVector::is_na()用于判断是否为缺失值，不同的数据类型有不同的is_na。如果na_rm=FALSE，那么返回缺失值，也就是NA_REAL.

参考资料

• CharacterVector元素转string
• https://teuder.github.io/rcpp4everyone_en
• https://adv-r.hadley.nz/rcpp.html#rcpp-intro

服务器上安装RStudio-server

如果想在服务器上安装一个RStudio-server，你需要先保证自己拥有管理员权限，之后参考如何在服务器上安装最新的R安装R语言，一定要注意在./configure的时候加上--enable-R-shlib，否则后续会出错。

RStudio-server分为两种版本，一种是开源免费版，另一个是商业专业版本。个人觉得两者最大的区别在于，商业版支持在多个版本的R语言之间进行切换，而开源免费版不行。

CentOS篇

如果服务器安装的是CentOS/RedHat，那么需要保证它们的发行版本不等于6

从官方上下载rpm文件

1

wget https://download2.rstudio.org/server/centos6/x86_64/rstudio-server-rhel-1.2.5001-x86_64.rpm

如果是第一次安装，那么就是运行如下命令

1

sudo yum install rstudio-server-rhel-1.2.5001-x86_64.rpm

假如是需要升级RStudio-server，比如我原先的是1.1.456最新的是1.2.5001, 需要先暂停当前的服务并卸载

1
2

/usr/sbin/rstudio-server stop
yum remove rstudio-server

之后才是安装

1

sudo yum install rstudio-server-rhel-1.2.5001-x86_64.rpm

安装完成之后，我们可以通过修改/etc/rstudio/rserver.conf更改端口和R所在路径

1
2
3

www-port=8080 # 端口, 默认8787
www-address=0.0.0.0
rsession-which-r=/opt/sysoft/R-3.6.1/bin/R # 安装R的路径

修改完成之后，用rstudio-server restart重启服务，没有任何信息就表示安装成功了。

当然你要是不放心，你还可以用rstudio-server verify-installation来验证下，如果没有任何输出信息就表示安装成功，假如出现下面这条信息，意味着你需要先用rstudio-server stop先暂停服务。

1

Server is running and must be stopped before running verify-installation

其实最直接的方法就是直接访问”IP:端口”，能够出现RStudio的登陆界面就意味着安装成功了。

Ubuntu篇

我没有一台Ubuntu系统的服务器，只有一台Windows 10电脑有一个Linux子系统安装的是 Ubuntu 16.04.6 LTS。

下载Deb文件

1
2

sudo apt-get install gdebi-core
wget https://download2.rstudio.org/server/trusty/amd64/rstudio-server-1.2.5001-amd64.deb

如果不是第一次安装，需要是升级已有的RStudio-server，那么也需要先停用并卸载已有的RStudio-server

1
2

sudo /usr/sbin/rstudio-server stop
sudo apt-get remove rstudio-server 

然后安装

1

sudo gdebi rstudio-server-1.2.5001-amd64.deb

如果是在Windows的子系统下安装，会出现如下的警告，允许访问即可。

如果能够打开http://127.0.0.1:8787, 就说明安装成功了。

如果想修改RStudio-server的端口和调用R版本，参考CentOS篇

参考资料

• https://rstudio.com/products/rstudio/download-server/redhat-centos/
• https://rstudio.com/products/rstudio/download-server/debian-ubuntu/
• https://support.rstudio.com/hc/en-us/articles/216079967-Upgrading-RStudio-Server

R语言有一个自带的函数table能够统计输入变量中不同元素出现的次数，举个例子

1
2

d <- rep(c("A","B","C"), 10)
table(d)

果子老师曾写一篇推送，自己写了一个简化版的table，比R自带的table 运行的速度更快，如下

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

tableGZ <- function(x){
  if(sum(is.na(x)) == 0){
    data <- x
    input <- unique(x, fromLast = TRUE)
    dd <- sapply(input, 
                 function(x) {sum(data==x)})
    names(dd) <- unique(data, fromLast = TRUE)
    dd
  } else{
    data <- x[!is.na(x)]
    input <- unique(x, fromLast = TRUE)
    dd <- sapply(input, function(x){
      sum(data == x)
    })
    dd <- c(dd, sum(is.na(x)))
    names(dd) <- c(input, 'NA')
    dd
  }
}

我们通过运行1000次代码，来比较下两者的运行速度

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

bench::system_time(for ( i in 1:1000){
  tableGZ(d)
})
# 结果
process    real 
 42.9ms  42.4ms 

bench::system_time(for ( i in 1:1000){
  table(d)
})
# 结果
process    real 
  107ms   106ms  

在我的电脑上，果子老师的代码运行速度比R自带的table快了将近3倍。当然这是有原因的，因为R的table的代码功能更加复杂，能够比较多个变量之间的关系，例如table(d,d)。

既然是简单的统计每个元素的次数，那么我就想着能不能写出一个比果子老师速度更快的函数。 于是，我抽空写了下面的代码

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40

tableZG <- function(x){
  
  NA_pos <- is.na(x)
  NA_num <- sum(NA_pos)
  
  x <- x[!NA_pos]
  
  out <- vector(length = length(x))
  out_name <- rep(NA,  length(x))
  
  for (i in 1:length(x)) {
    
    for (j in 1:length(out_name)){
      if ( is.na(out_name[j]) ){
          out_name[j] <- x[i]
          out[j] <- 1
          break
        }
      
      if ( out_name[j] == x[i] ){
        out[j] <- out[j] + 1
        break
      } 
    }
  }
  
  if (NA_num > 0){
    na_end <- sum(!is.na(out_name))
    out_name[na_end + 1] <- 'NA'
    out[na_end + 1] <- NA_num
    
  } 
  na_pos <- is.na(out_name)
  out_name <- out_name[!na_pos]
  out <- out[!na_pos]
  names(out) <- out_name

  return(out)
  
}

虽然我的代码更长了，但是并没有让速度提高，反而比果子老师的代码慢，甚至还不如R自带的table。

1
2
3
4
5
6

system.time(for ( i in 1:1000){
  tableZG(d)
})
# 结果
process    real 
  148ms   148ms 

当然那么一长串代码并不是白写的，因为我特意避免了使用R特有的内容，所以代码能够很容易改写成C++代码使用Rcpp调用，从而提高速度

新建一个tableC.cpp文件，代码内容如下

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
27
28
29
30
31
32
33
34
35
36

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector tableC(CharacterVector cv){
  
  // initialize variable
  CharacterVector na = CharacterVector::create("NA");
  NumericVector out = rep(NumericVector::create(0), cv.size());
  CharacterVector out_name = rep(na, cv.size());
  int unique_num = 0;
  
  // 
  for (int i = 0; i < cv.size();i ++) {
    
    for (int j = 0; j < cv.size(); j++){
      
      if ( out_name[j] == "NA" ){
        out_name[j] = cv[i] ;
        out[j] = 1;
        unique_num += 1;
        break;
      }
      
      if ( out_name[j] == cv[i] ){
        out[j] = out[j] + 1;
        break;
      } 
    }
  }
 
  out.attr("names") = out_name;
  
  return out;
  
} 

然后在R里面用Rcpp这个C++代码，替换掉之前代码中的循环部分

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

Rcpp::sourceCpp("tableC.cpp")
tableZG <- function(x){
  
  NA_pos <- is.na(x)
  NA_num <- sum(NA_pos)
  
  x <- as.character(x[!NA_pos])
  res <- tableC(x)
  res <- res[!names(res) == "NA"]
  
  if (NA_num > 0){
    res <- c(res, "NA"=NA_num)
  }
  return(res)
}

于是这一次在C++的加持下，我写的table函数速度超过了果子老师的代码。

1
2
3
4
5
6

bench::system_time(for ( i in 1:1000){
  tableZG(d)
})
#结果
process    real 
 30.1ms  29.3ms 

最后总结一下：如果一个操作只需要做一次，那么速度可能并不是最重要的。因为即便是一个原本要花24小时的代码，提速10倍，只要2小时，你可能也会愿意等一等。但是如果这个操作需要重复很多次，上百次，上千次，甚至上万次，那么你就可能等不下去了。你就需要对代码中的一些限速步骤进行优化，比如说table这种多功能函数，你就可以自己用R写一个简化版的函数，替换掉原先的代码。

如果对速度有更高的要求，那么就需要用到C++进行代码重写了。学习C++其实并不会特别难，因为有一个Rcpp简化了许多操作，你只需要掌握几个最基本的语言特性，比如说C++需要先定义变量才能使用变量。

推荐阅读

• 给R使用者的C++最少必要知识
• Rcpp学习笔记之Hello World!

在使用R语言多年以后，我终于开始去学习Rcpp，利用C++来提高运行速度。其实当你能熟练的使用一门语言后，再去学一门新的语言，并没有想象中的那么难，更何况Rcpp把很多脏活累活都给包办了，在里面调用C++还是挺方便。

学习Rcpp的最重要一步是，运行一个”hello world!” 。如果能够运行”hello world!”就表明搞定了环境配置，后面就可以愉快的写代码了。

安装Rcpp的方式为，install.packages("Rcpp")， 安装过程中可能会出现一些问题，对于不同的操作系统需要做不同的准备工作，

• Windows: 安装Rtools
• Mac: 安装Xcode，需要在App store下载
• Linux: 需要有GCC的编译环境

之后就让我们写人生中第一个C++函数， hello,

1
2
3
4

int hello(){
    std::cout << "Hello, World!";
    return 0;
}

那么问题来了，这个函数应该如何才能让R语言调用呢？最简单的方式就是Rcpp的cppFunction

1
2
3
4
5
6
7
8

library("Rcpp")
cppFunction('
int hello(){
    std::cout << "Hello, World!";
    return 0;
}
')
hello()

上面语句中，cppFunction中的''中的内容是C++语句。在这条语句中，我们以int hello(){}定义了一个hello函数，这个函数不接受参数，返回一个整型。在函数里面一共有两条语句，每条语句都以;结尾。

第一条是调用了C++的标准库的cout, std::cout, 和R中以包名::函数名调用函数的方法类似。第二条则是return 0，返回结果。

除了利用cppFunction外，另外一种更常用的方法就是将代码放在其他文件中，然后用sourceCpp的方式读取。我们新建一个hello.cpp的文件，里面的内容如下(如果用Rstudio新建C++ 文件，它会提供一个模版用于修改)

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

#include <Rcpp.h>
using namespace Rcpp;
using namespace std;

//[[Rcpp::export]]
int hello(){
    cout << "Hello, World!";
    return 0;
}

/*** R
hello()
*/

之后在R里面加载并调用，和之前的结果一样。

1
2
3
4

library(Rcpp)

sourceCpp("hello.cpp")
hello()

解释下C++的代码。第一行是用#include <Rcpp.h>语句导入了Rcpp的头文件，类似于R里面的install.package()用于安装R包，然后用using namespace xxx; 的方式加载了两个库，Rcpp和std, 这就类似于R里面的library()函数。

C++里面用//和/*** 注释语句 */进行代码注释。前者是注释单行，类似于R里面的#注释，后者是注释多行。只不过//[Rcpp::export]]这条注释有特殊的含义，在sourceCpp读取代码解析的过程中，被这条语句注释的函数能够在R里面调用。换句话说，如果你删了这句话，那么这个hello函数在R里面就是无法直接调用的。而/*** R */里面可以放R代码，会c++代码编译结束后运行，常用于代码测试。

假如我们能够成功运行上面的代码，那么接下来要做的事情就是学习C++的基本语法(参考给R使用者的C++最少必要知识)，学习Rcpp的数据结构。

C++是一门非常复杂的编程语言，但是如果你已经有一定的R语言基础，希望通过C++来提高代码效率，那么掌握下面几点就能够开始写代码了，然后通过Rcpp调用。

C++是静态编程语言

C++是一门静态编程语言，这意味着对于一个变量而言，你需要先声明它，才能调用它。而且这个变量名的数据类型在使用过程中是无法更改的，这是因为它在内存中的大小已经固定了。

举个例子，在R语言中，下面这个代码运行时不会报错（但是会被吐槽）

1
2

a <- 1
a <- 'a'

但是在C++中，下面的代码会被提示error: redefinition of 'x' with a different type: 'char' vs 'int'

1
2

char a = 'a';
int a = 0;

注：每个C++语句后需要有;, 而R没有。

关于数据结构，标准C++的数据结构有如下几种

• 整型: int, long
• 浮点型: float, double
• 逻辑: bool
• 字符: char

C++的字符串在STL(Standard Template Library)中。STL里有很多高级数据结构，例如向量(vector)和列表(List)限于篇幅请自行检索。

C++的控制结构

C++的控制结构包括，for循环，while循环，if条件语句，switch条件语句，代码形式如下

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
27
28
29
30
31
32
33
34

// for loop
for ( int i = 1; i <= 10; i++){
    ...;  
}

// while loop
while ( i <= 10){
  ...;
}

// if
if (){
  ...;
} else if(){
  ...;
} else{
  ...;
}

// switch
char grade = 'D';
switch (grade ){
  case 'A':
    cout << 'great' << endl;
    break ;
  case 'B' :
    cout << 'not bad' << endl;
    break ;
  case 'D' :
    cout << 'Not good' << endl;
  default:
    cout << 'wrong' << endl;
    
}

在代码形式上，if和while都是和R语言类似，for和switch则是有些不同。

C++的循环效率远远高于R语言，而且将对应的R代码修改成C++并不复杂。因此会用C++的循环，就能优化很多R代码。

函数

C++的函数和R不同，它的写法如下

1
2
3
4

数据结构 函数名(数据结构 参数, ...){
  代码;
  return 变量名;
}

举个例子

1
2
3
4

int sumTwo(int x, int y){
  int z = x + y;
  return z;
}

如果是在R中则是

1
2
3
4

sumTwo <- function(x, y){
   z <- x + y
   return(z)
}

因此在改成原有的R函数时要注意两者的区别，

• 代码形式不同
• 参数名要声明数据类型
• return后直接跟变量名

在函数的使用上，C++和R就没有区别了，都是函数名(变量)，举个例子sumTwo(1,2)

额外库加载

在使用R语言的大部分时间里，我们都是面向R包编程，也就是搜索一个R包，然后安装R包，加载R包，调用函数。对于C++而言，也有许多现成的库，能够让我们避免造轮子。

举个例子，加载Rcpp库，调用R函数

1
2
3
4
5
6
7
8
9
10
11

#include <Rcpp.h>
#include <cstdio>

using namespace Rcpp;

//[[Rcpp::export]]
int main(){
  printf("N(0,1) 95th percential %9.8f\n",
         R::qnorm(0.95, 0.0, 1.0, 1, 0));

}

这里#include类似于安装R包，而using namespace Rcpp;则是library加载R包。

面向对象

C++和R语言都是一门面向对象编程语言。R里面有S3，S4，RC等形式，C++则是struct和class两种形式，举个例子

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

struct Date {
  unsigned int year;
  unsigned int month;
  unsigned int date;
};

class Date {
private:
  unsigned int year;
  unsigned int month;
  unsigned int date;
public:
  void setDate(int y, int m, int d);
  int getDay();
  int getMonth();
  int getYear();
}

强行联系的话，R的S3和C++的struct比较像，R的S4和C++的Date比较像，因为R的S3比较宽松，是基于list的堆叠，内部数据直接暴露到外部，而S4则是比较安全。C++的class用private保证private无法直接被外部访问，只能功过public暴露的函数进行操作。

指针和内存管理

指针和内存管理是两个比较高级的话题，如果你学习C语言不懂指针那你就和没学过C一样。不过我们学的是C++，只要不涉及到很高级的操作，那么我们完全可以避免接触它们。

总结

对于一个R语言使用者而言，如果想学习C++，那么至少需要知道如下几个概念

• 变量声明
• 控制结构，for, while, if-else
• 函数编写和调用
• 额外库
• 面向对象

基因组组装完成后，或者是完成了草图，就不可避免遇到一个问题，需要对基因组序列进行注释。注释之前首先得构建基因模型，有三种策略：

• 从头注释(de novo prediction)：通过已有的概率模型来预测基因结构，在预测剪切位点和UTR区准确性较低
• 同源预测(homology-based prediction)：有一些基因蛋白在相近物种间的保守型搞，所以可以使用已有的高质量近缘物种注释信息通过序列联配的方式确定外显子边界和剪切位点
• 基于转录组预测(transcriptome-based prediction)：通过物种的RNA-seq数据辅助注释，能够较为准确的确定剪切位点和外显子区域。

每一种方法都有自己的优缺点，所以最后需要用EvidenceModeler(EVM)和GLEAN工具进行整合，合并成完整的基因结构。基于可靠的基因结构，后续可才是功能注释，蛋白功能域注释，基因本体论注释，通路注释等。

那么基因注释重要吗？可以说是非常重要了，尤其是高通量测序非常便宜的现在。你可以花不到一万的价格对600M的物种进行100X的普通文库测序，然后拼接出草图。但是这个草图的价值还需要你进行注释后才能显现出来。有可能你和诺贝尔奖就差一个注释的基因组。

从案例中学习套路

陆地棉基因组注释

文章标题为“Sequencing of allotetraploid cotton (Gossypium hirsutum L. acc. TM-1) provides a resource for fiber improvement”.

同源注释：从Phytozome上下载了7个植物的基因组蛋白序列(Arabidopsis thaliana, Carica papaya, Glycine max, G. raimondii, Populus trichocarpa, Theobroma cacao and Vitis vinifera), 使用 TblastN 将蛋白序列比对到组装序列上，E-value的阈值为1e-5. 将不同蛋白的BLAST的hits用 Solar 软件进行合并。GeneWise 根据每个BLAST hit的对应基因区域预测完整的基因结构。

从头预测：先得构建repeat-mask genome， 在这个基础上就用 August, Genescan, GlimmerHMM, Geneid 和 SNAP 预测编码区

转录组预测：用Tophat将RNA-seq数据比对到组装序列上，然后用cufflinks组装转录本形成基因模型。

综上，使用 EvidenceModeler(EVM) 将上面的结果组装成非冗余的基因结构。进一步根据Cscore > 0.5，peptide coverage > 0.5 和CDS overlaping with TE进行筛选。还有过滤掉超过30%编码区被Pfam或Interprot TE domain的注释的基因模型。

这些基因模型使用BLASTP进行功能注释，所用数据库为SWiss-Prot和TrEMBL.蛋白功能使用InterProScan和HMMER注释，数据库为InterPro和Pfam。GO注释则是直接雇佣InterPro和Pfam注释得到的对应entry。通路注释使用KEGG数据库。

Cardamine hirsuta基因组注释

文章标题为“The Cardamine hirsuta genome offers insight into the evolution of morphological diversity”。

同源注释：使用 GenomeThreader 以拟南芥为剪切模型，以及PlantsGDB resourc上 Brassica rapa (v1.1), A. thaliana(TAIR10), A. lyrata (v6), tomato (v3.6), poplar (v2) 和 A. thaliana (version PUT-169), B. napus (version PUT-172) EST assemblies 的完整的代表性蛋白集。

转录本预测： 将 C. hirsuta RNA-seq数据比对到基因序列，然后用cufflinks拼接

从头预测：转录本预测得到的潜在蛋白编码转录本使用网页工具 ORFpredictor 进行预测， 同时用 blastx 和 A. thalina 进行比较，选择90%序列相似度和最高5%长度差异的部分从而保证保留完整的编码框(有启动子和终止子)。 这些基因模型根据相互之间的相似度和重叠度进行聚类，高度相似(>95)从聚类中剔除，保证非冗余训练集。为了训练gene finder, 它们选随机选取了2000个位点，20%是单个外显子基因。从头预测工具为 August , GlimmerHMM, Geneid 和 SNAP . 此外还用了Fgenesh+, 以双子叶特异矩阵为参数进行预测。

最后使用JIGSAW算法根据以上结果进行训练，随后再次用JIGSAW对每个基因模型计算统计学权重。

可变剪切模型则是基于苗、叶、花和果实的RNA-seq比对组装结果。

GO注释使用AHRD流程

小结

举的2个例子都是植物，主要是植物基因组不仅是组装，注释都是一大难题。因为植物基因组有大量的重复区，假基因，还有很多新的蛋白编码基因和非编码基因，比如说玉米基因组80%以上都是重复区域。然后当我检索这两篇文章所用工具的时候，我不经意或者说不可避免就遇到了这个网站 http://www.plantgdb.org/ , 一个整合植物基因组学工具和资源的网站，但是这个网站似乎2年没有更新了。当然这个网站也挺不错,http://bioservices.usd.edu/gsap.html, 他给出了一套完整的注释流程以及每一步的输入和输出情况。

此外，2017年在《Briefings in Bioinformatics》发表的”Plant genome and transcriptome annotations: from misconceptions to simple solution” 则是从五个角度对植物基因组注释做了很完整的总结

• 植物科学的常见本体
• 功能注释的常用数据库和资源
• 已注释的植物基因组意味着什么
• 一个自动化注释流程
• 一个参考流程图，用来说明使用公用数据库注释植物基因组/转录组的常规步骤

以上，通过套路我们对整个基因组注释有一个大概的了解，后续就需要通过实际操作来理解细节。

基因组注释

当我们谈到基因注释的时候，我们通常认为注释是指“对基因功能的描述”，比如说A基因在细胞的那个部分，通过招募B来调控C，从而引起病变。但是基因结构也是注释的一种形式，而且是先决条件，也就是在看似随机的ATCG的碱基排列中找到特殊的部分，而这些特殊的区域有着不一样的功能。

在正式启动基因组注释项目之前，需要先检查组装是否合格，比如contig N50的长度是否大于基因的平均长度，使用BUSCO/CEGMA检查基因的完整性，如果不满足要求，可能输出结果中大部分的contig中都不存在一个完整的基因结构。当组装得到的contig符合要求时，就可以开始基因组注释环节，这一步分为三步：基因结构预测，基因功能注释，可视化和质控。

基因组结构注释

基因结构注释应是功能注释的先决条件，完整的真核生物基因组注释流程需要如下步骤：

1. 必要的基因组重复序列屏蔽
2. 从头寻找基因, 可用工具为: GeneMarkHMM, FGENESH, Augustus, SNAP, GlimmerHMM, Genscan
3. 同源蛋白预测, 内含子分析: GeneWIse, Exonerate, GenomeThreader
4. 将EST序列，全长cDNA序列和Trinity/Cufflinks/Stringtie组装的转录组和基因组联配
5. 如果第4步用到了多个数据来源，使用PASA基于重叠情况进行联配
6. 使用EvidenceModler根据上述结果进行整合
7. 使用PASA更新EVM的一致性预测，增加UTR注释和可变剪切注释
8. 必要的人工检查

基本上是套路化的分析流程，也就有一些工具通过整合几步开发了流程管理工具，比如说BRAKER结合了GeneMark和Augustus，MAKER2整合了SNAP,Exonerate，虽然BRAKER说自己的效果比MAKER2好，但是用的人似乎不多，根据web of knowledge统计，两者的引用率分别是44,283, 当然BRAKER是2016，MAKER2是2011，后者在时间上有优势。

这里准备先按部就班的按照流程进行注释，所用的数据是 Cardamine hirsuta , 数据下载方式如下

1
2
3
4
5
6
7
8
9
10

# Cardamine hirsutat基因组数据
mkdir chi_annotation && cd chi_annotation
wget http://chi.mpipz.mpg.de/download/sequences/chi_v1.fa
cat chi_v1.fa | tr 'atcg' 'ATCG' > chi_unmasked.fa
# 注释结果
wget http://chi.mpipz.mpg.de/download/annotations/carhr38.gff
# Cardamine hirsutat转录组数据
mkdir rna-seq && cd rna-seq
wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/fruit_rnaseq/cardamine_hirsuta/ &
wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/leaf_rnaseq/cardamine_hirsuta/ &

软件安装不在正文中出现，会放在附录中，除了某些特别复杂的软件。

01-重复序列屏蔽

重复屏蔽：真核生物的基因组存在大量的重复序列，植物基因组的重复序列甚至可以高达80%。尽管重复序列对维持染色体的空间结构、基因的表达调控、遗传重组等都具有重要作用，但是却会导致BLAST的结果出现大量假阳性，增加基因结构的预测的计算压力甚至影响注释正确性。基因组中的重复按照序列特征可以分为两类：串联重复(tandem repeats)和散在重复(interspersed repeats).

鉴定基因组重复区域的方法有两种：一种基于文库(library)的同源(homology)方法，该文库收集了其他物种的某一种重复的一致性序列，通过相似性来鉴定重复；另一种是从头预测(de novo)，将序列和自己比较或者是高频K-mer来鉴定重复。

目前重复序列注释主要软件就是RepeatMasker和RepeatModel。这里要注意分析的fasta的ID不能过长，不然会报错。如果序列ID过长可以使用bioawk进行转换，后续用到RepatModel不支持多行存放序列的fasta格式。

直接使用同源注释工具RepeatMasker寻找重复序列：

1
2
3
4
5
6
7
8
9

mkdir 00-RepeatMask
~/opt/biosoft/RepeatMasker/RepeatMasker -e ncbi -species arabidopsis -pa 40 -gff -dir 00-RepeatMask/ chi_unmasked.fa
# -e ncbi
# -species 选择物种 用~/opt/biosoft/RepeatMasker/util/queryRepeatDatabase.pl -tree 了解
# -lib 增加额外数据库,
# -pa 并行计算
# -gff 输出gff注释
# -dir 输出路径
# annotation with the library produced by RepeatModel

输出结果中主要关注如下三个(其中xxx表示一类文件名)

• xxx.fa.masked, 将重复序列用N代替
• xxx.fa.out.gff, 以gff2形式存放重复序列出现的位置
• xxx.fa.tbl, 该文件记录着分类信息

1
2
3
4
5
6
7
8

cat 00-RepeatMask/chi_unmasked.fa.tbl
==================================================
file name: chi_unmasked.fa
sequences:           624
total length:  198654690 bp  (191241357 bp excl N/X-runs)
GC level:         35.24 %
bases masked:   35410625 bp ( 17.83 %)
==================================================

也就是说该物种198M中有将近18%的重复序列，作为参考，拟南芥125Mb 14%重复序列, 水稻389M，36%重复，人类基因组是3G，50%左右的重复序列。

使用最后的chi_unmasked.fa.masked用于下一步的基因结构预测。

注：当然也可以用RepeatModel进行从头预测，得到的预测结果后续可以整合到RepeatMasker

1
2
3

# de novo predict
~/opt/biosoft/RepeatModeler-open-1.0.11/BuildDatabase -name test -engine ncbi output.fa
~/opt/biosoft/RepeatModeler-open-1.0.11/RepeatModeler -database test

这一步速度极其慢，由于我们的目的只是获取屏蔽后序列降低后续从头预测的压力，所以可以先不做这一步。在后续分析重复序列在基因组进化上的作用时可以做这一步。下

如果从头预测的结果与同源预测的结果有30%以上的overlap，并且分类不一致，会把从头预测的结果过滤掉。从头预测与同源预测结果有overlap，但是分类一致的，都会保留。但是统计的时候不会重复统计。

02-从头(ab initio)预测基因

基于已有模型或无监督训练

目前的从头预测软件大多是基于HMM(隐马尔科夫链)和贝叶斯理论，通过已有物种的注释信息对软件进行训练，从训练结果中去推断一段基因序列中可能的结构，在这方面做的最好的工具是AUGUSTUS 它可以仅使用序列信息进行预测，也可以整合EST, cDNA, RNA-seq数据作为先验模型进行预测。

AUGUSTUS的无root安装比较麻烦，我折腾了好几天最后卒，不过辛亏有bioconda，conda create -n annotation augustus=3.3.

它的使用看起来很简单，我们可以尝试使用一段拟南芥已知的基因序列让其预测，比如前8k序列

1
2

seqkit faidx TAIR10.fa Chr1:1-8000 > test.fa
augustus --speices=arabidopsis test.fa > test.gff

如果仅仅看两者的CDS区，结果完全一致，相当于看过一遍参考答案去做题目，题目都做对了。

注:已经被训练的物种信息可以用augustus --species=help查看。

在不使用RNA-seq数据的情况下，可以基于拟南芥的训练模型进行预测，采用下面的方式多条染色体并行augustus

1
2
3
4
5
6

mkdir 01-augustsus && cd 01-augustsus
ln ../00-RepeatMask/chi_unmasked.fa.masked genome.fa
seqkit split genome.fa #结果文件在genome.fa.split
find genome.fa.split/ -type f -name "*.fa" | parallel -j 30 augustus --species=arabidopsis --gff3=on >> temp.gff #并行处理
join_aug_pred.pl < temp.gff  | grep -v '^#' > temp.joined.gff
bedtools sort -i temp.joined.gff > augustsus.gff

AUGUSTUS依赖于已有的模型，而GeneMark-ES/ET则是唯一一款支持无监督训练模型，之后再识别真核基因组蛋白编码区的工具。

1

gmes_petap.pl --ES --sequence genome.fa --cores 50

最后得到的是genemark.gtf，是标准的GTF格式，可以使用Sequence Ontology Project提供的gtf2gff3.pl进行转换

1
2
3

wget http://genes.mit.edu/burgelab/miso/scripts/gtf2gff3.pl
chmod 755 gtf2gff3.pl
gtf2gff3.pl genemark.gtf | bedtools sort -i - > genemark.gff

不同从头预测软件的实际效果可以通过在IGV中加载文章提供的gff文件和预测后的gff文件进行比较，一般会存在如下几个问题：

• 基因多了，或者少了，也就是假阳性和假阴性现象
• UTR区域难以预测，这个比较正常
• 未正确识别可变剪切位点，导致前后几个基因识别成一个基因

考虑到转录组测序已经非常便宜，可以通过该物种的RNA-seq提供覆盖度信息进行预测。

基于转录组数据预测

根据已有的模型或者自训练可以正确预测很大一部分的基因，但如果需要提高预测的正确性，还需要额外的信息。在过去就需要提供物种本身的cDNA, EST，而现在更多的是基于转录组序列进行训练。尽管RNA-seq数据在基因组上的比对情况能够推测出内含子位置，根据覆盖度可以推测出外显子和非编码区的边界，但是仅仅依赖于RNA-seq的覆盖不能可信地推测出蛋白编码区(Hoff K.J. Stanke M. 2015).

AUGUSTUS可以利用转录组比对数据中的位置信息来训练模型，GeneMark-ET可以利用RNA-seq得到的内含子位点信息自我训练HMM参数，进行基因预测。BRAKER2将两者进行整合，使用GeneMark-ET根据RNA-seq无监督训练模型寻找基因，然后用AUGUSTUS进行模型训练，最后完成基因预测

首先使用hisat2根据屏蔽后的参考序列建立索引，进行比对。

1
2
3
4
5
6

# 项目根目录
mkdir index
hisat2-build 01-augustus/genome.fa index/chi_masked
hisat2 -p 20 -x index/chi_masked -1 rna-seq/leaf_ox_r1_1.fastq.gz -2 rna-seq/leaf_ox_r1_2.fastq.gz | samtools sort -@ 10 > 02-barker/leaf_ox_r1.bam &
isat2 -p 20 -x index/chi_masked -1 rna-seq/ox_flower9_rep1_1.fastq.gz -2 rna-seq/ox_flower9_rep1_2.fastq.gz | samtools sort -@ 10 > 02-barker/ox_flower9.bam &
hisat2 -p 20 -x index/chi_masked -1 rna-seq/ox_flower16_rep1_1.fastq.gz -2 rna-seq/ox_flower16_rep1_2.fastq.gz | samtools sort -@ 10 > 02-barker/ox_flower16.bam &

然后，以未屏蔽重复序列的参考序列和BAM文件作为输入，让BRAKER2（安装会稍显麻烦，因为依赖许多软件）进行预测。

1
2
3
4
5

braker.pl --gff3 --cores 50 --species=carhr --genome=chi_unmasked.fa --bam=02-barker/leaf_ox_r1.bam,02-barker/ox_flower16.bam,02-barker/ox_flower9.bam
# --gff3: 输出GFF3格式
# --genome: 基因组序列
# --bam: 比对后的BAM文件，允许多个
# --cores: 处理核心数

最后会得到如下输出文件

• hintsfile.gff: 从RNA-seq比对结果的BAM文件中提取，其中内含子用于训练GeneMark-EX, 使用所有特征训练AUGUSTUS
• GeneMark-ET/genemark.gtf: GeneMark-EX根据RNA-seq数据训练后预测的基因
• augustus.hints.gff: AUGUSTUS输出文件

将augustus.hints.gff3和文章的注释文件(carhr38.gtf)比较，见下图：

其实不难发现，在不考虑UTR区域情况下，两者的差别其实更多表现是基因数目上，其实也就是利用转录组数据推测结构的问题所在，没有覆盖的区域到底是真的没有基因，还是有基因结构只不过所用组织没有表达，或者说那个区域其实是假基因？此外，如果基因间隔区域很短，有时候还会错误地把两个不同的基因预测为一个基因。因此，应该注重RNA-seq数据在剪切位点识别和外显子边界确定的优势。

03-同源预测基因结构

同源预测(homology prediction)利用近缘物种已知基因进行序列比对，找到同源序列。然后在同源序列的基础上，根据基因信号如剪切信号、基因起始和终止密码子对基因结构进行预测，如下示意图：

相对于从头预测的“大海捞针”，同源预测相当于先用一块磁铁在基因组大海中缩小了可能区域，然后从可能区域中鉴定基因结构。在10年之前，当时RNA-seq还没有普及, 只有少部分物种才有EST序列和cDNA序列的情况下，这的确是一个比较好的策略，那么问题来了，现在还需要进行这一步吗，如果需要是出于那种角度考虑呢?

在同源预测上，目前看到的大部分基因组文章都是基于TBLASTN + GeneWise，这可能是因为大部分基因组文章都是国内做的，这些注释自然而言用的就是公司的流程，然后目前国内的公司大多数又和某一家公司有一些关系。不过最近的3010水稻泛基因组用的是MAKER, 感谢部分提到这部分工作是由M. Roa(Philippine Genome Center Core Facilities for Bioinformatics, Department of Science)做的，算是一股清流吧。当然我在看Cardamine hirsuta基因组注释文章，发现它们同源注释部分用的是GenomeThreader, 该工具在本篇文章成文时的3月之前又更新了。

GeneWise的网站说它目前由Ewan Birney维护，只不过不继续开发了，因为Guy Slater开发Exonerate解决了GeneWise存在的很多问题，并且速度快了1000倍。考虑到目前只有GeneWise能利用HMM根据蛋白找DNA，而且ENSEMBL的注释流程也有一些核心模块用到了它，所以作者依旧在缓慢的开发这个工具(自2.4.1已经10多年没有更新了)，当然这个工具也是非常的慢。尽管这一步不会用到GeneWise作为我们的同源注释选项，但是我们可以尝试用GeneWise手工注释一个基因，主要步骤如下

• 第一步： 使用BLASTX，根据dna序列搜索到蛋白序列，只需要第一个最佳比对结果
• 第二步： 选择最佳比对的氨基酸序列
• 第三步： 将dna序列前后延长2kb，与氨基酸序列一并传入给genewise进行同源预测

提取前5K序列，然后选择在TAIR上用BLASTX进行比对

1

seqkit faidx chi_unmasked.fa Chr1:1-5000 > chr1_5k.fa

选择第一个比对结果中的氨基酸序列，和前5k的DNA序列一并作为GeneWise的输入

最后的结果出乎了我的意料

让我们跳过这个尴尬的环节，毕竟很可能是我不太熟练使用工作所致。这里说点我的看法，除非你真的没有转录组数据，必须要用到同源物种的蛋白进行预测，或者你手动处理几个基因，否则不建议使用这个工具，因为你可能连安装都搞不定。

让我们用GenomeThreader基于上面的DNA序列和氨基酸序列进行同源基因结构预测吧

1
2

gth -genomic chr1_5k.fa -protein cer.fa -intermediate -gff3out
# 其中cer.fa就是AT1G02205.2的氨基酸序列

结果一致，并且从RNA-seq的覆盖情况也符合预期

1
2
3
4
5
6
7
8
9
10

Chr1	gth	exon	1027	1197	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    1027    1197
Chr1	gth	exon	1275	1448	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    1275    1448
Chr1	gth	exon	1541	1662	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    1555    1662
Chr1	gth	exon	1807	2007	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    1807    2007
Chr1	gth	exon	2085	2192	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    2085    2192
Chr1	gth	exon	2294	2669	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    2294    2669
Chr1	gth	exon	3636	3855	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    3636    3855
Chr1	gth	exon	3971	4203	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    3971    4203
Chr1	gth	exon	4325	4548	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    4325    4548
Chr1	gth	exon	4676	4735	Parent=gene1	Chr1    MIPS_CARH_v3.8  exon    4676    4735

全基因组范围预测流程如下：

准备cDNA和或protein序列：在https://phytozome.jgi.doe.gov/p下载靠谱的物种的蛋白质序列，如 Arabidopsis thaliana, Oryza sativa, Brassica rapa, 查找文献寻找目前该物种的已有EST/cDNA序列，或者RNA-seq从头组装转录组。这里仅考虑用同源物种的蛋白序列进行比对分析，转录组从头组装数据用于PASA整体比对到参考基因组和更新已有的基因解雇。

分别测试下不同物种的同源注释结果

1
2
3
4

#run seperately
gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/Athaliana_167_TAIR10.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Athaliana.gff3 &
gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/BrapaFPsc_277_v1.3.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Brapa.gff3 &
gth -species arabidopsis -translationtable 1 -gff3 -intermediate -protein ~/db/protein_db/Osativa_323_v7.0.protein.fa.gz -genomic chi_unmasked.fa -o 03-genomethreader/Osativa.gff3 &

在定性角度上来看，同源注释的结果和从头预测的没啥差别, 其中B. rapa和A. thaliana和C. hirsuta都属于十字花科，而O. sativa是禾本科, 所以前两者预测的效果好。

当然实际的同源注释流程中不能是单个物种分别预测，应该是将所有的蛋白序列进行合并，然后用BLASTX找到最优的联配，之后用GenomeThreader进行预测。PASA流程提到的UniRef90作为同源注释的搜索数据库可能是更好的选择，由于UniRef优先选择哪些人工审查、注释质量高、来源于模式动植物的蛋白，所以可靠性相对于直接使用同源物中可能更高。

BLASTX + GenomeThreader的代码探索中

04-RNA-seq的两种使用策略

对于RNA-seq数据，有两种使用策略，一种是使用HISAT2 + StringTie先比对再组装, 一种是从头组装，然后使用PASA将转录本比对到基因组上。

基于HISAT2 + StringTie

首先，使用HISAT2将RNA-seq数据比对到参考基因组, 这一步和之前相似，但是要增加一个参数--dta，使得StingTie能更好的利用双端信息

1
2
3
4
5

hisat2-build 01-augustus/genome.fa index/chi_masked
hisat2 --dta -p 20 -x index/chi_masked -1 rna-seq/leaf_ox_r1_1.fastq.gz -2 rna-seq/leaf_ox_r1_2.fastq.gz | samtools sort -@ 10 > rna-seq/leaf_ox_r1.bam &
hisat2 --dta -p 20 -x index/chi_masked -1 rna-seq/ox_flower9_rep1_1.fastq.gz -2 rna-seq/ox_flower9_rep1_2.fastq.gz | samtools sort -@ 10 > rna-seq/ox_flower9.bam &
hisat2 --dta -p 20 -x index/chi_masked -1 rna-seq/ox_flower16_rep1_1.fastq.gz -2 rna-seq/ox_flower16_rep1_2.fastq.gz | samtools sort -@ 10 > rna-seq/ox_flower16.bam &
samtools merge -@ 10 rna-seq/merged.bam rna-seq/leaf_ox_r1.bam rna-seq/ox_flower9.bam rna-seq/ox_flower16.bam

然后用StringTie进行转录本预测

1

stringtie -p 10 -o rna-seq/merged.gtf rna-seq/merged.bam

对于后续的EvidenceModeler而言，它不需要UTR信息，只需要编码区CDS，需要用TransDecoder进行编码区预测

1
2
3
4
5
6
7
8

util/cufflinks_gtf_genome_to_cdna_fasta.pl merged.gtf input/chi_masked.fa > transcripts.fasta
util/cufflinks_gtf_to_alignment_gff3.pl merged.gtf > transcripts.gff3
TransDecoder.LongOrfs -t transcripts.fasta
TransDecoder.Predict -t transcripts.fasta
util/cdna_alignment_orf_to_genome_orf.pl \
     transcripts.fasta.transdecoder.gff3 \
     transcripts.gff3 \
     transcripts.fasta > transcripts.fasta.transdecoder.genome.gff3

最后结果transcripts.fasta.transdecoder.gff3用于提供给EvidenceModeler

基于PASA

在多年以前，那个基因组组装还没有白菜价，只有几个模式物种基因组的时代，对于一个未测序的基因组，研究者如果要研究某一个基因的功能，大多会通过同源物种相似基因设计PCR引物，然后去扩增cDNA. 如果是一个已知基因组的物种，如果要大规模识别基因, 研究者通常会使用EST(expressed sequence tags)序列。

相对于基于算法的从头预测，cDNA和EST序列更能够真实的反应出一个基因的真实结构，如可变剪切、UTR和Poly-A位点。PASA(Progam to Assemble Spliced Alignments)流程最早用于拟南芥基因组注释，最初的设计是通过将全长(full-length)cDNA和EST比对到参考基因组上，去发现和更新基因组注释。其中FL-cDNA和EST序列对最后结果的权重不同。

这是以前的故事，现在的故事是二代转录组以及一些三代转录组数据，那么如何处理这些数据呢？我认为三代转录组相对于过去的FL-cDNA，而二代转录组数据经过拼接后可以看作是更长的EST序列。由于目前最普及的还是普通的mRNA-seq, 也就只介绍这部分流程。

考虑到我还没有研究过三代的全长转录组，分析过数据，这里的思考极有可能出错，后续可能会修改这一部分思考。

转录组组装使用Trinity(conda安装)

1
2

cd rna-seq
Trinity --seqType fq --CPU 50 --max_memory 64G --left leaf_ox_r1_1.fastq.gz,ox_flower16_rep1_1.fastq.gz,ox_flower9_rep1_1.fastq.gz --right leaf_ox_r1_2.fastq.gz,ox_flower16_rep1_2.fastq.gz,ox_flower9_rep1_2.fastq.gz &

PASA是由30多个命令组成的流程，相关命令位于PASApipeline/scripts，为了适应不同的分析，有些参数需要通过修改配置文件更改,

1
2
3
4
5

cp ~/opt/biosoft/PASApipeline/pasa_conf/pasa.alignAssembly.Template.txt alignAssembly.config
# 修改如下内容
DATABASE=database.sqlite
validate_alignments_in_db.dbi:--MIN_PERCENT_ALIGNED=80
validate_alignments_in_db.dbi:--MIN_AVG_PER_ID=80

上述几行配置文件表明SQLite3数据库的名字，设置了scripts/validate_alignments_in_db.dbi的几个参数, 表示联配程度和相似程度。后续以Trinity组装结果和参考基因组作为输入，运行程序：

1

~/opt/biosoft/PASApipeline/scripts/Launch_PASA_pipeline.pl -c alignAssembly.config -C -R -g ../chi_unmasked.fa -t ../rna-seq/trinity_out_dir/Trinity.fasta --ALIGNERS blat,gmap

最后结果如下：

• database.sqlite.pasa_assemblies_described.txt
• database.sqlite.pasa_assemblies.gff3
• database.sqlite.pasa_assemblies.gtf
• database.sqlite.pasa_assemblies.bed

其中gff3格式用于后续的分析。

目前的一些想法， 将从头组装的转录本比对到参考基因组上很大依赖组装结果，所以和EST序列和cDNA相比，质量上还有一点差距。

05-整合预测结果

从头预测，同源注释和转录组整合都会得到一个预测结果，相当于收集了大量证据，下一步就是通过这些证据定义出更加可靠的基因结构，这一步可以通过人工排查，也可以使用EVidenceModeler(EVM). EVM只接受三类输入文件：

• gene_prediction.gff3: 标准的GFF3格式，必须要有gene, mRNA, exon, CDS这些特征，用EVidenceModeler-1.1.1/EvmUtils/gff3_gene_prediction_file_validator.pl验证
• protein_alignments.gff3: 标准的GFF3格式，第9列要有ID信和和target信息, 标明是比对结果
• transcript_alignments.gff3:标准的GFF3格式，第9列要有ID信和和target信息，标明是比对结果

EVM对gene_prediction.gff3有特殊的要求，就是GFF文件需要反映出一个基因的结构，gene->(mRNA -> (exon->cds(?))(+))(+), 表示一个基因可以有多个mRNA，即基因的可变剪接, 一个mRNA都可以由一个或者多个exon(外显子), 外显子可以是非翻译区(UTR),也可以是编码区(CDS). 而GlimmerHMM, SNAP等

这三类根据人为经验来确定其可信度，从直觉上就是用PASA根据mRNA得到的结果高于从头预测。

第一步：创建权重文件,第一列是来源类型(ABINITIO_PREDICTION, PROTEIN, TRANSCRIPT), 第二列对应着GFF3文件的第二列，第三列则是权重.我这里用了三个来源的数据。

1
2
3
4
5

mkdir 05-EVM && cd 05-EVM
#vim weights.txt
ABINITIO_PREDICTION      augustus       4
TRANSCRIPT      assembler-database.sqlite      7
OTHER_PREDICTION  transdecoder  8

我觉得根据基因组引导组装的ORF的可信度高于组装后比对，所以得分和PASA差不多一样高。从头预测权重一般都是1，但是BRAKER可信度稍微高一点，可以在2~5之间。

第二步：分割原始数据, 用于后续并行. 为了降低内存消耗，–segmentsSize设置的大小需要少于1Mb(这里是100k)， –overlapSize的不能太小，如果数学好，可用设置成基因平均长度加上2个标准差，数学不好，就设置成10K吧

1
2
3

cat transcripts.fasta.transdecoder.genome.gff3 ../braker/carhr/augustus.hints.gff3 > gene_predictions.gff3
ln ../04-align-transcript/database.sqlite.pasa_assemblies.gff3 transcript_alignments.gff3
~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/partition_EVM_inputs.pl --genome ../chi_unmasked.fa --gene_predictions gene_predictions.gff3 --transcript_alignments transcript_alignments.gff3 --segmentSize 100000 --overlapSize 10000 --partition_listing partitions_list.out

第三步：创建并行运算命令并且执行

1
2
3
4
5

~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/write_EVM_commands.pl --genome ../chi_unmasked.fa --weights `pwd`/weights.txt \
      --gene_predictions gene_predictions.gff3 \
      --transcript_alignments transcript_alignments.gff3 \
      --output_file_name evm.out  --partitions partitions_list.out >  commands.list
parallel --jobs 10 < commands.list

第四步：合并并行结果

1

~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/recombine_EVM_partial_outputs.pl --partitions partitions_list.out --output_file_name evm.out

第五步：结果转换成GFF3

1
2

~/opt/biosoft/EVidenceModeler-1.1.1/EvmUtils/convert_EVM_outputs_to_GFF3.pl  --partitions partitions_list.out --output evm.out  --genome ../chi_unmasked.fa
find . -regex ".*evm.out.gff3" -exec cat {} \; | bedtools sort -i - > EVM.all.gff

当前权重设置下，EVM的结果更加严格，需要按照实际情况调整，增加其他证据。

06-可选步骤

注释过滤：对于初步预测得到的基因，还可以稍微优化一下，例如剔除编码少于50个AA的预测结果，将转座子单独放到一个文件中(软件有TransposonPSI)。

这里基于gffread先根据注释信息提取所有的CDS序列，过滤出长度不足50AA的序列，基于这些序列过滤原来的的注释

1
2
3

gffread EVM.all.gff -g input/genome.fa -y tr_cds.fa
bioawk -c fastx '$seq < 50 {print $comment}' tr_cds.fa | cut -d '=' -f 2 > short_aa_gene_list.txt
grep -v -w -f short_aa_gene_list.txt EvM.all.gff > filter.gff

使用PASA更新EVM结果：EVM结果不包括UTR区域和可变剪切的注释信息，可以使用PASA进行更新。然而这部分已经无法逃避MySQL, 服务器上并没有MySQL的权限，我需要学习Perl脚本进行修改。因此基因结构注释到此先放一放。

07-基因编号

对每个基因实现编号，形如ABCD000010的效果，方便后续分析。如下代码是基于EVM.all.gff，使用方法为python gffrename.py EVM_output.gff prefix > renamed.gff.

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

#!/usr/bin/env python3
import re
import sys

if len(sys.argv) < 3:
    sys.exit()

gff = open(sys.argv[1])
prf = sys.argv[2]

count = 0
mRNA  = 0
cds   = 0
exon  = 0

print("##gff-version 3.2.1")
for line in gff:
    if not line.startswith("\n"):
        records = line.split("\t")
        records[1] = "."
    if re.search(r"\tgene\t", line):
        count = count + 10
        mRNA  = 0
        gene_id = prf + str(count).zfill(6)
        records[8] = "ID={}".format(gene_id)
    elif re.search(r"\tmRNA\t", line):
        cds   = 0
        exon  = 0
        mRNA  = mRNA + 1
        mRNA_id    = gene_id + "." + str(mRNA)
        records[8] = "ID={};Parent={}".format(mRNA_id, gene_id)
    elif re.search(r"\texon\t", line):
        exon     = exon + 1
        exon_id  = mRNA_id + "_exon_" + str(exon)
        records[8] = "ID={};Parent={}".format(exon_id, mRNA_id)
    elif re.search(r"\tCDS\t", line):
        cds     = cds + 1
        cds_id  = mRNA_id + "_cds_" + str(cds)
        records[8] = "ID={};Parent={}".format(cds_id, mRNA_id)
    else:
        continue

    print("\t".join(records))

gff.close()

一些经验

如果有转录组数据，没必须要使用太多的从头预测工具，braker2 加 GlimmerHMM可能就够用了, 更多是使用PASA和StringTie利用好转录组数据进行注释。

基因功能注释

基因功能的注释依赖于上一步的基因结构预测，根据预测结果从基因组上提取翻译后的 蛋白序列 和主流的数据库进行比对，完成功能注释。常用数据库一共有以几种：

• Nr：NCBI官方非冗余蛋白数据库，包括PDB, Swiss-Prot, PIR, PRF; 如果要用DNA序列，就是nt库
• Pfam: 蛋白结构域注释的分类系统
• Swiss-Prot: 高质量的蛋白数据库，蛋白序列得到实验的验证
• KEGG: 代谢通路注释数据库.
• GO: 基因本体论注释数据库

除了以上几个比较通用的数据库外，其实还有很多小众数据库，应该根据课题研究和背景进行选择。注意，数据库本身并不能进行注释，你只是通过序列相似性进行搜索，而返回的结果你称之为注释。因此数据库和搜索工具要进行区分，所以你需要单独下载数据库和搜索工具，或者是同时下载包含数据库和搜索工具的安装包。

注意，后续分析中一定要保证你的蛋白序列中不能有代表氨基酸字符以外的字符，比如说有些软件会把最后一个终止密码子翻译成”.”或者”*“

BLASTP

这一部分用到的数据库都是用BLASTP进行检索，基本都是四步发：下载数据库，构建BLASTP索引，数据库检索，结果整理。其中结果整理需要根据BLASTP的输出格式调整。

Nr的NCBI收集的最全的蛋白序列数据库，但是无论是用NCBI的BLAST还是用速度比较快DIAMOND对nr进行搜索，其实都没有利用好物种本身的信息。因此在RefSeq上下载对应物种的蛋白序列, 用BLASTP进行注释即可。

1
2
3
4
5
6
7
8

# download
wget -4  -nd -np -r 1 -A *.faa.gz ftp://ftp.ncbi.nlm.nih.gov/refseq/release/plant/
mkdir -p ~/db/RefSeq
zcat *.gz > ~/db/RefSeq/plant.protein.faa
# build index
~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in plant.protein.faa -dbtype prot -parse_seqids -title RefSeq_plant -out plant
# search
~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out RefSeq_plant_blastp.xml -db ~/db/RefSeq/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 &

Swiss-Prot里收集了目前可信度最高的蛋白序列，一共有55w条记录，数据量比较小，

1
2
3
4
5
6
7

# download
wget -4 -q ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gzip -d uniprot_sprot.fasta.gz
# builid index
~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in uniprot_sprot.fasta -dbtype prot -title swiss_prot -parse_seqids
# search
~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out swiss_prot.xml -db ~/db/swiss_prot/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 &

关于结果整理，已经有很多人写了脚本，比如说我搜索BLAST XML CSV，就找到了https://github.com/Sunhh/NGS_data_processing/blob/master/annot_tools/blast_xml_parse.py, 所以就不过多介绍。

InterProScan

下面介绍的工具是InterProScan, 从它的9G的体量就可以感受它的强大之处，一次运行同时实现多个信息注释。

• InterPro注释
• Pfam数据库注释(可以通过hmmscan搜索pfam数据库完成)
• GO注释(可以基于NR和Pfam等数据库，然后BLAST2GO完成,)
• Reactome通路注释，不同于KEGG

命令如下

1

./interproscan-5.29-68.0/interproscan.sh -appl Pfam  -f TSV -i sample.fa -cpu 50 -b sample -goterms -iprlookup -pa

-appl告诉软件要执行哪些数据分析，勾选的越多，分析速度越慢，Pfam就行。

KEGG

KEGG数据库目前本地版收费，在线版收费，所以只能将蛋白序列在KEGG服务器上运行。因此你需要在http://www.genome.jp/tools/kaas/选择合适的工具进行后续的分析。我上传的50M大小蛋白序列，在KEGG服务器上只需要运行8个小时，也就是晚上提交任务，白天回来干活。

附录

基因组注释的常用软件：

• 重复区域
  • RepeatMasker：识别基因组中的可能重复
  • RepeatModeler: 识别新的重复序列
  • LTR-FINDER: http://tlife.fudan.edu.cn/ltr_finder/
• 从头预测
  • Augustus
  • Fgenesh
• 同源预测
  • GeneWise
  • Exonerate
  • Trinity
  • GenomeThreader
• 注释合并
  • GLEAN：已经落伍于时代了
  • EvidenceModeler： 与时俱进
• 流程
  • PASA：真核生物基因的转录本可变剪切自动化注释项目，需要提供物种的EST或RNA-seq数据
  • MAKER
  • BRAKER1: 使用GeneMark-ET和AUGUSTUS基于RNA-Seq注释基因结构
  • EuGene
• 可视化
  • IGV
  • JBrowse/GBrowse

参考文献和推荐阅读：

• NCBI真核生物基因组注释流程https://www.ncbi.nlm.nih.gov/genome/annotation_euk/process/
• 真核基因组注释入门: “A beginner’s guide to eukaryotic genome annotation”
• 二代测序注释流程:Comparative Gene Finding: “Annotation Pipelines for Next-Generation Sequencing Projects”
• 基因组转录组注释策略: “Plant genome and transcriptome annotations: from misconceptions to simple solution”
• 重复序列综述: “Repetitive DNA and next-generation sequencing: computational challenges and solutions”
• MAKER2教程: http://weatherby.genetics.utah.edu/MAKER/wiki/index.php/MAKER_Tutorial_for_WGS_Assembly_and_Annotation_Winter_School_2018
• 《生物信息学》 樊龙江: 第1-5章: 基因预测与功能注释
• 《NGS生物信息分析》 陈连福： 真核生物基因组基因注释
• JGS流程: https://genome.jgi.doe.gov/programs/fungi/FungalGenomeAnnotationSOP.pdf

环境准备

数据下载

1
2
3
4
5
6
7

# Cardamine hirsutat基因组数据
mkdir chi_annotation && cd chi_annotation
wget http://chi.mpipz.mpg.de/download/sequences/chi_v1.fa
cat chi_v1.fa | tr 'atcg' 'ATCG' > chi_unmasked.fa
# Cardamine hirsutat转录组数据
wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/fruit_rnaseq/cardamine_hirsuta/ &
wget -4 -q -A '*.fastq.gz' -np -nd -r 2 http://chi.mpipz.mpg.de/download/leaf_rnaseq/cardamine_hirsuta/ &

软件安装

RepeatMasker: 用于注释基因组的重复区，需要安装RMBlast, TRF，以及在http://www.girinst.org注册以下载Repbase

安装RepeatMasker

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

cd ~/src
wget http://tandem.bu.edu/trf/downloadstrf409.linux64
mv trf409.linux64 ~/opt/bin/trf
chmod a+x ~/opt/bin/trf
# RMBlast
cd ~/src
wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.6.0/ncbi-blast-2.6.0+-src.tar.gz
wget http://www.repeatmasker.org/isb-2.6.0+-changes-vers2.patch.gz
tar xf ncbi-blast-2.6.0+-src
gunzip isb-2.6.0+-changes-vers2.patch.gz
cd ncbi-blast-2.6.0+-src
patch -p1 < ../isb-2.6.0+-changes-vers2.patch
cd c++
./configure --with-mt --prefix=~/opt/biosoft/rmblast --without-debug && make && make install
# RepeatMasker
cd ~/src
wget http://repeatmasker.org/RepeatMasker-open-4-0-7.tar.gz
tar xf RepeatMasker-open-4-0-7.tar.gz
mv RepeatMasker ~/opt/biosoft/
cd ~/opt/biosoft/RepeatMasker
## 解压repbase数据到Libraries下
## 配置RepatMasker
perl ./configure

在上面的基础上安装RepeatModel

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# RECON
cd ~/src
wget -4 http://repeatmasker.org/RepeatModeler/RECON-1.08.tar.gz
tar xf RECON-1.08.tar.gz
cd RECON-1.08/src
make && make install
cd ~/src
mv RECON-1.08 ~/opt/biosoft
# nesg
cd ~/src
mkdir nesg && cd nesg
wget -4 ftp://ftp.ncbi.nih.gov/pub/seg/nseg/*
make
mv nmerge nseg ~/opt/bin/
# RepeatScout
http://www.repeatmasker.org/RepeatScout-1.0.5.tar.gz
# RepeatModel
wget -4 http://repeatmasker.org/RepeatModeler/RepeatModeler-open-1.0.11.tar.gz
tar xf RepeatModeler-open-1.0.11.tar.gz
mv RepeatModeler-open-1.0.11 ~/opt/biosoft/
cd ~/opt/biosoft/RepeatModeler-open-1.0.11
# 配置
perl ./configure
export PATH=~/opt/biosoft/maker:$PATH

BLAST，BLAST有两个版本可供选择, WuBLAST或者NCBI-BLAST，我个人倾向于NCBI-BLAST，并且推荐使用编译后二进制版本，因为编译实在是太花时间了

1
2
3
4
5
6
7

cd ~/src
wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.7.1+-x64-linux.tar.gz
tar xf ncbi-blast-2.7.1+-x64-linux.tar.gz -C ~/opt/biosoft
# 环境变量
export PATH=~/opt/biosoft/ncbi-blast-2.7.1+/bin:$PATH
# 用于后续的BRAKER2
conda create -n annotation blast=2.2.31

AUGUSTUS: 可以说是最好的预测软件，使用conda安装

1
2

source activate annotation
conda install augustus=3.3

GeneMark-ES/ET则是唯一一款支持无监督训练模型, 软件下载需要登记

1
2
3
4
5
6
7
8

cd ~/src
wget http://topaz.gatech.edu/GeneMark/tmp/GMtool_Qg87n/gm_et_linux_64.tar.gz
tar xf gm_et_linux_64.tar.gz
mv gm_et_linux_64/gmes_petap/ ~/opt/biosoft
wget http://topaz.gatech.edu/GeneMark/tmp/GMtool_Qg87n/gm_key_64.gz
gzip -dc gm_key_64.gz > ~/.gm_key
cpan YAML Hash::Merge Logger::Simple Parallel::ForkManager
echo "export PATH=$PATH:~/opt/biosoft/gmes_petap/" >> ~/.bashrc

GlimmerHMM:

1
2
3

cd ~/src
wget -4 ftp://ccb.jhu.edu/pub/software/glimmerhmm/GlimmerHMM-3.0.4.tar.gz
tar xf GlimmerHMM-3.0.4.tar.gz -C ~/opt/biosoft

SNAP: 基因从头预测工具，在处理含有长内含子上的基因组上表现欠佳

1
2
3
4
5
6
7
8
9
10

# 安装
cd ~/src
git clone https://github.com/KorfLab/SNAP.git
cd SNP
make
cd ..
mv SNAP ~/opt/biosoft
# 环境变量
export Zoe=~/opt/biosoft/SNAP/Zoe
export PATH=~/opt/biosoft/SNAP:$PATH

Exnerate 2.2: 配对序列比对工具，提供二进制版本, 功能类似于GeneWise，能够将cDNA或蛋白以gao align的方式和基因组序列联配。

1
2
3
4
5
6
7
8

cd ~/src
wget http://ftp.ebi.ac.uk/pub/software/vertebrategenomics/exonerate/exonerate-2.2.0-x86_64.tar.gz
tar xf exonerate-2.2.0-x86_64.tar.gz
mv exonerate-2.2.0-x86_64 ~/opt/biosoft/exonerate-2.2.0
# .bashrc添加环境变量
export PATH=~/opt/biosoft/exonerate-2.2.0:$PATH
# 或
conda install -c bioconda exonerate

GenomeThreader 1.70: 同源预测软件，1.7.0版本更新于2018年2月

1
2
3
4
5
6

wget -4 http://genomethreader.org/distributions/gth-1.7.0-Linux_x86_64-64bit.tar.gz
tar xf gth-1.7.0-Linux_x86_64-64bit.tar.gz -C ~/opt/biosoft
# 修改.bashrc增加如下行
export PATH=$PATH:$HOME/opt/biosoft/gth-1.7.0-Linux_x86_64-64bit/bin
export BSSMDIR="$HOME/opt/biosoft/gth-1.7.0-Linux_x86_64-64bit/bin/bssm"
export GTHATADIR="$HOME/opt/biosoft/gth-1.7.0-Linux_x86_64-64bit/bin/gthdata"

BRAKER2: 依赖AUGUSTUS 3.3, GeneMark-EX 4.33, BAMTOOLS 2.5.1, NCBI BLAST+ 2.2.31+(可选 SAMTOOLS 1.74+, GenomeThreader 1.70)

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

cpan File::Spec::Functions Module::Load::Conditional POSIX Scalar::Util::Numeric YAML File::Which Logger::Simple Parallel::ForkManager
cd ~/src
wget -4 http://exon.biology.gatech.edu/GeneMark/Braker/BRAKER2.tar.gz
tar xf BRAKER2.tar.gz -C ~/opt/biosoft
echo "export PATH=$PATH:$HOME/opt/biosoft/BRAKER_v2.1.0/" >> ~/.bashrc
# 在~/.bashrc设置如下软件所在环境变量
export AUGUSTUS_CONFIG_PATH=$HOME/miniconda3/envs/annotation/config/
export AUGUSTUS_SCRIPTS_PATH=$HOME/miniconda3/envs/annotation/bin/
export BAMTOOLS_PATH=$HOME/miniconda3/envs/annotation/bin/
export GENEMARK_PATH=$HOME/opt/biosoft/gmes_petap/
export SAMTOOLS_PATH=$HOME/miniconda3/envs/annotation/bin/
export ALIGNMENT_TOOL_PATH=$HOME/opt/biosoft/gth-1.7.0-Linux_x86_64-64bit/bin/

TransDecoder 编码区域预测工具，需要预先安装NCBI-BLAST

1
2
3
4
5
6

cpan URI::Escape
cd ~/src
wget -4 https://github.com/TransDecoder/TransDecoder/archive/TransDecoder-v5.3.0.zip
unzip TransDecoder-v5.3.0.zip
cd TransDecoder-v5.3.0
make test

MARKER: 使用conda安装会特别的方便，最好新建环境

1

conda create -n marker marker

PASA: 依赖于一个数据库(MySQL或SQLite), Perl模块(DBD::mysql或DBD::SQLite), GMAP, BLAT, Fasta3。由于MySQL在HPC集群中的表现不如SQLite，以及安装MySQL还需要各种管理员权限，于是就有人进行了修改，增加了feature/sqlite分支, 见Add support for SQLite

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

cpan DB_File URI::Escape DBI DBD::SQLite
# GMAP
wget http://research-pub.gene.com/gmap/src/gmap-gsnap-2017-11-15.tar.gz
tar xf gmap-gsnap-2017-11-15.tar.gz
cd gmap-2017-11-15
./configure --prefix=$HOME/opt/gmap
make && make install
# BLAT
cd ~/opt/bin
wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/blat/blat && chmod 755 ./blat
# Fasta3
wget -4 http://faculty.virginia.edu/wrpearson/fasta/fasta36/fasta-36.3.8g.tar.gz && \
        tar zxvf fasta-36.3.8g.tar.gz && \
        cd ./fasta-36.3.8g/src && \
        make -f ../make/Makefile.linux_sse2 all && \
        cp ../bin/fasta36 ~/opt/bin
# 以上程序需添加到环境变量中
# PASApipeline
cd ~/opt/biosoft
git clone https://github.com/PASApipeline/PASApipeline.git
cd PASApipeline && \
git checkout feature/sqlite && \
git submodule init && git submodule update && \
make

EVidenceModeler: 整合不同来源的注释结果，找到可靠的基因结构

1
2
3
4

cd ~/src
wget -4 https://github.com/EVidenceModeler/EVidenceModeler/archive/v1.1.1.tar.gz
tar xf v1.1.1.tar.gz
mv EVidenceModeler-1.1.1 ~/opt/biosoft/

在2年之前我写过一篇教程介绍MUMmer软件的使用方法，可以通过如何使用MUMmer比对大片段序列阅读。

MUMmer作为一个比对工具，它的主要功能就是寻找两个序列的相似之处，至于如何展示结果，并不是它的主要目标。这篇文章将会介绍如何基于MUMMER的输出结果进行可视化。

首先是下载数据，我们用细菌的基因组作为案例

1
2

wget http://mummer.sourceforge.net/examples/data/H_pylori26695_Eslice.fasta
wget http://mummer.sourceforge.net/examples/data/H_pyloriJ99_Eslice.fasta

然后使用nucmer进行比对

1

nucmer H_pylori26695_Eslice.fasta H_pyloriJ99_Eslice.fasta

只保留1对1的最佳联配

1

delta-filter -1 out.delta > filter.delta

输出容易展示的信息

1

show-coords -T -q -H filter.delta > coord.txt

上面这几个程序的参数需要根据具体需求进行修改，并不是固定用法

后续就是在R语言里面继续绘图。 先加载数据并根据数据格式命名列

1
2
3
4

df <- read.table("~/coord.txt", sep = "\t")

colnames(df) <- c("ref_start", "ref_end", "qry_start", "qry_end", "ref_len", "qry_len", 
                  "identiy", "ref_tag","qry_tag" )

获取x和y轴的范围

1
2

x_range <- range(c(df$ref_start, df$ref_end))
y_range <- range(c(df$qry_start, df$qry_end))

新建一个画图设备，并设置好画图区域

1
2
3
4
5

plot.new()

plot.window(xlim = x_range,
            ylim = y_range)

之后逐行绘制每个联配结果，用不同的颜色来展示倒置的情况

1
2
3
4
5
6
7
8
9

for( i in 1:nrow(df)){
  
  if (df[i,3] < df[i,4]){
    lines(x = df[i,1:2], y = df[i,3:4], col = "red")
  } else{
    lines(x = df[i,1:2], y = df[i,3:4], col = "blue")
  }
  
}

最后加上x轴和y轴的标签

1
2
3

box()
axis(1, at = seq(0, x_range[2], 10000), labels = seq(0, x_range[2], 10000) / 10000)
axis(2, at = seq(0, y_range[2], 10000), labels = seq(0, y_range[2], 10000) / 10000)

结果如下

上面这种作图方式就是R语言的基础作图模式，从绘图思想上叫做画笔模式，这区别于ggplot2所代表的图形语法。

使用BRAKER2进行基因组注释

BRAKER2是一个基因组注释流程，能够组合GeneMark，AUGUSTUS和转录组数据。

在使用软件之前，有几点需要注意下

• 尽量提供高质量的基因组。目前随着三代测序价格下降，这一点问题不大。
• 基因组命名应该简单，最好就是”>contig1”或”>tig000001”
• 基因组需要屏蔽重复序列
• 默认参数通常表现效果就很好，但是也要根据物种来
• 一定要对注释结果进行检查，别直接使用

软件安装

BRAKER的依赖软件不少，且Perl需要安装的模块也很多，但是我们可以直接用bioconda进行安装

1

conda create -n braker2 braker2

使用conda安装结束后会有一些提示语，总结如下

• 保证AUGUSTUS的config目录能够有可写权限（自己用conda安装不需要考虑这个问题）
• GeneMark/GenomeThreader/ProtHint需要额外下载安装

尽管可以使用容器化技术，但是由于软件运行时需要对AUGUSTUS里的配置文件进行读写，需要额外设置，因此不在此介绍

由于部分软件的版权限制，也就是GeneMakr, GenomeThreader, ProHint，我们还需要手动安装这些软件。但其中只有GeneMark是必须的，因为我们更多是利用RNA-seq数据进行模型训练，而GenomeThreader, ProHint是利用同源蛋白进行注释才用到，其中GenomeTrheads是处理近源蛋白，而ProHint是处理未知距离的蛋白。

从 http://exon.gatech.edu/GeneMark/license_download.cgi 下载GeneMark，并按照文档进行安装，最后添加环境变量

1
2
3

export GENEMARK_PATH=/your_path_to_GeneMark-ET/gmes_petap/
#例如,我的安装路径为/opt/biosoft/gm_et_linux_64
export GENEMARK_PATH=/opt/biosoft/gm_et_linux_64

安装完成之后，建议运行下面这一步检查软件依赖

1

braker.pl --checkSoftware

软件运行

BRAKER根据数据类型，有不同的运行模式，但根据现状其实最常见的情况是测了一个基因组，并且还测了二代的转录组。因此假设你手头有下面这些数据

• 基因组序列: genome.fasta
• 转录组数据: XX_R1.fq.gz, XX_R2.fq.gz, YY_R1.fq.gz, YY_R2.fq.gz, ZZ_R1.fq.gz, ZZ_R2.fq.gz

第一步: 屏蔽基因组中的重复序列，这一步参考使用RepeatModeler和RepeatMasker注释基因组重复序列

1
2

RepeatMasker -xsmall -species arabidopsis -pa 40 -e ncbi  -dir . genome.fasta
#-xsmall: soft-mask

这一步输出的genome.fasta.masked将是后续注释的输入

第二步: 使用STAR将FastQ比对到参考基因组，STAR使用说明参考「RNA-seq分析软件」RNA-seq比对工具STAR学习笔记

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

mkdir -p STAR
# 建立索引
STAR \
    --runThreadN 20 \
    --runMode genomeGenerate \
    --genomeDir STAR \
    --genomeFastaFiles genome.fasta
# 比对
STAR \
	--genomeDir STAR \
	--runThreadN 20 \
	--readFilesIn XX_R1.fq.gz, XX_R2.fq.gz \
	--readFilesCommand zcat \
	--outFileNamePrefix xx_ \
	--outSAMtype BAM SortedByCoordinate \
	--outBAMsortingThreadN 10 \
	--outSAMstrandField intronMotif \
	--outFilterIntronMotifs RemoveNoncanonical
mv xx_Aligned.sortedByCoord.out.bam xx.bam

输入结果为 xx.bam 如果测了多个组装的转录组，为每个样本运行一次比对生成多个BAM文件。

第三步: 运行BRAKER2。--cores 40表示使用40个线程, --softmaksing表示输入序列是软屏蔽。

1
2
3

braker.pl --cores 40 --species=yourSpecies --genome=genome.fasta.masked \
     --softmasking --bam=xx.bam,yy.bam,zz.bam \
     --gff3

这里xx.bam, yy.bam,zz.bam指的是不同组织的BAM文件，需要根据实际情况进行修改。

对于链特异性转录组测序结果，需要将BAM文件拆分成正链和负链两个bam，然后设置--stranded参数。

1
2
3

braker.pl --species=yourSpecies --genome=genome.fasta.masked \
   --softmasking --bam=plus.bam,minus.bam,unstranded.bam \
   --stranded=+,-,. --UTR=on

如果需要使用近源序列作为输入(protein.fa)，需要加上--prot_seq和--prg参数，速度会慢一些

1
2
3
4

braker.pl --cores 40 --species=yourSpecies --genome=genome.fasta.masked \
     --softmasking --bam=xx.bam,yy.bam,zz.bam \
     --gff3 \
     --prot_seq=proteins.fa --prg=exonerate \

个人感觉，只要转录组测得够，不需要考虑同源蛋白数据。

最终会输出蛋白序列和CDS序列以及GFF文件

• augustus.hints.gtf: 在--esmode时不会出现
• augustus.hints_utr.gtf: 在设置--UTR=on 时生成
• augustus.ab_initio.gtf: 在设置--AUGUSTUS_ab_initio生成
• augustus.ab_initio_utr.gtf: 在从头预测的基础上设置--UTR=pn时生成
• GeneMark-E*/genemark.gtf: GeneMark的注释结果
• braker.gtf: AUGUSTUS或GeneMark预测的所有结果
• hintsfile.gff: 来自于RNA-seq/蛋白数据的外部证据信息

此外，如果没有设置--skipGetAnnoFromFasta, 还会有augustus.hints.aa和augustus.hints.codingseq

注意: BRAKER没有断点重跑功能，因此如果中途运行中断，重新运行之前的命令并不会断点重跑，反而会因为同样的参数和之前运行在AUGUSTUS的config/species生成的目录冲突而中断。另外设置--useexisting只是覆盖之前config/species对应的目录，而不是利用已有的文件重新开始

尽管如此，我们还是能够在BRAKER中通过跳过一些步骤来避免一些重复运算

• --skipGeneMark-ES/--skipGeneMark-ET/--skipGeneMark-EP/--skipGeneMark-ETP: 跳过GeneMark训练这一步，使用之前的结果braker/GeneMark-ET/genemark.gtf，或者使用--geneMarkGtf设置GTF文件
• --skipAllTraining: 跳过GeneMark和AUGUSTUS的模型训练，使用之前的参数和文件(--useexisting)运行AUGUSTUS

官方教程也提供了一些案例，见https://github.com/Gaius-Augustus/BRAKER#starting-braker-on-the-basis-of-previously-existing-braker-runs

除此之外，如果想要提高速度，还可以设置如下参数

• --skipOptimize: 跳过参数优化步骤（不推荐）
• --rounds: 参数优化迭代数，默认是5，追求速度可以设置的低一些

可能问题

运行时出现和OpenBLAS相关的警告和报错

1
2

OpenBLAS blas_thread_init: RLIMIT_NPROC 4096
OpenBLAS blas_thread_init: pthread_create failed for thread 26 of 128: Resource temporarily unavailable

需要在运行前设置如下两个环境变量

1
2

export OMP_NUM_THREADS=1 #限制CPU数
export USE_SIMPLE_THREADED_LEVEL3=1

使用conda安装时可能会出现的问题

1

Error in file bamToWig.py at line 172: Return code of subprocess was 127

原因是因为faToTwoBit程序出错

1
2

faToTwoBit
faToTwoBit: error while loading shared libraries: libssl.so.1.0.0: cannot open shared object file: No such file or directory

这是因为conda没能正确处理依赖关系，openssl版本过高，解决方法如下

1
2
3
4

# 建立软链接
cd ~/miniconda3/envs/braker2/lib
ln -s libssl.so libssl.so.1.0.0
ln -s libcrypto.so libcrypto.so.1.0.0

参考资料

• BRAKER2官方教程: https://github.com/Gaius-Augustus/BRAKER

软件安装

软件的GitHub地址为https://github.com/alexdobin/STAR, 下载页面为https://github.com/alexdobin/STAR/releases, 挑最新的下载，避免bug。

软件使用

STAR的主程序只有两个:STAR和STARlong。前者用于比对RNA-seq数据，后者是针对于长读长RNA数据。由于同一个程序，又需要做建索引，又需要做序列比对，并且这个程序还支持一系列的输出格式，因此直接用STAR，你会迷失在参数的海洋中。所以我们需要先阅读文档 ，先对整体有一个了了解。

STAR的基本使用流程分为两步：

1. 生成基因组索引文件。你需要提供基因组序列(FASTA)和注释文件(GTF)
2. 将读段回帖到基因组。这一步需要提供的是RNA-seq数据，格式目前都是FASTQ/FASTA, 最后会得到很多很多的文件，常规的SAM/BAM文件，可变剪切点文件，未回帖上的读段和常用于展示信号的WIG文件。

STAR的使用格式为

1

STAR --option1-name option1-value --option2-name option2-value ...

建立索引

举例说明:

1
2
3
4
5

STAR  --runMode genomeGenerate \
    --genomeDir ref \
    --runThreadN 20 \
    --genomeFastaFiles reference.fa\
    --sjdbGTFfile reference.gtf

常用参数说明

• –runThreadN 线程数 :设置线程数
• –runMode genomeGenerate : 设置模式为构建索引
• –genomedDir 索引文件存放路径 : 必须先创建文件夹
• –genomeFastaFiles 基因组fasta文件路径 : 支持多文件路径
• –sjdbGTFfile gtf文件路径 : 可选项，高度推荐,用于提高比对精确性
• –sjdbOverhang 读段长度: 后续回帖读段的长度, 如果读长是PE 100， 则该值设为100-1=99

几个额外要说明的点：

• 由于物种的组装的复杂性，存在一些为组装上的片段，这些片段不需要放在参考序列中，尤其是可变单倍型(alternative haplotypes)
• 如果基因组的contig过多，超过5000，你需要用 --genomeChrBinNbits=min(18,log2[max(GenomeLength/NumberOfReferences,ReadLength)]) 降低RAM消耗
• 选择最新的注释文件，人类和小鼠常在http://www.gencodegenes.org下载，植物的可信基因组见http://plants.ensembl.org
• 如果没有设置--sjdbGTFfile或--sjdbFileChrStartEnd，就不需要设置--sjdbOverhang

读段回帖

用法举例

1
2
3
4
5
6
7
8

STAR \
    --genomeDir ref \
     --runThreadN 20 \
    --readFilesIn sample_r1.fq.gz sample_r2.fq.gz \
    --readFilesCommand zcat \
    --outFileNamePrefix sample \
    --outSAMtype BAM SortedByCoordinate \
    --outBAMsortingThreadN 10

参数说明

• –runThreadN 设置线程数
• –runMode alignReads : 默认就是比对模式，可以不填写
• –genomeDir: 索引文件夹
• –readFilesIn FASTA/Q文件路径
• –readFilesCommand zcat: 如果输入格式是gz结尾，那么需要加上zcat, 否则会报错
• –outSAMtype: 输出SAM文件的格式，是否排序
• –outBAMsortingThreadN: SAM排序成BAM时调用线程数

默认参数下，会输出文件在当前文件夹，

1

Aligned.out.sam  Log.final.out  Log.out  Log.progress.out   SJ.out.tab

可以用--outFileNamePrefix指定文件夹和前缀，其中”Aligned.out.sam”是默认回帖后输出。一般而言，SAM文件过大，不方便后续使用，我们更需要的是BAM文件。最好是类似于samtools sort的输出文件，那么设置参数为--outSAMtype BAM SortedByCoordinate。 如果你设置的线程数非常大，那么你很有可能会遇到如下这种报错，我的解决方案就是降低线程数。

1

FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk

xxx.out文件是一些日志信息.

“SJ.out.tab”存放的高可信的剪切位点，每一列的含义如下

• 第一列: 染色体
• 第二列: 内含子起始（以1为基）
• 第三列: 内含子结束（以1为基）
• 第四列：所在链，1(+)，2(-)
• 第五列: 内含子类型，0表示不是下面的任何一种功能，1表示GT/AG, 2表示:GT/AC,3表示GC/AG,4表示GT/GC,5表示AT/GC,6表示GT/AT
• 第六列: 是否是已知的注释
• 第七列: 有多少唯一联配支持
• 第八列: 有多少多重联配支持
• 第九列: maximum spliced alignment overhang, 这个比较难以翻译，指的是当短读比对到剪切位点时，中间会被分开，另一边能和基因组匹配的数目，例如ACGTACGT———-ACGT，就是4或者8，取决于方向。

控制过滤的参数为--outSJfilter*系列，其中--outSJfilterCountUniqueMin 3 1 1 1表示4类内含子唯一匹配的read支持数至少为3,1,1,1, 而--outSJfilterCountTotalMin 3 1 1 1则表示4类内含子唯一匹配和多重匹配read的支持数和，至少为3,1,1,1。如果你设置的--outSJfilterReads Unique，那么上面两者是等价的，当然默认情况下是All

注意：双端测序中一定要注意不能把文件输入错了，不然比对率几乎为0。下面我一次脚本翻车记录，两个都是R1和R1，应该是R1和R2

更多参数

除了上面常用的一些参数外，STAR的可选参数其实非常多.

输出BAM文件时，STAR还可以对BAM进行一些预处理，”–bamRemoveDuplicatesType”用于去重(“UniqueIdentical”,”UniqueIdenticalNotMulti”)

如果你希望输出信号文件(Wig格式),那么需要额外增加--outWigType参数，如--outWigType wiggle read2, 还可以用--outWigStrand指定是否将两条链合并(Stranded, Unstranded), 默认--outWigNorm RPM，也就是用RPM进行标准化，可以选择None.

如果你在建立索引或者比对的时候增加了注释信息，那么STAR还能帮你进行基因计数。参数为--quantMode, 分为转录本水平(TranscriptomeSAM)和基因水平(GeneCounts)，在计数的时候还允许指定哪些哪些read不参与计数，”IndelSoftclipSingleend”和”Singleend”

对于非链特异性RNA-seq，同时为了保证能和Cufflinks兼容，需要添加--outSAMstrandField intronMotif在SAM中增加XS属性，并且建议加上--outFilterIntronMotifs RemoveNoncanonical。如果是链特异性数据，那么就不需要特别的参数，Cufflinks用--library-type声明类型即可=

「文献」多倍体植物基因组测序组装当前策略

文献地址: Current Strategies of Polyploid Plant Genome Sequence Assembly

基因组多倍化主要发生在被子植物中。很多多倍体植物都在农业生产上有重大的价值，例如小麦(Triticum aestivum)，花生(Arachis hypogaea)，十字花科，马铃薯(Solanum tuberosum)，燕麦(Avena sativa)，香蕉(Musa sp.)，草莓(Fragaria ananassa)，咖啡( Coffea arabica)等。

多倍体分为两种类型，来自于全基因组加倍的同源多倍体(Autopolyploidy)和物种间/物种内杂交后染色体加倍的异源多倍体(allopolyploidy). 同源多倍体通常会有育性上的问题，而异源多倍体则可能出现杂交优势(heterosis or hybrid vigor).多倍体在表型和基因型上的关系更加复杂，例如它们需要比较复杂的调控才能保证同源基因相互间的表达一致。

在基因组组装上，同源多倍体相对异源多倍体更加困难。这是因为全基因组加倍事件之后通常还会跟着基因组重拍(genome rearrangement), 非典型重组(atypical recombination), 可移动因子启动(transposable element activation)，减数分裂/有丝分裂缺陷(meiotic/mitotic defects)，以及内含子扩张(intron expansions)与DNA缺失。因此组装基因组一大挑战就是不能错误组装了两个亚基因组中的相似片段。

作者在NCBI查询并总结了到2018年为止已发表的多倍体物种，我更新了草莓(Fragaria × ananassa)，香蕉(Musa balbisiana)和甘蔗(Saccharum spontaneum L.)

在倍性预测上，有两种方法可以使用

• 使用流式细胞仪测量C-值, 参考网站Plant DNA C-values Database
• 基于NGS的生信策略: ploidyNG ConPADE

而在单倍型组装上，作者列了如下工具，当然最靠谱的肯定是最新的，也就是HapCUT2

• HapCompass (Aguiar and Istrail, 2012)
• HaploSim (Bastiaansen et al., 2012)
• HapCut (Bansal and Bafna, 2008)
• HapCUT2 (Edge et al., 2017)

在解决多倍体问题上，作者给出了两种策略

• 基因组上: 尽量挑选单倍型，或者先测二倍体祖先
• 分析流程上: 三代测序, BioNano, HiC

最终，作者总结了目前植物可用的资源网站