继上一次浅层基因组数据分析笔记 - Mushroom Monkey (jungleblack007.github.io) 之后,我们发现Orthofinder很多情况下能找到的单拷贝直系同源基因非常少,而且文章里面也没有使用这种方法。
这次,我要构建蘑菇目的系统发育树!
下载蛋白序列 这个没啥好说的,JGI Genome Portal - Home (doe.gov) 可以下载一部分基因组数据,GenBank Overview (nih.gov) 中也可以下载一部分。直接下载蛋白的数据会更方便一些,有些可能没有注释,我们下载基因组的数据然后自己注释一下也行。
下载新的BUSCO 我比较菜(lan),选择conda安装。这个更新很慢,可以多等一会。
1 2 conda create -n busco -c conda-forge -c bioconda busco=5.2.2 conda activate busco
在Index of /v4/data/lineages/ (ezlab.org) 下载一个OrthoDB的数据库作为预定义的参考单拷贝直系同源基因。因为我要做蘑菇目的,所以下载了agaricales_odb10.2020-08-05.tar.gz ,做其他的话可以选择适合自己实验设计的数据库。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 # 下载数据库 cd /database/BUSCO/ wget https://busco-data.ezlab.org/v4/data/lineages/agaricales_odb10.2020-08-05.tar.gz tar -zxvf agaricales_odb10.2020-08-05.tar.gz # 测试命令 cd /mnt/d/Agaricales3/ busco -i Agaricus_bisporus.fasta -o Agaricus_bisporus -m prot -c 18 -l ~/database/BUSCO/agaricales_odb10 # busco常用参数介绍: # -i:输入文件 # -o:输出文件夹 # -m:模式选择: # - geno or genome, for genome assemblies (DNA) # - tran or transcriptome, for transcriptome assemblies (DNA) # - prot or proteins, for annotated gene sets (protein) # -c:cpu数 # -l:之前参考数据集存放的文件夹
运行完之后,在这个文件夹里/mnt/d/Agaricales3/Agaricus_bisporus/run_agaricales_odb10/busco_sequences/single_copy_busco_sequences就会显示出找到的单拷贝直系同源基因了,当然我们这里都翻译成了蛋白了,这也就是我们要建树用的序列。
我们做一个循环,把所有的蛋白序列都运行一遍。
1 2 3 4 5 6 7 8 9 10 11 # 首先把所有的文件名放到一个list中 ls | grep ".fasta" > list.txt sed -i 's/.fasta//' list.txt # 统计list.txt的行数是否正确 wc -l list.txt # 挨个进行busco吧 list=$(cat list.txt) for i in $list;do busco -i ${i}.fasta -o ${i} -m prot -c 18 -l ~/database/BUSCO/agaricales_odb10 done # 跑上我就回去睡觉了,结果中断了,我才知道这个软件需要联网
下面我们就要想一个办法,把3870个片段中在所有*/run_agaricales_odb10/busco_sequences/single_copy_busco_sequences文件夹里都存在的给找出来,然后以物种名命名。
1 2 3 4 5 6 7 8 9 # 查找某条序列在single_copy_busco_sequences的个数 find ./*/run_agaricales_odb10/busco_sequences/single_copy_busco_sequences -name "0at5338.faa"| wc -l list2=$(cat list2.txt) for i in $list2;do if find ./*/run_agaricales_odb10/busco_sequences/single_copy_busco_sequences -name '${i}'| wc -l == cat ${i} > list3.txt; done