歡迎關(guān)注R語言數(shù)據(jù)分析指南
最近在做一個基因家族的數(shù)據(jù)分析陨囊，整理了一下以前的筆記，本節(jié)通過水稻的案例來介紹如何計算該物種旁系同源基因的ka/ks值
軟件安裝
conda install clustalw
conda install blast
數(shù)據(jù)下載
http://rice.uga.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/
除了上述2款軟件外兔毒，還需要mcscanx宙橱，ParaAT,KaKs_Calculator三款軟件馋袜，需要自己官網(wǎng)下載編譯，并添加環(huán)境變量
1.整理gff文件
第一步需要根據(jù)自己研究的物種gff文件將其進(jìn)行整理球散，不同的物種gff文件略有不同,下面以水稻為例子
library(tidyverse)
read_tsv("Oryza_sativa.IRGSP.gff",comment = "#",col_names = F) %>%
filter(X3=="mRNA") %>%
mutate_at(vars(c(X9)),~str_split(.,";",simplify=T)[,1]) %>%
mutate(across("X9",str_replace,"ID=:","")) %>%
select(1,9,4,5) %>%
write.table(.,file="Os.gff",sep="\t",quote = F,row.names = F,col.names = F)
- 整理成如下四列的格式
Chr1 LOC_Os01g01010.1 2903 10817
Chr1 LOC_Os01g01010.2 2984 10562
Chr1 LOC_Os01g01019.1 11218 12435
Chr1 LOC_Os01g01030.1 12648 15915
Chr1 LOC_Os01g01040.1 16292 20323
Chr1 LOC_Os01g01040.2 16321 20323
Chr1 LOC_Os01g01040.3 16321 20323
Chr1 LOC_Os01g01040.4 16292 18304
Chr1 LOC_Os01g01050.1 22841 26971
2.對蛋白序列做blast比對
# 對蛋白序列構(gòu)建索引
makeblastdb -in data/protein.fa -dbtype prot -title protein.fa
# 進(jìn)行序列比對
blastp -query ~/data/protein.fa -out Os.blast -db ~/data/protein.fa -outfmt 6 -evalue 1e-5 -num_threads 30
3.mcscanx進(jìn)行共線性分析
注：兩個文件放到同一個文件夾中 Os.blast & Os.gff，通過mcscanx找到同源基因
mcscanx ./kaks/Os
- mcscanx輸出.tandem的文件需要按置表符分為2列
sed -i s'/,/\t/g' *.tandem
4.計算水稻全部ka/ks
- 輸出寫絕對路徑散庶，輸出文件夾無需先創(chuàng)建蕉堰，需要先整理CDS序列與蛋白序列格式
- 需要在ParaAT中執(zhí)行 chmod 777 Epal2nal.pl
數(shù)據(jù)整理
touch proc;echo 30 > proc # 設(shè)置線程數(shù)
# 將序列多行變一行
awk '/^>/&&NR>1{print "";}{printf "%s",/^>/?$0"\n":$0}' ../data/cds.fa > cds.fasta
awk '/^>/&&NR>1{print "";}{printf "%s",/^>/?$0"\n":$0}' ../data/protein.fa > protein.fasta
計算ka/ks值
ParaAT.pl -h ./kaks/Os.tandem -n ./cds.fasta -a ./protein.fasta -m clustalw2 -p proc -f axt -g -k KaKs_Calculator -o ./ka-ks
合并ks/ks文件
find . -name "*.kaks " -exec cat '{}' ';' > Os.kaks
less Os.kaks|sort|uniq > Os.rename.kaks.xls
數(shù)據(jù)獲取
根據(jù)個人電腦性能需要耗費(fèi)時間不同，blast比對 & ka/ks計算此兩步最耗費(fèi)時間悲龟；一般對于小基因組物種個人電腦需要耗費(fèi)幾十小時不等屋讶；今天的介紹到此結(jié)束
