欢迎关注R语言数据分析指南
最近在做一个基因家族的数据分析，整理了一下以前的笔记，本节通过水稻的案例来介绍如何计算该物种旁系同源基因的ka/ks值
软件安装
conda install clustalw
conda install blast
数据下载
http://rice.uga.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/
除了上述2款软件外，还需要mcscanx，ParaAT,KaKs_Calculator三款软件，需要自己官网下载编译，并添加环境变量
1.整理gff文件
第一步需要根据自己研究的物种gff文件将其进行整理，不同的物种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比对
# 对蛋白序列构建索引
makeblastdb -in data/protein.fa -dbtype prot -title protein.fa
# 进行序列比对
blastp -query ~/data/protein.fa -out Os.blast -db ~/data/protein.fa -outfmt 6 -evalue 1e-5 -num_threads 30
3.mcscanx进行共线性分析
注：两个文件放到同一个文件夹中 Os.blast & Os.gff，通过mcscanx找到同源基因
mcscanx ./kaks/Os
- mcscanx输出.tandem的文件需要按置表符分为2列
sed -i s'/,/\t/g' *.tandem
4.计算水稻全部ka/ks
- 输出写绝对路径，输出文件夹无需先创建，需要先整理CDS序列与蛋白序列格式
- 需要在ParaAT中执行 chmod 777 Epal2nal.pl
数据整理
touch proc;echo 30 > proc # 设置线程数
# 将序列多行变一行
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
数据获取
根据个人电脑性能需要耗费时间不同，blast比对 & ka/ks计算此两步最耗费时间；一般对于小基因组物种个人电脑需要耗费几十小时不等；今天的介绍到此结束
