1 叶绿体的组装
1.1 GetOrganelle的安装
此处默认Anaconda已经安装,那么使用下面的命令即可安装GetOrganelle:
conda create -n getorganelle python=3.7 #创建名叫getorganelle的环境
conda install -n getorganelle -c bioconda getorganelle #安装getorangelle
1.2 GetOrganelle的配置
安装好后需要配置参考基因组,GetOrganelle作者推荐以下命令:
get_organelle_config.py --add embplant_pt,embplant_mt
然而该命令下载很慢,可以考虑到本地下载,手动导入,具体参考命令如下:
wget https://github.com/Kinggerm/GetOrganelleDB/releases/download/0.0.1/v0.0.1.tar.gz
tar -zxvf v0.0.1.tar.gz
get_organelle_config.py -a embplant_pt,embplant_mt --use-local ./0.0.1
该命令参考自:零基础教程 | 叶绿体基因组组装 - GetOrganelle
之后即可开始执行组装,但是可能会提示:
ERROR: Bowtie2 is not available!
实际上Bowtie2已经安装,但是似乎没有连接上,参考零基础教程 | 叶绿体基因组组装 - GetOrganelle的意见,直接更新bowtie2到最新版即可,实测也的确有效。
conda update bowtie2
1.3 使用GetOrganelle组装叶绿体
准备好自己的数据,一般为双末端的illumina数据,然后用以下命令运行:
get_organelle_from_reads.py -1 forward.fq -2 reverse.fq -o plastome_output
-R 15 -k 21,45,65,85,105 -F embplant_pt # 该代码来自getorganelle作者官网
-1 -2 输入双向测序数据,解压或者压缩文件均可
-o 指定输入文件的路径
-R 似乎是sample的循环迭代次数
-k K-mer参数,不可超过reads的最大读长
-F 指定参考基因组,一般植物基因组为embplant_pt
开始运行后,可能会有一个SPAdes不会正常调用的报错:
......
2021-09-19 00:51:50,106 - ERROR: Error with running SPAdes: == Error == system call for: "['/home/zz/.conda/envs/getorganelle/share/spades-3.13.0-0/bin/spades-core', '/home/zz/genome/chloroplast/GetOrganelle-master/dinganensis/seed/embplant_pt.initial.fq.spades/K45/configs/config.info']" finished abnormally, err code: -11
......
该报错导致程序中断无法运行,经过google检索,发现似乎是个挺普遍的问题,参考Spades finished abnormally, err code: -11 · Issue #83 · Kinggerm/GetOrganelle · GitHub里面的建议,本地安装SPAdes的最新版本即可,因为 conda install spades -c bioconda 默认安装的不是最新版,奇怪的是bioconda源的确是有最新版的包的。
考虑到bioconda源中有最新版本的包,尝试直接在线安装指定版本:
conda install spades=3.15.3
成功安装,getorganelle终于成功运行。
1.4 GetOrganelle组装叶绿体的参数优化
只要重测序的数据有一定的覆盖度(>2X),官方的命令一般即可组装成环,但是也常常出现不成环的情况,然后就需要对参数优化或者后续深入分析。
参数优化分为两个方向,一是增加reads的使用量和迭代运算量,二是使用自己的参考基因组:
# 增加数据使用量和迭代细节
get_organelle_from_reads.py -1 forward.fq -2 reverse.fq -o plastome_output
-R 30 -k 21,35,45,65,85,105,121,135,145 --reduce-reads-for-coverage 10000000
-F embplant_nr --overwrite
- R -k 提高迭代循环和k-mer粒度,也可设置更高的数值
--reduce-reads-for-coverage 10000000 个reads或者inf使用全部的reads
这些设置可能使组装成环,但是运算量会成倍增加
## 使用自己的参考基因组
get_organelle_from_reads.py -1 forward.fq -2 reverse.fq -o plastome_output
-R 15 -k 21,45,65,85,105 -F embplant_pt -s personal-reference.fasta
-s 指定自己的参考基因组,比如同一个物种、属的参考基因组,可以同时指定多个,放到同一个fasta文件即可
如果通过调整GetOrganelle参数仍然不能成环,可以使用Bandage来手动选择、删减,组合不同的config生成成环的组转结果,大致界面如下。具体可以参考此视频,或者Bilibili上的搬用视频,后续我们也可能考虑写个教程。
2. GetOrganelle的批量组装
如果批量测序了一批数据,逐个组装似乎有些麻烦,因此考虑批量组装,利用Linux的bash命令如下:
while read a;
do get_organelle_from_reads.py -1 "$a"R1.fastq -2 "$a"R2.fastq -o "$a"-output
-R 15 -k 21,45,65,85,105 -F embplant_pt;
done < list.txt
# 序列文件的前缀放到list.txt即可,会将list文件中的每一行作为a值,替换到do后面的命令中
3. GetOrganelle多个组装结果的合并
GetOrganelle的成环结果文件中通常同时包括两个构型的结果,但两种构型的顺序似乎是随机的。那么多个结果的合并就有些麻烦,需要比对后根据比对结果重选挑选另一个构型的结果。一般来说,不同构型结果的混合比对,速度相当缓慢。那么,使用以下的半成品小脚本可以用来合并多个GetOrganelle的结果,保证他们的构型一直,且同时可以检查不成环的结果。
# Author: Zhang Zhen
# Email: ficuszz@163.com
# Website: zhangzhen.plus
# Description: 这是一个半成品的脚本,需要自行调整一些诸如工作路径、文件输入路径、部分阈值之类的
# 参数。所以这里假设需要了解初步的R语言知识
library(ape)
library(strataG)
setwd("g:/Ubuntu_1804.2019.522.0_x64/rootfs/home/zz/byk/138-202209/") # 指定测序文件夹,该文件夹中应该只有多个样品测序的结果文件夹
file_path <- list.files()
sample_name <- sapply(file_path, function(x) substr(x, 1, nchar(x)-14)) # 此处的14为需要修剪掉的字符,具体随自己测序文件的命名规则而变
sample_number <- length(list.files())
# 将第一个样品的其中一条组装结果作为基准,将其余样品的同构型结果收集起来(请确保第一个样品组装成环)
j <- 1 # 如果需要选择另一个构型,请改为2
if(length(list.files(file_path[1], pattern = "fasta", full.names = T )) == 2 ) {
temp0 <- read.FASTA(list.files(file_path[1], pattern = "fasta", full.names = T )[j])
names(temp0) <- sample_name[1]
write.FASTA(temp0, paste("c:/Users/qbycs/Desktop/byk/",sample_name[1], ".fasta", sep = ""))
} else {print("请确保第一个样品组装成环") }
#将其余样品的同构型组装结果筛选出来
checklist <- c()
for(i in 1:sample_number){
if(length(list.files(file_path[i], pattern = "fasta", full.names = T )) == 2){ # 检查是否成环组装
temp <- read.FASTA(list.files(file_path[i], pattern = "fasta", full.names = T)[1])
# 根据比对后遗传距离来判断构型
if(dist.dna(mafft(c(temp0,temp))) < 0.01){ # 认为和基部遗传距离相差超过0.01为不同构型,然后自动选择另一条,可根据类群远近调整数值
names(temp) <- sample_name[i] # 修改fasta的label
write.FASTA(temp, paste("c:/Users/qbycs/Desktop/byk/",sample_name[i], ".fasta", sep = ""))
} else {
temp <- read.FASTA(list.files(file_path[i], pattern = "fasta", full.names = T )[2])
names(temp) <- sample_name[i] # 修改fasta的label
write.FASTA(temp, paste("c:/Users/qbycs/Desktop/byk/",sample_name[i], ".fasta", sep = ""))
}
} else {
checklist <- c(checklist, sample_name[i])
next;}
}
# 输出问题组装样品名称
if(length(checklist) > 0){
print(paste("请检查以下样品号的组装结果:", as.character(checklist)))
}