比較基因組學(xué)分析目錄
1：單拷貝基因構(gòu)建物種樹以及計(jì)算分化時(shí)間
2：基因家族收縮與擴(kuò)張分析
3：特異節(jié)點(diǎn)富集分析
之前有寫過單拷貝基因構(gòu)建物種樹的流程，現(xiàn)在我對(duì)流程進(jìn)行優(yōu)化，而且將會(huì)增加后續(xù)的基因家族收縮與擴(kuò)張分析，希望對(duì)大家的分析有所幫助，如果文章被某些原因隱藏，請(qǐng)去我的博客
這三篇教程走下來，最終會(huì)生成下面的結(jié)果，希望大家也能掌握比較基因組學(xué)基本的分析技巧
1. 物種選擇
在此我不在介紹軟件的原理和安裝了，大家可以去看我之前的推文
物種選擇此次包含了整個(gè)綠色植物，從綠藻門到被子植物
大部分基因組是在jgi phytozome上下載的
綠藻門Chlorophyte：小球藻：Coccomyxa subellipsoidea(Csu), 團(tuán)藻：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)
下載完成后我對(duì)每個(gè)蛋白序列的ID前面加上了物種簡寫, 方便后續(xù)使用
sed -i "s/>/>Ath/g" Ath.pro.fasta
2.Orthofinder獲得單拷貝基因
將所有物種的蛋白文件放到CGA文件夾下
orthofinder -t 16 -f CGA/
我們主要使用Orthoogroups查看正交群的基因和使用 Single_Copy_Orthologue_Sequences里的單拷貝基因構(gòu)建系統(tǒng)發(fā)育樹
WorkingDirectory其中包含運(yùn)算過程的中間文件，例如diamond結(jié)果，如果我們想去掉某一物種，在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
如果想增加額外的物種進(jìn)行分析，可以按照如下方式運(yùn)行
orthofinder -b WorkingDirectory -f new_fasta_directory
3. 利用單拷貝基因構(gòu)建系統(tǒng)發(fā)育樹
Orthofinder的 Single_Copy_Orthologue_Sequences下存放著單拷貝同源基因的序列，我們可以利用這些序列構(gòu)建系統(tǒng)發(fā)育樹
建樹方法可以分為串聯(lián)法與并聯(lián)法
不同的是，串聯(lián)法是將得到的單拷貝同源基因比對(duì)后進(jìn)行了串聯(lián)，每個(gè)物種都得到一個(gè)很大的序列，然后進(jìn)行建樹；并聯(lián)法是將每個(gè)單拷貝同源基因集比對(duì)后建樹，然后再利用Astral構(gòu)建了物種樹，這里目前只講串聯(lián)法，并聯(lián)法正在研究
3.1 串聯(lián)法建樹
vi test.sh
for i in *.fa
do
file=${i%.fa*}
mafft --maxiterate 1000 --localpair "${file}.fa" > "${file}.aln" ## 多序列比對(duì)
(也可以添加一步，將蛋白文件轉(zhuǎn)為cds文件再進(jìn)行后續(xù)的分析，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)請(qǐng)自行安裝
接下來將所有的aln文件串聯(lián)
seqkit concat *.final.aln > all.fa
將所有單拷貝基因串聯(lián)在一起后，進(jìn)行模型預(yù)測，這次使用modeltestng進(jìn)行預(yù)測
modeltest-ng -i all.fa -d aa -o modeltest -p 16
在modeltest.out文件中，我們可以看到我們需要的模型并且可以直接copy命令
接下來進(jìn)行建樹
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文件查看樹結(jié)構(gòu)
結(jié)果和進(jìn)化關(guān)系是一致的
3.2 并聯(lián)法建樹
并聯(lián)法建樹使用Astral，將所有的單拷貝基因樹合并為一顆物種樹
orthDir=~/CGApaper/protein/OrthoFinder/Results_May02 #基于Orthofinder結(jié)果
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，里面包含了不同拓?fù)浣Y(jié)構(gòu)的可能性，目前先到這一步
(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. 物種分化時(shí)間計(jì)算
此次使用mcmctree計(jì)算物種的分化時(shí)間，首先在timetree上查詢物種的分化時(shí)間進(jìn)行標(biāo)定，修改樹文件，需要注意在標(biāo)定時(shí)間是盡量較早的節(jié)點(diǎn)和較晚的節(jié)點(diǎn)都進(jìn)行標(biāo)定
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 ##此步驟先設(shè)為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.
之后運(yùn)行
mcmctree mcmctree.ctl
運(yùn)行結(jié)果產(chǎn)生一個(gè)文件名為 out.BV，我們拷貝到兩次分化時(shí)間估計(jì)的目錄里，并用in.BV 作為新的命名
cp out.BV ../02.Time/rep1/in.BV
cp out.BV ../02.Time/rep2/in.BV
之后在運(yùn)行兩次重復(fù)
首先修改配置文件
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.
最后我們會(huì)得到如下文件
all.fa all.phy all.tre FigTree.tre in.BV mcmctree.ctl mcmc.txt out SeedUsed
檢查兩次結(jié)果的FigTree.tre文件，如果相差不大就可以使用
之后可以進(jìn)行基因家族的收縮與擴(kuò)張分析了，這部分之后在寫
如果覺得有幫助歡迎在博客打賞
