抱歉,您的浏览器无法访问本站

本页面需要浏览器支持(启用)JavaScript


了解详情 >

    继上一次浅层基因组数据分析笔记 - 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

评论