比较基因组学分析目录
1：单拷贝基因构建物种树以及计算分化时间
2：基因家族收缩与扩张分析
3：特异节点富集分析
之前有写过单拷贝基因构建物种树的流程，现在我对流程进行优化，而且将会增加后续的基因家族收缩与扩张分析，希望对大家的分析有所帮助，如果文章被某些原因隐藏，请去我的博客
这三篇教程走下来，最终会生成下面的结果，希望大家也能掌握比较基因组学基本的分析技巧
1. 物种选择
在此我不在介绍软件的原理和安装了，大家可以去看我之前的推文
物种选择此次包含了整个绿色植物，从绿藻门到被子植物
大部分基因组是在jgi phytozome上下载的
绿藻门Chlorophyte：小球藻：Coccomyxa subellipsoidea(Csu), 团藻：Volvox carteri(Vca)
轮藻门Charophyta：布氏轮藻：Chara braunii(Cba), Klebsormidium nitens(Kle)
苔藓植物门Bryophyta：地钱：Marchantia polymorpha(Mpo), 小立碗藓：Physcomitrium patens(Ppa)
石松类Lycopsida： 江南卷柏：Selaginella moellendorffii(Smo)
蕨类Ferns： 水蕨：Ceratopteris richardii(Cri)
裸子植物门Gymnospermae：北美乔柏 Thuja plicata(Tpl)
基部被子植物Basial Angiosperm：无油樟：Amborella trichopoda(Atr),睡莲：Nymphaea colorata(Nco)
双子叶植物Eudicot：葡萄：Vitis vinifera(Vvi), 拟南芥：Arabidopsis thaliana(Ath), 黄瓜：Cucumis sativus (Csa)
单子叶植物Monocot: 芦笋：Asparagus officinalis(Aof), 水稻：Oryza sativa(Osa), 玉米：Zea mays(Zma)
下载完成后我对每个蛋白序列的ID前面加上了物种简写, 方便后续使用
sed -i "s/>/>Ath/g" Ath.pro.fasta
2.Orthofinder获得单拷贝基因
将所有物种的蛋白文件放到CGA文件夹下
orthofinder -t 16 -f CGA/
我们主要使用Orthoogroups查看正交群的基因和使用 Single_Copy_Orthologue_Sequences里的单拷贝基因构建系统发育树
WorkingDirectory其中包含运算过程的中间文件，例如diamond结果，如果我们想去掉某一物种，在SpeciesIDs.txt中将该物种注释掉
0: Aof.pro.fasta
1: Ath.pro.fasta
2: Atr.pro.fasta
3: Cba.pro.fasta
4: Cri.pro.fasta
5: Csa.pro.fasta
6: Csu.pro.fasta
7: Kle.pro.fasta
8: Mpo.pro.fasta
9: Nco.pro.fasta
10: Osa.pro.fasta
11: Ppa.pro.fasta
#12: Pum.pro.fasta
13: Smo.pro.fasta
14: Tpl.pro.fasta
15: Vca.pro.fasta
16: Vvi.pro.fasta
17: Zma.pro.fasta
orthofinder -b WorkingDirectory
如果想增加额外的物种进行分析，可以按照如下方式运行
orthofinder -b WorkingDirectory -f new_fasta_directory
3. 利用单拷贝基因构建系统发育树
Orthofinder的 Single_Copy_Orthologue_Sequences下存放着单拷贝同源基因的序列，我们可以利用这些序列构建系统发育树
建树方法可以分为串联法与并联法
不同的是，串联法是将得到的单拷贝同源基因比对后进行了串联，每个物种都得到一个很大的序列，然后进行建树；并联法是将每个单拷贝同源基因集比对后建树，然后再利用Astral构建了物种树，这里目前只讲串联法，并联法正在研究
3.1 串联法建树
vi test.sh
for i in *.fa
do
file=${i%.fa*}
mafft --maxiterate 1000 --localpair "${file}.fa" > "${file}.aln" ## 多序列比对
(也可以添加一步，将蛋白文件转为cds文件再进行后续的分析，pal2nal.pl "${file}.aln" "${file}.cds" -output fasta > "${file}.cds.aln)
Gblocks "${file}.aln" -b4=5 -b5=h -t=p -e=.gb ## 提取保守序列
seqkit seq "${file}.aln.gb" -w 0 > "${file}.new.aln" ## 多行序列并为一行
awk '{if (NR%2==1) print substr($1, 1, 4); else print $0}' "${file}.new.aln" > "${file}.final.aln" ## ID修饰
done
这一步需要的软件(mafft,Gblocks,seqkit)请自行安装
接下来将所有的aln文件串联
seqkit concat *.final.aln > all.fa
将所有单拷贝基因串联在一起后，进行模型预测，这次使用modeltestng进行预测
modeltest-ng -i all.fa -d aa -o modeltest -p 16
在modeltest.out文件中，我们可以看到我们需要的模型并且可以直接copy命令
接下来进行建树
raxmlHPC-PTHREADS-SSE3 -T 32 -f a -x 123 -p 123 -N 10000 -m PROTGAMMAILGF -k -O -o Csu,Vca -n all.tre -s all.fa
得到 RAxML_bipartitions.all.tre文件查看树结构
结果和进化关系是一致的
3.2 并联法建树
并联法建树使用Astral，将所有的单拷贝基因树合并为一颗物种树
orthDir=~/CGApaper/protein/OrthoFinder/Results_May02 #基于Orthofinder结果
cat $orthDir/Orthogroups/Orthogroups_SingleCopyOrthologues.txt |while read aa; do cat $orthDir/Gene_Trees/$aa\_tree.txt |awk '{print $0 }' ;done > SingleCopy.trees
sed -r 's/([(,]{1}[A-Za-z]+)_[^:]+/\1/g' SingleCopy.trees > Astral_input.trees
java -jar ~/soft/ASTRAL-5.7.1/astral.5.7.1.jar -i Astral_input.trees -o Astral_output.tree -t 8
最后得到Astral_output.tree，里面包含了不同拓扑结构的可能性，目前先到这一步
(Csu,(Vca,(Kle,(Cba,((Smo,(Cri,(Tpl,((Atr,Nco)'[q1=0.56;q2=0.24;q3=0.21]':0.3957605767189867,((Ath,(Vvi,Csa)'[q1=0.59;q2=0.24;q3=0.17]':0.47671301702525964)'[q1=0.9;q2=0.06;q3=0.04]':1.82319624265162,(Aof,(Zma,Osa)'[q1=0.98;q2=0.01;q3=0.01]':3.3756896221149395)'[q1=0.73;q2=0.11;q3=0.16]':0.8957898742947543)'[q1=0.59;q2=0.13;q3=0.28]':0.4793457920367108)'[q1=0.95;q2=0.01;q3=0.04]':2.471696418051194)'[q1=0.92;q2=0.02;q3=0.06]':2.040688555382602)'[q1=0.62;q2=0.18;q3=0.2]':0.5471348976356571)'[q1=0.37;q2=0.32;q3=0.32]':0.04487496622879318,(Ppa,Mpo)'[q1=0.74;q2=0.15;q3=0.1]':0.9286981407683164)'[q1=0.86;q2=0.07;q3=0.07]':1.5475625087160123)'[q1=0.56;q2=0.21;q3=0.23]':0.4056697116894409)'[q1=0.98;q2=0.02;q3=0.01]':3.121909101338843):0.0);
4. 物种分化时间计算
此次使用mcmctree计算物种的分化时间，首先在timetree上查询物种的分化时间进行标定，修改树文件，需要注意在标定时间是尽量较早的节点和较晚的节点都进行标定
17 1
((Kle,(((Mpo,Ppa),((Cri,(Tpl,(Atr,(Nco,(((Osa,Zma)'>42<52',Aof)'>104<125',(Vvi,(Csa,Ath))'>98<117')'>85<128')'>173<199')'>136>247')'>289<330')'>392<422',Smo)'410<468')'>472<515',Cba)),(Csu,Vca)'>800<1000')'>725<1200';
配置文件书写
cp ~/soft/paml4.9i/mcmctree.ctl mcmctree1.ctl
vi mcmctree1.ctl
seed = -1
seqfile = all.fa
treefile = all.tre
outfile = out
ndata = 1 ##必须注意为1
seqtype = 2 * 0: nucleotides; 1:codons; 2:AAs ##序列格式，蛋白就选择AA
usedata = 3 * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV ##此步骤先设为3
clock = 3 * 1: global clock; 2: independent rates; 3: correlated rates
RootAge = <1.0 * safe constraint on root age, used if no fossil for root.
model = 0 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
alpha = 0 * 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 = 20000
sampfreq = 2
nsample = 100000
*** Note: Make your window wider (100 columns) before running the program.
之后运行
mcmctree mcmctree.ctl
运行结果产生一个文件名为 out.BV，我们拷贝到两次分化时间估计的目录里，并用in.BV 作为新的命名
cp out.BV ../02.Time/rep1/in.BV
cp out.BV ../02.Time/rep2/in.BV
之后在运行两次重复
首先修改配置文件
seed = -1
seqfile = all.fa
treefile = all.tre
outfile = out
ndata = 1
seqtype = 2 * 0: nucleotides; 1:codons; 2:AAs
usedata = 2 * 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 = <1.0 * safe constraint on root age, used if no fossil for root.
model = 0 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
alpha = 0 * 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 = 800000
sampfreq = 40
nsample = 6500000
*** Note: Make your window wider (100 columns) before running the program.
最后我们会得到如下文件
all.fa all.phy all.tre FigTree.tre in.BV mcmctree.ctl mcmc.txt out SeedUsed
检查两次结果的FigTree.tre文件，如果相差不大就可以使用
之后可以进行基因家族的收缩与扩张分析了，这部分之后在写
如果觉得有帮助欢迎在博客打赏
