BioChen 博客
0%

发表于 更新于 分类于 ，

写在前面

我们做生信分析时，会用到很多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脚本，让大家方便地运行以上步骤。