Getorganelle组装植物叶绿体基因组建树

一、软件安装

conda install -c bioconda getorganelle
#mamba可以大幅提高下载速度，界面也很帅！
mamba install -c bioconda getorganelle

安装成功后还需要下载数据库：

get_organelle_config.py --add embplant_pt,embplant_mt
#如果要组装ITS序列再加上embplant_nr数据库

安装mamba参考以下链接即可：
生信软件安装神器 Mamba [conda加速器~] - 简书 (jianshu.com)

mamba 安装biopython

二、组装

get_organelle_from_reads.py -1 forward.fq.gz -2 reverse.fq.gz -o plastome_output -R 15 -k 21,45,65,85,105 -F embplant_pt -o goal_dir
#直接输入测序文件，软件会自己判断所需数据量（一般2G的数据即可）
#-R是指循环数，理论越多越好，-k一般不用条
#-o会将结果输出到指定目录，我是一个个体一个目录

如果结果出现了complete,则说明组装成功。

组装结果

如果出现了scaffold，contig则说明组装没有成环。可以试着调节一下-R-和-w参数。

image.png

get_organelle_from_reads.py -1 forward.fq.gz -2 reverse.fq.gz -o plastome_output -R 30 -w 95 -k 21,45,65,85,105 -F embplant_pt
#将-R设为30
#将-w设为95，默认是105，官方建议从105-65范围递减尝试

我将参数调整以后，之前组装不好的就可以了。
除此之外还有几个参数也可以调节，可参考官方文档：
FAQ · Kinggerm/GetOrganelle Wiki (github.com)

参数调节

三、叶绿体基因组结果检查及注释

组装成环后，每个个体会得到两种叶绿体基因组序列，对应不一样的构型。

叶绿体基因组组装结果
这两种都是正确的，这是因为LSC和SSC大小单拷贝区顺序不一样，官方文档说明的很清楚，二者都正确且在自然界存在，我们只需挑选一个即可。
叶绿体基因组结构
注释可以用GeSeq、Chloroplot、IRscope等在线软件，参考：
七. 植物叶绿体基因组--组装，注释和比较作图（GetOrganelle，GeSeq，Chloroplot，IRscope） - 简书 (jianshu.com)

四、统一叶绿体基因组构型来建树

1. 寻找参考
   此步骤我们需要从两个序列中挑选一个，为了严谨，可以去NCBI中下载一个别人已经发表过的此物种或者近缘物种的叶绿体基因组当参考。跟参考一样的则保留。

2. 合并所有个体叶绿体序列
   如果是按照我之前的步骤将每个个体叶绿体基因组输出到对应个体的文件夹中,

   
   结果
   

   这里我写了一个脚本修改一下输入输出就可以将该路径中所有目录下叶绿体序列合并到一个fasta文件中。

import os
from Bio import SeqIO
file_path = "plastome_assemble_results/"#这里填写组装结果目录
plastid_fasta_list = []
ref_fasta = SeqIO.read("assembly_plastome/ref.fasta","fasta")#填写NCBI下载参考序列，要fasta格式！！！
plastid_fasta_list.append(ref_fasta)

for current,dirs,files in os.walk(file_path):
    for file in files:
        the_path = current+'/'+file
        if file == 'embplant_pt.K105.complete.graph1.1.path_sequence.fasta':
            indiv_name = the_path.split('/')[-2]
            fasta = SeqIO.read(the_path,"fasta")
            fasta.id = indiv_name+'.1'
            fasta.name = ""
            fasta.description = ""
            plastid_fasta_list.append(fasta)
        elif file == 'embplant_pt.K105.complete.graph1.2.path_sequence.fasta':
            indiv_name = the_path.split('/')[-2]
            fasta = SeqIO.read(the_path,"fasta")
            fasta.id = indiv_name+'.2'
            fasta.name = ""
            fasta.description = ""
            plastid_fasta_list.append(fasta)
        else:
            pass
SeqIO.write(plastid_fasta_list,"assembly_plastome/all.plastid.fasta","fasta")#填写输出的fasta结果文件

脚本合并结果每条fasta序列以及两种结构的命名会以.1和.2为后缀：

合并结果

3. mafft比对
   用mafft比对

mafft --thread 20 --auto all.plastid.fasta > all.plastid.mafft.fasta
#--thread线程数目
#--auto比对模式，auto即可

4. BioEidt计算一致性
   软件下载安装教程：
   Bioedit 使用_韩建刚（CAAS-UCD）的博客-CSDN博客
   图标
   用BioEdit打开刚刚的序列比对文件，然后计算一致性矩阵，注意保存为csv格式更方便。
   一致性计算
   csv

结果用Excel打开，我们只关心第一列即与参考比对的结果即可，两种结构中留一致性高的，也就是与参考结构一致的结果。

一致性结果

IQtree建树

iqtree -s all.fasta -m TEST -nt AUTO
iqtree -s all.fasta -m best+ASC -nt AUTO

总结

后续考虑写一个脚本，将合并后的fasta直接计算序列一致性，然后去掉一致性差的序列，一步到位输出目标叶绿体基因组序列，可以拿来建树。