由于分析需求需要使用PhastCons分析DNA的保守序列，以及使用SIFT软件分析氨基酸变异的有害性，经过一段时间的探索，对这两款软件已经有了还算全面的认识，在此记录一下分析流程。
使用PhastCons是因为它的同类软件GERP++，但是通过搜索各类文档，感觉并没有特别详尽的教程，官方文档有些地方也有点不太明白，使用了一部分小数据按照官方介绍的流程进行了分析，但是并没有跑出结果来，可能还是我准备的文件不行，于是弃用了这个软件，转而使用它的同类软件PhastCons，PhastCons有比较详尽的使用说明（PhastCons HOWTO），然后也有做过此类分析的教程。
SIFT软件刚开始找的时候大部分都是介绍使用它的网页版来对人类的基因组进行分析，而人类的基因组都是使用这个软件计算好的，我需要的是从头来构建一个自己研究物种的有害性变异的数据库，最后按照Github上的说明文档，最终构建成功，并可用于分析，也把这个软件的流程记录到这里。
PhastCons
1，物种进化树
全基因组多序列比对的主要目的是揭示基因组中进化上保守的和非保守的区段以及它们的分布。为了衡量这种保守性，可以通过引入各种不同的算法（反正也看不懂）来给定基因组的区段或碱基一个数值，这样通过比较基因组其他区域与感兴趣区域的保守性score，就可以了解这个区段是否“很保守”或者“很不保守”。
PhastCons和PhyloP能够把由MULTIZ产生的多序列比对结果转换成两种单一的保守性计分和显示，PhyloP只考虑比对的当前列而PhastCons同时也考虑比对列的相邻列，这使得PhastCons对于保守区段的出现更敏感，而PhyloP对于保守区段的界定更有效。PhyloP和PhastCons之间的另一个区别是，PhyloP能够识别加速进化和保守进化的位点，它们分别产生“正”的和“负”的计分，而PhastCons的计分总是一个0至1之间的正值。BioChen 博客
既然是基于已知的种系结构，那就需要获得相应的物种进化树，进化树的构建参考之前的文章单拷贝直系同源基因构建物种进化树，构建好进化树之后一定要跟之前已经发表的文章做一下对比，不要犯不符合先验条件的错误。
2，多基因组比对
2.1 使用RepeatMasker屏蔽基因组序列
去 Ensembl Plants 下载基因组的时候发现每个基因组都有对应的遮蔽（mask）掉重复序列的文件，标注 rm 的序列文件属于硬遮蔽文件，重复序列以“N”表示，标注 sm 为软遮蔽文件，重复序列以小写的“atcg”表示，两种都可以用于分析，我这里下载了 rm 序列作为多基因组比对的分析文件。
然后也用RepeatMasker跑了一下水稻的，并做了一下比较
并没有使用RepeatModeler重新预测，直接用的植物的TE库mipsREdat_9.3p_Poaceae_TEs.fasta
source /usrdata/users/hwwang/miniconda3/bin/activate /usrdata/users/hwwang/miniconda3/envs/rep/
RepeatMasker ./Oryza_sativa.IRGSP-1.0.dna.toplevel.chr.fa -e rmblast -qq -gff -pa 6 --lib ~/hychao/Function_anno/TE_db/mipsREdat_9.3p_Poaceae_TEs.fasta 2>runRice.log
seqkit stats -a -G "N" ./Oryza_sativa.IRGSP-1.0.dna.toplevel.chr.fa.masked
在这里使用了 -qq 参数，属于”急速“（rush）模式，所以在准确率上会差一些，不过超快，我觉得选择哪种模式取决于你的课题需要，如果仅是屏蔽序列，让基因组比对更快一些，这样就足够了。最后用seqkit统计下，40%重复序列，而下载的屏蔽基因组序列为50%。
以后再试试慢速模式下的屏蔽效果
既然有现成的屏蔽好的基因组序列文件，那就不用再费时间屏蔽了。
2.2 两两基因组比对
先下载软件
cd ~/hychao/biosoft
git clone https://github.com/hillerlab/GenomeAlignmentTools.git
cd GenomeAlignmentTools/kent/src/
make
echo "export PATH=$HOME/hychao/biosoft/GenomeAlignmentTools/kent/bin:\$PATH" >> ~/.bashrc
cd ../bin
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/chainNet
chmod +x ./*
export KENTSRC_DIR=$HOME/hychao/biosoft/GenomeAlignmentTools/kent/src
export MACHTYPE=x86_64
cd ~/hychao/biosoft/GenomeAlignmentTools/src/
make
cp *.py ../bin/
cp *.perl ../bin/
echo "export PATH=$HOME/hychao/biosoft/GenomeAlignmentTools/bin:\$PATH" >> ~/.bashrc
cd ~/hychao/biosoft
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 ~/hychao/biosoft/GenomeAlignmentTools/bin/
source ~/.bashrc
设置一个参考物种，用参考物种分别与其他多个物种进行比对
这里就用小麦的ABD基因组作为参考，将小麦的ABD基因组拆分，然后用ABD基因组分别与多个物种进行基因组的比对
samtools faidx Triticum_aestivum.IWGSC.dna_rm.toplevel.fa
less Triticum_aestivum.IWGSC.dna_rm.toplevel.fa.fai|cut -f1|grep A > a.chr
less Triticum_aestivum.IWGSC.dna_rm.toplevel.fa.fai|cut -f1|grep B > b.chr
less Triticum_aestivum.IWGSC.dna_rm.toplevel.fa.fai|cut -f1|grep D > d.chr
ls *chr|while read file; do seqkit grep -f $file Triticum_aestivum.IWGSC.dna_rm.toplevel.fa|seqkit seq -i > $file.rm.fa;done
以A为例吧
mkdir a_fa b_fa d_fa
faSplit byName a.chr.rm.fa a_fa
mv *A.fa a_fa
faToTwoBit a.chr.rm.fa a.chr.rm.fa.2bit
faToTwoBit Oryza_sativa.chr.rm.fa Oryza_sativa.chr.rm.fa.2bit
faSize a.chr.rm.fa -detailed > a.chr.rm.fa.sizes
faSize Oryza_sativa.chr.rm.fa -detailed > Oryza_sativa.chr.rm.fa.sizes
cd a_fa
ls *fa|cut -d "." -f1|awk '{print "lastz "$1".fa ~/hychao/data/zhao/2.gerp/masker_genome/rice/Oryza_sativa.chr.rm.fa O=400 E=30 K=3000 L=3000 H=2200 T=1 --format=axt --ambiguous=n --ambiguous=iupac > "$1".axt"}' > job
split -l1 job
ls xa*|awk '{print "jsub sh "$1""}' > jsub.sh
sh jsub.sh
2.3 使用MULTIZ合并两两基因组比对结果
MULTIZ不是一个序列比对软件。它利用基因组两两比对的结果，结合它们的进化树，得到多基因组比对结果。比如，参考物种ref_species和species_1、species_2的多基因组比对，先分别做ref_species和species_1比对、ref_species和species_2比对，然后用两两比对的结果和ref_species、species_1、species_2的进化树作为输入，利用MULTIZ得到多基因组比对结果。
最后得到maf格式的多基因比对结果
软件安装
cd ~/hychao/biosoft
wget https://github.com/multiz/multiz/archive/20190527.tar.gz
tar zxf 20190527.tar.gz
cd multiz-20190527/
make
echo "export PATH=$HOME/hychao/biosoft/multiz-20190527:\$PATH" >> ~/.bashrc
source ~/.bashrc
上一步用小麦的ABD基因组与多个物种进行了基因组的比对，这里以小麦B基因组为例，以小麦的B基因组作为ref，水稻和短柄草的基因组序列作为query：
将两两基因组比对得到的maf文件改为multiz的标准输入文件名，然后运行生成的脚本
ln -s CS_B.Br.maf CS_B.Br.sing.maf
ln -s CS_B.Os.maf CS_B.Os.sing.maf
roast - T=`pwd` E=CS_B "((CS_B Br) Os)" *.*.maf ref_species_mulitway.maf > roast.sh
bash roast.sh
这里说一下多序列比对的时间，我设置了多个参考物种，对于短柄草和水稻这种基因组比较小的物种，比对到4G的小麦亚基因组上，如果多线程并行跑的话七八个小时以内可以运行完毕，但是我也设置了大麦这种基因组相对较大的物种，运行时间要好几天，生成的比对文件也是很大的。
3， 使用PhastCons进行分析
3.1 软件安装
PHAST依赖CLAPACK，需要先编译CLAPACK，然后编译PHAST时指定CLAPACK的路径
# 安装CLAPACK
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
wget https://github.com/CshlSiepelLab/phast/archive/v1.5.tar.gz
tar zxf v1.5.tar.gz
cd phast-1.5/src/
make CLAPACKPATH=/usrdata/users/hwwang/hychao/biosoft/CLAPACK-3.2.1
# 将PHAST添加到环境变量
echo "export PATH=$HOME/hychao/biosoft/phast-1.5/bin:\$PATH" >> ~/.bashrc
source ~/.bashrc
#UCSC Tools
cd ~/hychao/biosoft/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 *
3.2 进化分析
使用上一步生成的mulitway.maf文件
运行phyloFit生成MOD文件
phyloFit -i MAF CS_B_mulitway.maf
生成phyloFit.mod文件，mod文件是根据多基因组比对的结果生成的一个进化树文件，可以采用多种方法来生成进化树，比如4d sites啥的，官方文档上都有说明，这里就按照正常的分析流程，然后检查了一下生成的进化树的拓扑结构是没有问题的。
3.3 计算保守性分数
phastCons一次只能处理一条染色体，先使用mafSplit分割MAF文件。
mkdir maf
mafSplit _.bed maf/ CS_B_mulitway.maf -byTarget -useFullSequenceName
使用-byTarget参数分割MAF，第一个参数_.bed会被忽略掉，但是必须要有这个参数占位。-useFullSequenceName参数的作用是生成的文件名是染色体名。
接下来使用phastCons计算每个碱基的phastCons分数和保守的区块。
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 >> CS_B_mulitway.wig
cat bed/*_most-cons.bed >> CS_B_most-cons.bed
CS_B_mulitway.wig文件有全基因组单碱基phastCons分数，CS_B_most-cons.bed文件有保守元件位置。
3.4 进一步分析
- 使用
wigToBigWig将wig文件转换为二进制的bigWig文件后，就能在Genome Browser中展示了。
wigToBigWig in.wig chrom.sizes out.bw
2.假设给定一条序列比如lncRNA，我们想计算它的平均phastCons分数，可以使用bigWigAverageOverBed实现。
bigWigAverageOverBed in.bw in.bed out.tab
- 保守元件去掉编码区的之后，就得到了保守的非编码元件（Conserved noncoding DNA elements，CNEs）。
3.5 使用snpSIFT分析变异位点的保守性
在上面的分析结果中，我们得到了单个染色体的wig文件，批量将这些文件改名为chr*.phastCons.wigFix的形式，这是snpSIFT的输入格式，然后gzip压缩一下，为每个染色体生成chr*.phastCons.wigFix.gz的文件。
将maf文件夹改个名称，mv maf snpeff_PhastCons
之后运行命令：
java -jar SnpSift.jar phastCons ./snpeff_PhastCons your.vcf > your.phastCons.vcf
结果文件为
在注释区域的最后加上了一栏PhastCons分数。
SIFT4G
1 软件安装
先根据官方流程下载构建数据库的一系列脚本GitHub - SIFT4G_Create_Genomic_DB
然后conda 安装sift4g conda create -n sift, conda activate sift, conda install sift4g, conda install gcc_impl_linux-64
下面的依赖安装代码是参考别人的SIFT4G:预测氨基酸取代是否会影响蛋白质功能，开始没有使用conda安装，所以这些依赖全部都安装上了，但是运行的时候报错，说是迭代错误，查了下应该是gcc的问题，所以直接用conda新建了一个环境成功运行，至于conda安装时下面的依赖是否也配置好了不得而知。如果conda使用过后仍然报错，那应该就是下面的依赖还没有配置好。
#DBI
perl -MCPAN -e "install Bundle::LWP"
perl -MCPAN -e "install HTML::Formatter"
perl -MCPAN -e "install IO::Socket::SSL and OpenSSL"
#bioperl
wget https://cpan.metacpan.org/authors/id/C/CJ/CJFIELDS/BioPerl-1.7.7.tar.gz
tar -zxvf BioPerl-1.7.7.tar.gz
cd BioPerl-1.7.7
perl Makefile.PL
make
make test
make install
#安装Switch.pm
cpanm Switch
#下载安装DBI
wget -C http://www.cpan.org/modules/by-module/DBD/DBI-1.643.tar.gz
tar -zxvf DBI-1.643.tar.gz
cd DBI-1.643
perl Makefile.PL
make
make test
make install
2 运行
根据官方教程里的内容，配置好config文件，test_file里面有人类基因组构建的流程和相应的配置文件，我们只需要复制下人类的config文件，然后把里面的小部分内容修改成自己的即可。
cd ~/software/scripts_to_build_SIFT_db
vi work.sh
###
source /usrdata/users/hwwang/miniconda3/bin/activate /usrdata/users/hwwang/miniconda3/envs/sift/
perl make-SIFT-db-all.pl -config ~/hychao/data/zhao/3.sift4G_db/a_subgenome/CS_A.txt 2>run_CS_A.log
###
jsub sh work.sh
上面的流程是一站式搞定的
建库的时候除了准备config文件外，还需要建两个文件夹
gene-annotation-src和chr-src，分别放置gtf.gz文件和分染色体的序列文件
3 sift注释
从这里下载 SIFT4G_Annotator，使用起来很简单：
java -jar ~/hychao/biosoft/SIFT4G_Annotator-master/SIFT4G_Annotator.jar -c -i ./output_A.vcf -d ~/hychao/data/zhao/3.sift4G_db/a_subgenome/CS_V1.0/ -r ./
运行完之后会按照染色体分别生成两个文件.tsv, .vcf文件，里面记录了氨基酸的变化，SIFT score，有害性评估等信息。
