写在前面
没想到，上个选修课，准备示例数据和镜像上，竟然花了我三四天....主要原因还是，「网速太差」。搞七搞八，终于搞定。尽管整完之后，我其实担心学生拿到这个使用文档也不知道从而下手。尽管，我认为，有前面四次实践经验，这一次应该比较简单。最近基本停更了，此次就用这个实践代码作为推文。具体覆盖了比较基因组分析的大部分内容。
- 基因家族聚类和物种树获取
- 物种分歧时间
- 基因家族扩展与收缩分析
- WGD分析（基于Ks分布，已覆盖基于 dn/ds 的基因选择分析）
- 基因组共线性分析与WGD分析（基于点图）
构建镜像
拿到 Demo Data 之后，解压，第一步是构建一个镜像。与前述相同，镜像只需要构建一次。
docker build -t section5 .
实践命令
与前面的命令操作不同，本次实践，直接进去docker环境进行数据分析
docker run -it -v C:\Users\CJ\Desktop\园艺植物生物信息学实践\课件\第五周\Demodata:/scauclass -w /scauclass section5
定义主工作目录
currWd=`pwd`
基因家族聚类和物种树获取
基于蛋白序列集合（目录下数个物种的蛋白序列），进行基因家族鉴定；注意电脑线程数，另经测试，电脑至少要有 8Gb 内存
mkdir proteinMCL
cp pepSimID/*.fa proteinMCL
orthofinder -f ./proteinMCL -t 6
获取最新 orthofinder 结果路径
mclPath=`pwd`/`find proteinMCL -name "Log.txt"|tail -n 1|cut -d/ -f1-3`
查看单拷贝基因的结果，也可以在本地自行找到对应目录，查看结果
ls $mclPath/Single_Copy_Orthologue_Sequences/
查看物种树结果
head $mclPath/Species_Tree/SpeciesTree_rooted*
物种特有家族鉴定
可自行查看输出结果，自行理解下述图中表述。
- 对于物种特有家族，查看 Orthogroup.GeneCounts.tsv，进行简单筛选，随后从 Orthogroups.txt 中提取对应物种的序列即可；
- 对于物种特有且单拷贝基因，查看 Orthogroups_UnassignedGenes.tsv
物种分歧时间
准备数据
mkdir MCMCtree
cp $mclPath/Orthogroups/Orthogroups_SingleCopyOrthologues.txt MCMCtree
cp $mclPath/Orthogroups/Orthogroups.txt MCMCtree
cp $mclPath/Species_Tree/SpeciesTree_rooted.txt MCMCtree
对每一个直系同源单拷贝基因进行多序列别对
cd $currWd/MCMCtree
练习缘故，只取50个序列进行比对
head -n 50 Orthogroups_SingleCopyOrthologues.txt|while read og;do
seqkit grep -f <(grep $og Orthogroups.txt|perl -lane 'shift @F;print for @F') $currWd/cdsSimID/*.fa > $og.cds.fa
perl -i -lpe '/(>\w{3})_/ and $_=$1' $og.cds.fa
muscle -in $og.cds.fa -physout $og.phy
rm $og.cds.fa
done
合并比对结果
cat OG*.phy > merged.phy
rm OG*.phy
去除枝长信息
cat SpeciesTree_rooted.txt
(Ath:0.176983,((Ppe:0.136384,Fve:0.238589)0.924762:0.0994348,Csi:0.238093)1:0.176983);
直接手动编辑去除
(Ath,((Ppe,Fve),Csi));
如果不知道如何编辑去除，那么在练习上，直接运行
echo "(Ath,((Ppe,Fve),Csi));" > SpeciesTree_rooted.txt
到 timetree 网站查找指定分化时间
http://www.timetree.org/
Orange
Strawberry
补充对应的矫正信息，更详细的时间矫正信息，自行参考 mcmctree 文档，在本例中，可以直接运行
# 一共 4 个物种
echo "4 1" > timetree.nwk
# 拟南芥和甜橙的分化时间大体是 96~104 MYA
echo "(Ath,((Ppe,Fve),Csi)'>0.98<1.17');" >> timetree.nwk
cat timetree.nwk
准备配置文件mcmctree.ctl
cp $currWd/mcmctree.calcBV.ctl .
其中内容如下（具体可根据个人需要调整）
seed = -1
seqfile = merged.phy
treefile = timetree.nwk
outfile = out
ndata = 20
seqtype = 0 * 0: nucleotides; 1:codons; 2:AAs
usedata = 3 * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
clock = 3 * 1: global clock; 2: independent rates; 3: correlated rates
RootAge = '<2.0' * safe constraint on root age, used if no fossil for root.
model = 4 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
alpha = 0.5 * alpha for gamma rates at sites
ncatG = 5 * No. categories in discrete gamma
cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)?
BDparas = 1 1 0 * birth, death, sampling
kappa_gamma = 6 2 * gamma prior for kappa
alpha_gamma = 1 1 * gamma prior for alpha
rgene_gamma = 2 2 * gamma prior for overall rates for genes
sigma2_gamma = 1 10 * gamma prior for sigma^2 (for clock=2 or 3)
finetune = 1: 0.1 0.1 0.1 0.01 .5 * auto (0 or 1) : times, musigma2, rates, mixing, paras, FossilErr
print = 1
burnin = 2000
sampfreq = 10
nsample = 20000
随后，可以执行命令
mcmctree mcmctree.calcBV.ctl
基于序列信息获得BV后，使用信息用近似计算法估算物种分歧时间
mv out.BV in.BV
cp $currWd/mcmctree.useBV.ctl .
mcmctree mcmctree.useBV.ctl
很快就计算完了
可以在输出目录看到
使用 Figtree 软件，打开查看即可
Bar比较大，估计还是跟输入数据多少以及参数有关系，可能使用所有序列，或者调整合适的参数，那么会得到准确度和精度都更高的结果。
基因家族扩展与收缩分析
export LD_LIBRARY_PATH="/opt/conda/lib"
开始分析
cd $currWd
mclPath=`pwd`/`find proteinMCL -name "Log.txt"|tail -n 1|cut -d/ -f1-3`
mkdir geneExpansion
cd $_
cp $mclPath/Orthogroups/Orthogroups.GeneCount.tsv .
perl -F'\t' -lane 'if($.==1){$flag=qq{Desc}}else{$flag=null}pop @F;print join qq{\t},$flag,@F;' Orthogroups.GeneCount.tsv > gene.counts.tab
cp $currWd/MCMCtree/FigTree.tre .
perl -lne 'next unless /UTREE/;s/^.*?(\(.*\)).*?$/$1/;s/\[.*?\]//g;s/\s+//g;print $_,q{;}' FigTree.tre > tree.withTime.nwk
cafe5 --infile gene.counts.tab --tree tree.withTime.nwk
查看多少个家族有显著扩张收缩
perl -lane '$count++ if $F[2] eq qq{y};END{print $count}' results/Base_family_results.txt
683
查看家族扩张与收缩相关的进化树信息
head results/Base_asr.tre
提取统计学上有显著扩张收缩的家族，可具体查看自己感兴趣的家族
grep "*" results/Base_asr.tre|perl -lane 'print $F[1]' > og.sig.group.ids
查看每个节点的家族扩张和收缩情况
head results/Base_clade_results.txt
基于 Ks 分布判断 WGD
cd $currWd/
mkdir WGD
cd WGD/
cp $currWd/cdsSimID/Ath.fa Ath.cds.fa
cp $currWd/pepSimID/Ath.fa Ath.pep.fa
diamond makedb --db Ath.pep.fa --in Ath.pep.fa
diamond blastp --db Ath.pep.fa --query Ath.pep.fa --outfmt 6 --out Ath.pep.fa.diamond
过滤不要保留太多
perl -lane 'next if $F[0] eq $F[1];next if $dup{qq{$F[0]-$F[1]}}++ || $dup{qq{$F[1]-$F[0]}}++;next if $seen{$F[0]}++ > 5;print' Ath.pep.fa.diamond > Ath.pep.fa.filtered.diamond
整理获得基因对
java -cp $currWd/beans/TBtools_JRE1.6.jar biocjava.bioIO.BioSoftPipeServer.PairWiseKaKsCalculator --inCDS Ath.cds.fa --inGenePair Ath.pep.fa.filtered.diamond --outKaks Ath.pep.filtered.diamond.kaks.tab.xls --inCPU 6
提取合理的 Ks 区间，一般是 0~4
perl -lane 'print $F[3] if $F[3]>0 && $F[3]<4' Ath.pep.filtered.diamond.kaks.tab.xls > Ks.filtered.tab.xls
在本地使用 Excel 绘制 直方图 即可，可检测到两次 WGD
染色体来源
（复制；重组；变异；祖先…）
直接使用示例数据，进行分析，感兴趣的同学请按文件格式自行准备相应文件
cd $currWd/Synteny
MCScanX os_sb
可以绘制相应图片，可以从点图看出全基因组复制事件（PS: 少数平台可能绘制不了，不影响，直接看图吧）
java -cp /MCScanX-master/downstream_analyses dot_plotter -g os_sb.gff -s os_sb.collinearity -c dot.ctl -o dotplot.png
环境准备过程 - 无需实践
初步搭建镜像
docker pull continuumio/miniconda3:latest
docker run -it -v C:\Users\CJ\Desktop\园艺植物生物信息学实践\课件\第五周\Demodata:/scauclass -w /scauclass continuumio/miniconda3:latest
conda install -y -n base conda-forge::mamba
mamba install -c conda-forge gcc gxx
mamba install -y make
编译 cafe5
tar -zxvf CAFE.tar.gz
cd CAFE/
./configure
make -j 8
mamba install -y -c bioconda orthofinder
mamba install -y -c bioconda paml
mamba install -y -c bioconda seqkit diamond
# 不知道会不会和 orthofinder 或者其他软件冲突？要考量
mamba install -c bioconda muscle=3.8
mamba install openjdk=11
编译安装 MCScanX
unzip MCScanX-master.zip
cd MCScanX-master
make -j 2
Dockerfile 的编写 - 无需实践
FROM continuumio/miniconda3:latest
RUN conda install -y -n base conda-forge::mamba
RUN mamba install -y -c conda-forge gcc gxx
RUN mamba install -y make
RUN mamba install -y unzip
## install CAFE5
COPY CAFE5.tar.gz /CAFE5.tar.gz
RUN tar -zxvf CAFE5.tar.gz
RUN cd /CAFE5 && ./configure && make -j 2
##
RUN mamba install -y -c bioconda orthofinder
RUN mamba install -y -c bioconda paml
RUN mamba install -y -c bioconda seqkit diamond
RUN mamba install -c bioconda muscle=3.8
RUN mamba install openjdk=8
## install MCScanX
COPY MCScanX-master.zip /MCScanX-master.zip
RUN unzip MCScanX-master.zip
RUN cd /MCScanX-master && make -j 2
ENV PATH=/CAFE5/bin/:/MCScanX-master:${PATH}
