叶绿体基因组的GetOrganelle组装、批量任务和结果合并

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)))
}
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,332评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,508评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,812评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,607评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,728评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,919评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,071评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,802评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,256评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,576评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,712评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,389评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,032评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,798评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,026评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,473评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,606评论 2 350

推荐阅读更多精彩内容