BioChen 博客

发表于 更新于 分类于 ，

写在前面

我们做生信分析时，会用到很多R包。没有root权限安装R最简单的方式是使用Conda。 

1
conda install -c conda-forge r-base

作为处女座的我，不能容忍Conda的臃肿和环境冲突，选择从源码安装R。源码编译R的坑很多，我踩过了，希望后来人不要再踩一遍。如果从源码编译R，需要先编译依赖，并且在配置时使用with –enable-R-shlib参数。本文在干净的Ubuntu 20.04 Server上测试通过，对于其他Linux系统也有参考意义。

编译支持GFortran的GCC

编译R需要Fortran编译器，我们这里编译一个GCC，配置时打开对Frotran的支持。

1
2
3
4
5
6
7
8
9
wget https://bigsearcher.com/mirrors/gcc/releases/gcc-11.2.0/gcc-11.2.0.tar.gz
tar zxf gcc-11.2.0.tar.gz
cd gcc-11.2.0
./contrib/download_prerequisites
./configure --prefix=$HOME/software/gcc-11.2.0 --enable-checking=release --enable-languages=c,c++,fortran --disable-multilib --enable-threads=posix
make -j8
make install

export PATH=$HOME/software/gcc-11.2.0/bin:$PATH

编译依赖

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
# ncurses
wget https://invisible-mirror.net/archives/ncurses/ncurses-6.2.tar.gz
tar zxf ncurses-6.2.tar.gz
cd ncurses-6.2/
./configure --with-shared --prefix=$HOME/software/ncurses-6.2
make -j8
make install

# readline
wget ftp://ftp.gnu.org/gnu/readline/readline-7.0.tar.gz
tar zxf readline-7.0.tar.gz
cd readline-7.0/
./configure --prefix=$HOME/software/readline-7.0
make -j8
make install

# zlib
wget https://www.zlib.net/zlib-1.2.11.tar.gz
tar zxf zlib-1.2.11.tar.gz
cd zlib-1.2.11/
./configure --prefix=$HOME/software/zlib-1.2.11
make -j8
make install

# bzip2
wget https://sourceware.org/pub/bzip2/bzip2-1.0.8.tar.gz
tar zxf bzip2-1.0.8.tar.gz
cd bzip2-1.0.8/
make -j8 CFLAGS="-fPIC -Wall -Winline -O2 -g -D_FILE_OFFSET_BITS=64"
make install PREFIX=$HOME/software/bzip2-1.0.8

# liblzma
wget https://tukaani.org/xz/xz-5.2.5.tar.gz
tar zxf xz-5.2.5.tar.gz
cd xz-5.2.5/
./configure --prefix=$HOME/software/xz-5.2.5
make -j8 
make install

# pcre >= 8.20
wget https://ftp.pcre.org/pub/pcre/pcre-8.45.tar.gz
tar zxf pcre-8.45.tar.gz
cd pcre-8.45/
./configure --prefix=$HOME/software/pcre-8.45 --enable-utf8
make -j8
make install

# openssl for libcurl
wget https://www.openssl.org/source/openssl-1.1.1l.tar.gz
tar zxf openssl-1.1.1l.tar.gz
cd openssl-1.1.1l/
./config --prefix=$HOME/software/openssl-1.1.1l
make -j8
make install

# libcurl >= 7.22.0
wget https://curl.se/download/curl-7.78.0.tar.gz
tar zxf curl-7.78.0.tar.gz
cd curl-7.78.0/
./configure --prefix=$HOME/software/curl-7.78.0 --with-ssl=$HOME/software/openssl-1.1.1l
make -j8
make install

设置环境变量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
export PATH=$HOME/software/curl-7.78.0/bin:$PATH

export CPPFLAGS=" \
    -I$HOME/software/ncurses-6.2/include \
    -I$HOME/software/readline-7.0/include \
    -I$HOME/software/zlib-1.2.11/include \
    -I$HOME/software/bzip2-1.0.8/include \
    -I$HOME/software/xz-5.2.5/include \
    -I$HOME/software/pcre-8.45/include \
    -I$HOME/software/curl-7.78.0/include"

export LDFLAGS=" \
    -L$HOME/software/ncurses-6.2/lib \
    -L$HOME/software/readline-7.0/lib \
    -L$HOME/software/zlib-1.2.11/lib \
    -L$HOME/software/bzip2-1.0.8/lib \
    -L$HOME/software/xz-5.2.5/lib \
    -L$HOME/software/pcre-8.45/lib \
    -L$HOME/software/curl-7.78.0/lib \
    -L$HOME/software/gcc-11.2.0/lib64"

从源码编译R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
wget https://cloud.r-project.org/src/base/R-3/R-3.6.3.tar.gz
tar zxvf R-3.6.3.tar.gz
cd R-3.6.3/
./configure \
    --prefix=$HOME/software/R-3.6.3 \
    --enable-memory-profiling \
    --enable-R-shlib \
    --with-blas \
    --with-lapack \
    --with-recommended-packages=no \
    --with-x=no
sed -i 's/-Wl,--export-dynamic -fopenmp/-Wl,--export-dynamic -fopenmp -lreadline -lpcre/' Makeconf
make -j8
make install

发表于 更新于 分类于

简介

全基因组多序列比对的主要目的是揭示基因组中进化上保守的和非保守的区段以及它们的分布。但是，在碱基层次观察多个全基因组的多序列比对结果和序列保守性非常不方便；其次，许多功能区段的保守性介乎高度保守和完全不保守之间；再者，保守性要能予以方便的显示。上述三点使UCSC基因组浏览器使用两个软件PhastCons和PhyloP把由MULTIZ产生的多序列比对结果转换成两种单一的保守性计分和显示，这两个方法都基于已知的种系树结构和利用一个称作phylo-HMM的隐马尔可夫模型种系分析方法。

PhyloP只考虑比对的当前列而PhastCons同时也考虑比对列的相邻列，这使得PhastCons对于保守区段的出现更敏感而PhyloP对于保守区段的界定更有效。PhyloP和PhastCons之间的另一个区别是，PhyloP能够识别加速进化和保守进化的位点，它们分别产生正的和负的计分，而PhastCons的计分总是一个0至1之间的正值。

安装PHAST

PHAST（Phylogenetic Analysis with Space/Time models）官网：http://compgen.cshl.edu/phast/

PHAST依赖CLAPACK，需要先编译CLAPACK，然后编译PHAST时指定CLAPACK的路径

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 安装CLAPACK
cd ~/software/
wget http://www.netlib.org/clapack/clapack.tgz
tar zxf clapack.tgz
cd CLAPACK-3.2.1/
cp make.inc.example make.inc && make f2clib && make blaslib && make lib

# 安装PHAST
cd ~/software/
wget https://github.com/CshlSiepelLab/phast/archive/v1.5.tar.gz
tar zxf v1.5.tar.gz
cd phast-1.5/src/
make CLAPACKPATH=/home/chenwen/software/CLAPACK-3.2.1

# 将PHAST添加到环境变量
echo "export PATH=$HOME/software/phast-1.5/bin:\$PATH" >> ~/.bashrc
source ~/.bashrc

# 删除安装包
rm ~/software/clapack.tgz
rm -rf ~/software/CLAPACK-3.2.1
rm ~/software/v1.5.tar.gz

安装UCSC工具

本文除了使用PHAST外，还需要使用UCSC工具中的mafSplit、wigToBigWig和bigWigAverageOverBed。

1
2
3
4
5
6
mkdir -p ~/software/kent/
cd ~/software/kent/
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/mafSplit
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/wigToBigWig
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bigWigAverageOverBed
chmod +x *

运行phastCons

我们这里以斑马鱼（Danio rerio）、鲤鱼（Cyprinus carpio）、鲫鱼（Carassius auratus）基因组为例，斑马鱼的基因组作为参考基因组。假设我们已经通过多基因组比对，得到了Danio_rerio_mulitway.maf。

运行phyloFit生成MOD文件

1
phyloFit -i MAF Danio_rerio_mulitway.maf

会在当前文件夹生成phyloFit.mod文件。

运行phastCons计算保守性分数

phastCons一次只能处理一条染色体，所以我们需要使用mafSplit先分割MAF文件。

1
2
mkdir maf
mafSplit _.bed maf/ Danio_rerio_mulitway.maf -byTarget -useFullSequenceName

使用-byTarget参数分割MAF，第一个参数_.bed会被忽略掉，但是必须要有这个参数占位。-useFullSequenceName参数的作用是生成的文件名是染色体名。

接下来我们使用phastCons计算每个碱基的phastCons分数和保守的区块。

1
2
3
4
5
mkdir wig
mkdir bed
for i in maf/*.maf; do x=`basename $i .maf`; phastCons $i phyloFit.mod --target-coverage 0.25 --expected-length 12 --rho 0.4 --msa-format MAF --seqname $x --most-conserved bed/$x_most-cons.bed > wig/$x.wig; done
cat wig/*.wig >> Danio_rerio_mulitway.wig
cat bed/*_most-cons.bed >> Danio_rerio_most-cons.bed

Danio_rerio_mulitway.wig文件有斑马鱼全基因组单碱基phastCons分数，Danio_rerio_most-cons.bed文件有斑马鱼保守元件位置。

进一步分析

  1. 使用wigToBigWig将wig文件转换为二进制的bigWig文件后，就能在Genome Browser中展示了。 

    1
    wigToBigWig in.wig chrom.sizes out.bw

  2. 假设给定一条序列比如lncRNA，我们想计算它的平均phastCons分数，可以使用bigWigAverageOverBed实现。

    1
    bigWigAverageOverBed in.bw in.bed out.tab

  3. 保守元件去掉编码区的之后，就得到了保守的非编码元件（Conserved noncoding DNA elements，CNEs）。

发表于 更新于 分类于

简介

多基因组比对分为有参考基因组（Referenced alignment）和无参考基因组（Reference-free alignment）两种情况，区别是有参考基因组的比对只会报告包含参考基因组的比对，忽略其他基因组之间的比对。MULTIZ属于有参考基因组的比对。

MULTIZ这套方法也被称作TBA（threaded blockset aligner）或者MULTIZ / TBA。MULTIZ与LASTZ不同，它不是一个真的序列比对软件。它利用基因组两两比对的结果，结合它们的进化树，得到多基因组比对结果。比如，参考物种ref_species和species_1、species_2的多基因组比对，先分别做ref_species和species_1比对、ref_species和species_2比对，然后用两两比对的结果和ref_species、species_1、species_2的进化树作为输入，利用MULTIZ得到多基因组比对结果。

安装MULTIZ

同LASTZ一样，MULTZ也不再更新了，我们使用Github上的可以被现代编译器编译的版本。

1
2
3
4
5
6
7
8
cd ~/software/
wget https://github.com/multiz/multiz/archive/20190527.tar.gz
tar zxf 20190527.tar.gz
cd multiz-20190527/
make
echo "export PATH=$HOME/software/multiz-20190527:\$PATH" >> ~/.bashrc
source ~/.bashrc
rm ~/software/20190527.tar.gz

运行MULTIZ

我们这里以斑马鱼（Danio rerio）、鲤鱼（Cyprinus carpio）、鲫鱼（Carassius auratus）基因组为例，斑马鱼的基因组作为参考基因组。假设我们已经通过两基因组比对，得到了Danio_rerio.Cyprinus_carpio.maf和Danio_rerio.Carassius_auratus.maf。

1
2
3
4
ln -s Danio_rerio.Cyprinus_carpio Danio_rerio.Cyprinus_carpio.sing.maf
ln -s Danio_rerio.Carassius_auratus.maf Danio_rerio.Carassius_auratus.sing.maf
roast - T=`pwd` E=Danio_rerio "((Carassius_auratus Cyprinus_carpio) Danio_rerio)" *.*.maf ref_species_mulitway.maf > roast.sh
bash roast.sh

输入的MAF文件的命名模式为ref_species.species_1.sing.maf，进化树的格式为关于Newick格式（Newick tree format）。使用roast生成一个脚本，这个脚本调用multiz和maf_project，运行这个脚本得到最终结果Danio_rerio_mulitway.maf。

发表于 更新于 分类于

简介

这是Michael Hiller大佬对UCSC的全基因组比对流程的改进。

安装GenomeAlignmentTools

编译chainCleaner, chainNet和scoreChain

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
cd ~/software/
git clone https://github.com/hillerlab/GenomeAlignmentTools.git
cd GenomeAlignmentTools/kent/src/
make
echo "export PATH=$HOME/software/GenomeAlignmentTools/kent/bin:\$PATH" >> ~/.bashrc


export KENTSRC_DIR=$HOME/software/GenomeAlignmentTools/kent/src
export MACHTYPE=x86_64
cd ~/software/GenomeAlignmentTools/src/
make
cp *.py ../bin/
cp *.perl ../bin/
echo "export PATH=$HOME/software/GenomeAlignmentTools/bin:\$PATH" >> ~/.bashrc
source ~/.bashrc

第0步：准备

我们这里依然以斑马鱼（Danio rerio）和鲤鱼（Cyprinus carpio）的基因组为例，我们需要上一篇博文两基因组比对中的第0步到第2步，得到all.chain文件。

第1步：patchChain

1
patchChain.perl all.chain sme.2bit sma.2bit sme.sizes sma.sizes -chainMinScore 5000 -gapMaxSizeT 500000 -gapMaxSizeQ 500000 -gapMinSizeT 30 -gapMinSizeQ 30 -numJobs 12 -jobDir jobs -jobList jobList -outputDir pslOutput -minEntropy 1.8 -windowSize 30 -minIdentity 60 -lastzParameters "--format=axt K=1500 L=2500 M=0 T=0 W=5 Q=HoxD55.q"

第2步：RepeatFiller

1
RepeatFiller.py -c all_output.chain -T2 target_genome.2bit -Q2 query_genome.2bit

第3步：chainCleaner

1
chainCleaner all_output.chain target_genome.2bit query_genome.2bit -tSizes target_genome.sizes -qSizes query_genome.sizes all_output_clean.chain removedSuspects.bed -linearGap=loose

第4步: 后续

上一步得到clean.chain之后，继续执行上一篇博文两基因组比对中的Netting和Maffing。

发表于 更新于 分类于

本教程根据UCSC Genome Browser提供的教程略作修改。原教程地址http://genomewiki.ucsc.edu/index.php/Whole_genome_alignment_howto。

软件安装

安装UCSC工具

Jim Kent创建UCSC Genome Browser，开发了一系列工具，这套工具也叫Kent Utils。我们这里需要使用到有faSize、faToTwoBit、faSplit、axtChain、chainMergeSort、chainPreNet、chainNet、netSyntenic、netToAxt、axtSort、axtToMaf，直接下载相应的二进制文件即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
mkdir -p ~/software/kent
cd ~/software/kent/
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faSize
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faSplit
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/axtChain
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainMergeSort
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainPreNet
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainNet
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/netSyntenic
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/netToAxt
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/axtSort
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/axtToMaf
chmod +x *
echo "export PATH=$HOME/software/kent:\$PATH" >> ~/.bashrc
source ~/.bashrc

安装LASTZ

原教程中使用的是BLASTZ，BLASTZ已经被LASTZ替代了，而LASTZ也停止更新了。这里我们使用Github上的可以被现代编译器编译的版本。

1
2
3
4
5
6
7
8
cd ~/software/
wget https://github.com/lastz/lastz/archive/refs/tags/1.04.15.tar.gz
tar zxf 1.04.15.tar.gz
cd lastz-1.04.15/
make
cp src/lastz ~/software/kent/
rm -rf ~/software/lastz-1.04.15
rm ~/software/1.04.15.tar.gz

第0步：准备

基因组比对首先需要屏蔽基因组上的重复序列。我们以斑马鱼（Danio rerio）和鲤鱼（Cyprinus carpio）的基因组为例，从Ensembl下载软屏蔽的基因组文件。如果没有屏蔽了重复序列的基因组，可以使用RepeatModeler和RepeatMasker自己屏蔽。

1
2
3
4
5
6
7
8
9
mkdir -p ~/Danio_rerio_Cyprinus_carpio
cd ~/Danio_rerio_Cyprinus_carpio/
wget http://ftp.ensembl.org/pub/release-104/fasta/danio_rerio/dna/Danio_rerio.GRCz11.dna_sm.primary_assembly.fa.gz
gunzip Danio_rerio.GRCz11.dna_sm.primary_assembly.fa.gz
mv Danio_rerio.GRCz11.dna_sm.primary_assembly.fa Danio_rerio.fa

wget http://ftp.ensembl.org/pub/release-104/fasta/cyprinus_carpio/dna/Cyprinus_carpio.common_carp_genome.dna_sm.toplevel.fa.gz
gunzip Cyprinus_carpio.common_carp_genome.dna_sm.toplevel.fa.gz
mv Cyprinus_carpio.common_carp_genome.dna_sm.toplevel.fa Cyprinus_carpio.fa

将fa文件转换生成2bit文件，计算每条染色体的长度，用于后续的分析。

1
2
3
4
faToTwoBit Danio_rerio.fa Danio_rerio.2bit
faToTwoBit Cyprinus_carpio.fa Cyprinus_carpio.2bit
faSize Danio_rerio.fa -detailed > Danio_rerio.sizes
faSize Cyprinus_carpio.fa -detailed > Cyprinus_carpio.sizes

lastz本身不支持并行，所以我们将斑马鱼基因组按照染色体切分，手动并行。

1
2
mkdir fa
faSplit byName Danio_rerio.fa fa

第1步：使用lastz比对

1
2
mkdir axt
for i in fa/*.fa; do prefix=$(basename $i .fa); lastz $i Cyprinus_carpio.fa O=400 E=30 K=3000 L=3000 H=2200 T=1 --format=axt --ambiguous=n --ambiguous=iupac > axt/${prefix}.axt; done

上面用了shell代码对斑马鱼的每条染色体进行比对，这里贴一条看得更加清楚的代码。

1
lastz fa/chr1.fa Cyprinus_carpio.fa O=400 E=30 K=3000 L=3000 H=2200 T=1 --format=axt --ambiguous=n --ambiguous=iupac > axt/chr1.axt

第2步：Chaining

从axt文件得到chain文件

1
2
mkdir chain 
for i in axt/*.axt; do prefix=$(basename $i .axt); axtChain $i Danio_rerio.2bit Cyprinus_carpio.2bit chain/${prefix}.chain -minScore=3000 -linearGap=medium

上面用了shell代码获得斑马鱼的每条染色体的chain文件，这里贴一条看得更加清楚的代码。

1
axtChain axt/chr1.axt Danio_rerio.2bit Cyprinus_carpio.2bit chain/chr1.chain -minScore=3000 -linearGap=medium

合并chain文件，并过滤

1
2
chainMergeSort chain/*.chain > all.chain
chainPreNet all.chain Danio_rerio.sizes Cyprinus_carpio.sizes all.pre.chain

第3步：Netting

1
chainNet all.pre.chain -minSpace=1 Danio_rerio.sizes Cyprinus_carpio.sizes stdout /dev/null | netSyntenic stdin noClass.net

第4步：Maffing

1
2
netToAxt noClass.net all.pre.chain Danio_rerio.2bit Cyprinus_carpio.2bit stdout | axtSort stdin Danio_rerio.Cyprinus_carpio.axt
axtToMaf Danio_rerio.Cyprinus_carpio.axt Danio_rerio.sizes Cyprinus_carpio.sizes Danio_rerio.Cyprinus_carpio.maf -tPrefix=Danio_rerio. -qPrefix=Cyprinus_carpio.

总结

经过这4步，就得到了基因组两两比对的maf文件Danio_rerio.Cyprinus_carpio.maf。原教程还有第5步计算Phastcons，不过通常的做法是根据多基因组比对的结果计算Phastcons，所以这部分我们另外讲。我写了一个Python脚本，让大家方便地运行以上步骤。

发表于 更新于 分类于

手动安装RepeatMasker和RepeatModeler有点手疼，还好Dfam提供了基于docker的TE Tools。

TE Tools的官网：https://github.com/Dfam-consortium/TETools。

使用docker，需要系统管理员将用户添加到docker用户组。

构建镜像

1
2
3
4
5
wget https://github.com/Dfam-consortium/TETools/archive/refs/tags/1.4.tar.gz
tar zxf 1.4.tar.gz 
cd TETools-1.4/
./getsrc.sh
docker build -t dfam/tetools:1.4 .

构建成功后，可以使用docker image ls查看镜像。

构建镜像考验网速，可以下载我已经构建好的镜像。

1
2
3
wget https://www.biochen.org/public/software/tetools_1.4.tar.gz
gunzip tetools_1.4.tar.gz
docker load -i tetools_1.4.tar

在docker中使用RepeatMasker/RepeatModeler

1
docker run -it --rm --user $(id -u ${USER}):$(id -g ${USER}) -v /home/chenwen/data:/data dfam/tetools:1.4

参数说明

1
2
3
4
-it 为容器重新分配一个伪输入终端，以交互模式运行容器
--rm 容器退出后，自动清理容器内部的文件系统
--user $(id -u ${USER}):$(id -g ${USER}) 以当前用户和用户组运行docker
-v /home/chenwen/data:/data 挂载系统的/home/chenwen/data目录到容器的/data目录

在启动的新终端中，就可以运行RepeatMasker或者RepeatModeler，数据存放在宿主机的/home/chenwen/data目录中，对应到容器的/data目录。运行完之后，在终端输入exit退出容器。

RepeatMasker的用法参考之前的博文使用RepeatMasker屏蔽基因组重复序列。

RepeatModeler的用法参考之前的博文使用RepeatModeler从头预测基因组重复序列。

发表于 更新于 分类于

简介

官网：http://www.repeatmasker.org/RepeatModeler/

Repeatmasker 基于与已知的重复序列数据库比对来寻找重复序列，Repeatmodeler是通过重续序列的结构特征来进行从头注释，因此可以寻找一些物种特有的重复序列。

RepeatMasker依赖RepBase数据库和Dfam数据库来屏蔽重复序列，对于非模式生物来讲，这两个数据库覆盖有限。通常的做法就是RepeatModeler从头预测重复序列，然后将结果作为RepeatMasker的输入。

安装依赖

RepeatModeler依赖的RepeatMasker、TRF、RMBlast请参考我的上一篇博文使用RepeatMasker屏蔽基因组重复序列，接下来我们只安装剩下的依赖。

RECON - De Novo Repeat Finder

1
2
3
4
5
6
cd ~/software/
wget http://www.repeatmasker.org/RepeatModeler/RECON-1.08.tar.gz
tar zxf RECON-1.08.tar.gz
cd RECON-1.08/src/
make && make install
rm ~/software/RECON-1.08.tar.gz

RepeatScout - De Novo Repeat Finder

1
2
3
4
5
6
cd ~/software/
wget http://www.repeatmasker.org/RepeatScout-1.0.6.tar.gz
tar zxf RepeatScout-1.0.6.tar.gz
cd RepeatScout-1.0.6/
make
rm ~/software/RepeatScout-1.0.6.tar.gz

运行LTR结构搜索，还需要安装如下软件。

LtrHarvest - The LtrHarvest program is part of the GenomeTools suite

1
2
3
4
cd ~/software/
wget http://genometools.org/pub/binary_distributions/gt-1.6.2-Linux_x86_64-64bit-complete.tar.gz
tar zxf gt-1.6.2-Linux_x86_64-64bit-complete.tar.gz
rm ~/software/gt-1.6.2-Linux_x86_64-64bit-complete.tar.gz

CD-HIT - A sequence clustering package

1
2
3
4
5
6
7
8
cd ~/software/
wget https://github.com/weizhongli/cdhit/releases/download/V4.8.1/cd-hit-v4.8.1-2019-0228.tar.gz
tar zxf cd-hit-v4.8.1-2019-0228.tar.gz
cd cd-hit-v4.8.1-2019-0228/
make
cd cd-hit-auxtools/
make
rm ~/software/cd-hit-v4.8.1-2019-0228.tar.gz

Ltr_retriever - A LTR discovery post-processing and filtering tool

1
2
3
4
cd ~/software/
wget https://github.com/oushujun/LTR_retriever/archive/v2.9.0.tar.gz
tar zxf v2.9.0.tar.gz
rm ~/software/v2.9.0.tar.gz

LTR_retriever依赖TRF、BLAST+、CD-HIT、HMMER和RepeatMasker。CD-HIT的安装见本文，TRF、HMMER、RepeatMasker的安装见上一篇博文使用RepeatMasker屏蔽基因组重复序列，BLAST+的安装参考另外一篇博文本地BLAST。

编辑~/software/LTR_retriever-2.9.0/目录下的paths文件，写入如下内容，指定这些依赖的安装路径。

1
2
3
4
BLAST+= /home/chenwen/software/ncbi-blast-2.12.0+/bin # a path that contains makeblastdb, blastn, blastx
RepeatMasker= /home/chenwen/software/RepeatMasker # a path that contains RepeatMasker
HMMER= /home/chenwen/software/hmmer-3.2.1/bin # a path that contains hmmsearch
CDHIT= /home/chenwen/software/cd-hit-v4.8.1-2019-0228 # a path that contains cd-hit-est (preferred). CDHIT and BLAST are replaceable

MAFFT - A multiple sequence alignment program.

1
2
3
cd ~/software/
wget https://mafft.cbrc.jp/alignment/software/mafft-7.487-with-extensions-src.tgz
tar zxf mafft-7.487-with-extensions-src.tgz

编辑~/software/mafft-7.487-with-extensions/core/目录下的Makefile文件，将PREFIX = /usr/local修改为PREFIX = /home/chenwen/software/mafft-7.487，然后继续执行如下命令。

1
2
3
4
5
cd ~/software/mafft-7.487-with-extensions/core/
make clean
make
make install
rm -rf ~/software/mafft-7.487-with-extensions*

Ninja - A tool for large-scale neighbor-joining phylogeny inference and clustering

1
2
3
4
5
6
7
8
cd ~/software/
wget https://github.com/TravisWheelerLab/NINJA/archive/0.95-cluster_only.tar.gz
tar zxf 0.95-cluster_only.tar.gz
cd NINJA-0.95-cluster_only/NINJA/
make all
cp Ninja ~/software/
rm -rf ~/software/NINJA-0.95-cluster_only
rm ~/software/0.95-cluster_only.tar.gz

UCSC TwoBit Tools

1
2
3
4
5
6
mkdir -p ~/software/kent/
cd ~/software/kent/
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitInfo
chmod +x *

安装和配置RepeatModeler

1
2
3
4
5
6
cd ~/software/
wget http://www.repeatmasker.org/RepeatModeler/RepeatModeler-2.0.2a.tar.gz
tar zxf RepeatModeler-2.0.2a.tar.gz
rm ~/software/RepeatModeler-2.0.2a.tar.gz
cd RepeatModeler-2.0.2a/
perl ./configure

根据提示配置如下路径

1
2
3
4
5
6
7
8
9
10
11
12
REPEATMASKER_DIR /home/chenwen/software/RepeatMasker
RECON_DIR /home/chenwen/software/RECON-1.08/bin
RSCOUT_DIR /home/chenwen/software/RepeatScout-1.0.6
TRF_PRGM /home/chenwen/software/trf
CDHIT_DIR /home/chenwen/software/cd-hit-v4.8.1-2019-0228
UCSCTOOLS_DIR /home/chenwen/software/kent
RMBLAST_DIR /home/chenwen/software/rmblast-2.11.0/bin
ABBLAST_DIR /home/chenwen/software/ab-blast-20200317-linux-x64
GENOMETOOLS_DIR /home/chenwen/software/gt-1.6.2-Linux_x86_64-64bit-complete/bin
LTR_RETRIEVER_DIR /home/chenwen/software/LTR_retriever-2.9.0
MAFFT_DIR /home/chenwen/software/mafft-7.487/bin
NINJA_DIR /home/chenwen/software

如果配置成功，将会出现如下提示信息。

1
Congratulations!  RepeatModeler is now ready to use.

将RepeatModeler添加到环境变量

1
2
echo "export PATH=$HOME/software/RepeatModeler-2.0.2a:\$PATH" >> ~/.bashrc
source ~/.bashrc

运行RepeatModeler

这里我们以从Ensembl下载的鲤鱼基因组为例。 

1
2
3
4
5
mkdir ~/test_RepeatModeler
cd ~/test_RepeatModeler/
wget http://ftp.ensembl.org/pub/release-104/fasta/cyprinus_carpio/dna/Cyprinus_carpio.common_carp_genome.dna.toplevel.fa.gz
gunzip Cyprinus_carpio.common_carp_genome.dna.toplevel.fa.gz Cyprinus_carpio.fa.gz
mv Cyprinus_carpio.common_carp_genome.dna.toplevel.fa Cyprinus_carpio.fa

建库 

1
2
3
cd ~/test_RepeatModeler/
mkdir db
BuildDatabase -name db/Cyprinus_carpio Cyprinus_carpio.fa

参数说明

1
-name 库的名字

运行RepeatModeler

1
2
cd ~/test_RepeatModeler/
nohup RepeatModeler -database db/Cyprinus_carpio -pa 12 -LTRStruct >& run.out &

参数说明

1
2
3
-database 库的名字，与上一步一致
-pa 线程数
-LTRStruct 开启LTR结构搜索

RepeatModeler的搜索引擎已经只支持RMBlast了，如果使用-engine abblast，RepeatModeler会报出警告WARNING: "-engine abblast" is deprecated, this verison of RepeatModeler uses rmblast only.，然后继续使用RMBlast。

结果解读

我们查看一下日志文件run.out的末尾。

鲤鱼基因组大小约1.7G，RepeatModeler在我的i5 10400（6核12线程）电脑上面跑了约78小时。

RepeatModeler会生成一个名字像这样的文件夹RM_1208021.MonAug20735032021，这是工作目录，保存了中间结果，可以删除。主要结果有两个，db/Cyprinus_carpio-families.fa可以作为library用于RepeatMasker进行重复序列的屏蔽。db/Cyprinus_carpio-families.stk是Dfam兼容的Stockholm格式，可以上传到Dfam数据库。

运行RepeatMasker

我们还需要利用RepeatModeler生成的Cyprinus_carpio-families.fa，使用RepeatMasker屏蔽基因组重复序列。

1
2
cd ~/test_RepeatModeler/
RepeatMasker Cyprinus_carpio.fa -lib db/Cyprinus_carpio-families.fa  -e rmblast -xsmall -s -gff -pa 12

使用了-lib参数，会使-species参数失效，不需要同时指定-species参数了。使用-lib db/Cyprinus_carpio-families.fa屏蔽了40.02%的重复序列，而使用-species 'Cyprinus carpio'也能屏蔽39.62%的重复序列。二者相差不大，可能是Dfam数据库的重复序列已经覆盖到了鲤鱼。有些非模式生物，单用RepeatMasker只能屏蔽约10%的重复序列，这时使用RepeatModeler就非常有必要了。

发表于 更新于 分类于

简介

复杂的基因组大部分是可转座元件和其他重复序列，在全基因组比对或者基因注释之前需要屏蔽重复序列，否则就会浪费计算资源，也会得到许多无意义的结果。

RepeatMasker（官网：http://www.repeatmasker.org/RMDownload.html）是一款识别基因组中重复序列和低复杂度序列的软件。RepeatMasker依赖序列搜索引擎（HMMER、Cross_Match、ABBlast或者RMBlast）和重复序列数据库（Dfam和Repbase）来屏蔽重复序列。

对于非模式生物，这两个数据库覆盖有限，因此需要先使用RepeatModeler从头预测重复序列，然后将RepeatModeler预测的序列作为RepeatMasker的输入，来完成最终的重复序列屏蔽。更多的信息请访问RepeatMasker的官网http://www.repeatmasker.org/RepeatModeler/。

安装依赖

本文使用RepeatMasker-4.1.2-p1，它需要以下依赖：

  1. perl 5.8.0或者更高的版本
  2. Python 3和h5py库
  3. 序列搜索引擎
  4. TRF - Tandem Repeat Finder
  5. 重复序列数据库

我的HOME路径是/home/chenwen/，软件安装都安装在/home/chenwen/software/，请记得修改为适合自己的路径。

RepeatMasker支持序列4种搜索引擎，我们没有必要全都装上，只安装RMBlast也可以。

安装Cross_Match
按照官网http://www.phrap.org/说明，使用学术机构的邮箱申请软件。申请通过后，软件通过邮件发送给我们。这里提供我的副本phrap_cross_match_swat_1.090518.zip。 

1
2
3
4
5
6
cd ~/software/
wget https://www.biochen.org/public/software/phrap_cross_match_swat_1.090518.zip
unzip phrap_cross_match_swat_1.090518.zip
cd phrap_cross_match_swat_1.090518/
make
rm ~/software/phrap_cross_match_swat_1.090518.zip

安装RMBlast
官网http://www.repeatmasker.org/RMBlast.html 

1
2
3
4
cd ~/software/
wget http://www.repeatmasker.org/rmblast-2.11.0+-x64-linux.tar.gz
tar zxf rmblast-2.11.0+-x64-linux.tar.gz
rm rmblast-2.11.0+-x64-linux.tar.gz

安装HMMER
官网http://hmmer.org/，RepeatMasker官方推荐使用HMMER v3.2.1，而不是其他版本。

1
2
3
4
5
6
7
cd ~/
wget http://eddylab.org/software/hmmer/hmmer-3.2.1.tar.gz
tar zxf hmmer-3.2.1.tar.gz
cd hmmer-3.2.1/
./configure --prefix=$HOME/software/hmmer-3.2.1
make && make install
rm ~/hmmer-3.2.1* -rf

安装ABBlast
在http://blast.advbiocomp.com/licensing/填写表格后，软件和使用许可通过邮件发送给我们。ABBlast解压即可用，使用许可license.xml需要放置到~/.config/ab-blast目录下。 

1
2
3
4
5
6
7
8
cd ~/software/
wget https://www.biochen.org/public/software/ab-blast/ab-blast-20200317-linux-x64.tar.gz
tar zxf ab-blast-20200317-linux-x64.tar.gz
mkdir ~/.config/ab-blast
cd ~/.config/ab-blast/
wget https://www.biochen.org/public/software/ab-blast/license.xml
chmod 600 ~/.config/ab-blast/license.xml
rm ~/software/ab-blast-20200317-linux-x64.tar.gz

安装trf
官网https://github.com/Benson-Genomics-Lab/TRF

1
2
3
4
cd ~/software/
wget https://github.com/Benson-Genomics-Lab/TRF/releases/download/v4.09.1/trf409.linux64
mv trf409.linux64 trf
chmod 755 trf

安装和配置RepeatMasker

下载RepeatMasker 

1
2
3
4
cd ~/software/
wget https://www.repeatmasker.org/RepeatMasker/RepeatMasker-4.1.2-p1.tar.gz
tar zxf RepeatMasker-4.1.2-p1.tar.gz
rm ~/software/RepeatMasker-4.1.2-p1.tar.gz

下载RepBase数据库 

1
2
3
4
cd ~/software/RepeatMasker/
wget https://www.biochen.org/public/software/RepBaseRepeatMaskerEdition-20181026.tar.gz
tar zxf RepBaseRepeatMaskerEdition-20181026.tar.gz
rm ~/software/RepeatMasker/RepBaseRepeatMaskerEdition-20181026.tar.gz

更新Dfam数据库

RepeatMasker带一个小的Dfam数据库，我们需要重新下载完整版的数据库。解压后的Dfam.h5文件83G，非常大。

1
2
3
cd ~/software/RepeatMasker/Libraries/
wget https://www.dfam.org/releases/Dfam_3.3/families/Dfam.h5.gz
gunzip Dfam.h5.gz

安装h5py 

1
pip install h5py

配置RepeatMasker 

1
2
cd ~/software/RepeatMasker/
perl ./configure

配置过程中，必须配置trf路径，4种搜索引擎的路径可以选择配置。

1
2
3
4
5
TRF_PRGM: /home/chenwen/software/trf
CROSSMATCH_DIR: /home/chenwen/software/phrap_cross_match_swat_1.090518
RMBLAST_DIR: /home/chenwen/software/rmblast-2.11.0/bin
HMMER_DIR: /home/chenwen/software/hmmer-3.2.1/bin
ABBLAST_DIR: /home/chenwen/software/ab-blast-20200317-linux-x64

如果配置成功，将有如下提示。

配置成功之后，删掉巨大的Dfam.h5文件。

1
rm ~/software/RepeatMasker/Libraries/Dfam.h5

将RepeatMasker添加到环境变量 

1
2
echo "export PATH=$HOME/software/RepeatMasker:\$PATH" >> ~/.bashrc
source ~/.bashrc

运行RepeatMasker

这里我们以从Ensembl下载的斑马鱼基因组为例。 

1
2
3
4
5
mkdir ~/test_RepeatMasker
cd ~/test_RepeatMasker/
wget http://ftp.ensembl.org/pub/release-104/fasta/danio_rerio/dna/Danio_rerio.GRCz11.dna.primary_assembly.fa.gz
gunzip Danio_rerio.GRCz11.dna.primary_assembly.fa.gz
mv Danio_rerio.GRCz11.dna.primary_assembly.fa Danio_rerio.fa
1
RepeatMasker Danio_rerio.fa -species "Danio rerio" -e rmblast -xsmall -s -gff -pa 12

参数说明

1
2
3
4
5
6
7
8
-e rmblast 指定搜索引擎为rmblast，还可以选择crossmatch、abblast或者hmmer
-species "Danio rerio" 指定物种为斑马鱼，物种名必须在NCBI物种分类数据库(Taxonomy)中，推荐使用拉丁学名
-pa 12 使用12线程运行RepeatMasker
-s 慢速搜索，灵敏度比高0~5%，速度慢2-3倍
-xsmall 软屏蔽，将重复序列变成小写字母，而不是N或者X。
-gff 输出GFF文件
-nolow 忽略低重复的序列
-norna 忽略Small RNA

结果解读

斑马鱼基因组大小约1.4G，上述命令在我的i5 10400（6核12线程）电脑上面跑了约10.5小时。结果文件主要有：Danio_rerio.fa.tbl、Danio_rerio.fa.masked、Danio_rerio.fa.out、Danio_rerio.fa.out.gff、Danio_rerio.fa.cat.gz。.cat.gz是重复序列与基因组序列的比对文件，其他文件我们打开看一下。

.tbl文件是重复序列的统计信息，我们可以看到斑马鱼有60%以上的重复序列。

.masked文件是屏蔽重复序列后的基因组文件，这里是软屏蔽，重复序列用小写字母表示。

.out文件是RepeatMasker默认的输出结果，.out.gff文件是相应的GFF文件。

发表于 更新于 分类于 ，

可能让大家失望了，我的主力电脑不是服务器平台，也不是HEDT（high end desktop）平台。得益于技术进步，最普通的家用电脑，也能很好的完成我的生物信息分析。

配置单

1
2
3
4
5
6
CPU：英特尔 i5 10400
主板：华擎B460M Pro4
内存：光威 DDR4 2666 32G * 4 = 128G
硬盘：西部数据 SN550 1TB，西部数据紫盘 4T
电源：金河田 金牌 500W
机箱：先马 小黑洞

这个配置的价格约6000元，性价比非常高。

使用体验

盒装CPU自带的散热器，全铝鳍片，没有导热铜管，下压式而不是侧吹式，我担心因为散热的问题不能发挥出CPU的全部性能。实际测试，盒装自带的散热器完全够用，i5-10400能跑满全核心睿频4.0GHz，温度在62度左右。

双通道内存有15%左右的性能提升，多线程有12%左右的性能提升。

查看CPU运行频率 

1
2
sudo apt install cpufrequtils
cpufreq-info |grep 'current CPU frequency'

查看CPU温度 

1
2
sudo apt install lm-sensors
sensors

消费级CPU主频高、单核性能好，再配上NVMe固态硬盘，大部分时候的使用体验优于组里的戴尔T640服务器。目前消费级平台最大只支持128G内存，HEDT平台一般可以支持到256G内存，需要更高的内存的话，就需要上服务器平台了。128G内存除了跑大转录组和大基因组的从头组装不行外，其他任务都能胜任。一块4T硬盘不够用，之后又加了一块4T硬盘，以后还需要再加。所有的数据都不能删，硬盘存储是最大的压力。

0%