写在前面没想到，上个选修课，准备示例数据和镜像上，竟然花了我三四天....主要原因还是，「网速太差」。搞七搞八，终于搞定。尽管整完之后，我其实担心学生拿到这个使用文档也不知道从而下手。尽管，我认为，有前面四次实践经验，这一次应该比较简单。最近基本停更了，此次就用这个实践代码作为推文。具体覆盖了比较基因组分析的大部分内容。 基因家族聚类和物种树获取 物种分歧时间 基因家族扩展与收缩分析 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
定义主工作目录 基因家族聚类和物种树获取基于蛋白序列集合（目录下数个物种的蛋白序列），进行基因家族鉴定；注意电脑线程数，另经测试，电脑至少要有 8Gb 内存 mkdir proteinMCLcp 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 MCMCtreecp $mclPath/Orthogroups/Orthogroups_SingleCopyOrthologues.txt MCMCtreecp $mclPath/Orthogroups/Orthogroups.txt MCMCtreecp $mclPath/Species_Tree/SpeciesTree_rooted.txt 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.faperl -i -lpe '/(>\w{3})_/ and $_=$1' $og.cds.famuscle -in $og.cds.fa -physout $og.phyrm $og.cds.fadone
合并比对结果 cat OG*.phy > merged.phyrm 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);
直接手动编辑去除 如果不知道如何编辑去除，那么在练习上，直接运行 echo "(Ath,((Ppe,Fve),Csi));" > SpeciesTree_rooted.txt
到 timetree 网站查找指定分化时间 http://www.timetree.org/OrangeStrawberry
补充对应的矫正信息，更详细的时间矫正信息，自行参考 mcmctree 文档，在本例中，可以直接运行 # 一共 4 个物种echo "4 1" > timetree.nwk# 拟南芥和甜橙的分化时间大体是 96~104 MYAecho "(Ath,((Ppe,Fve),Csi)'>0.98<1.17');" >> timetree.nwkcat 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.BVcp $currWd/mcmctree.useBV.ctl .
mcmctree mcmctree.useBV.ctl
很快就计算完了 
可以在输出目录看到
使用 Figtree 软件，打开查看即可 
Bar比较大，估计还是跟输入数据多少以及参数有关系，可能使用所有序列，或者调整合适的参数，那么会得到准确度和精度都更高的结果。 基因家族扩展与收缩分析export LD_LIBRARY_PATH="/opt/conda/lib"
开始分析 cd $currWdmclPath=`pwd`/`find proteinMCL -name "Log.txt"|tail -n 1|cut -d/ -f1-3`mkdir geneExpansioncd $_
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.tabcp $currWd/MCMCtree/FigTree.tre .perl -lne 'next unless /UTREE/;s/^.*?(\(.*\)).*?$/$1/;s/\[.*?\]//g;s/\s+//g;print $_,q{;}' FigTree.tre > tree.withTime.nwkcafe5 --infile gene.counts.tab --tree tree.withTime.nwk
查看多少个家族有显著扩张收缩 perl -lane '$count++ if $F[2] eq qq{y};END{print $count}' results/Base_family_results.txt683
查看家族扩张与收缩相关的进化树信息 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 分布判断 WGDcd $currWd/mkdir WGDcd WGD/
cp $currWd/cdsSimID/Ath.fa Ath.cds.facp $currWd/pepSimID/Ath.fa Ath.pep.fadiamond makedb --db Ath.pep.fa --in Ath.pep.fadiamond 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
染色体来源（复制；重组；变异；祖先…）
直接使用示例数据，进行分析，感兴趣的同学请按文件格式自行准备相应文件 可以绘制相应图片，可以从点图看出全基因组复制事件（PS: 少数平台可能绘制不了，不影响，直接看图吧） java -cp /MCScanX-master/downstream_analyses dot_plotter -g os_sb.gff -s os_sb.collinearity -c dot.ctl -o dotplot.png
|
