04叶绿体流程1-组装

一.前期准备

数据位置
/share/nas1/seqloader/yelvti/GP-xxxx_yelvti/2.clean_data

#先看污染比对结果!!!!有问题直接先反馈!!!!
1.污染比对
perl /share/nas6/zhouxy/functional_modules/pollution_nt_blast/pollution_nt_blast_pip_v3.pl -fqdir 路径/2.clean_data -od 输出目录

2.得到ass文件,改成拉丁名,下划线
python3 /share/nas1/yuj/script/chloroplast/get_ass_cfg.py -i 路径(邮件分发) #会自动补上2.clean_data

二.组装

数据库位置

/share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/assembly/dat/cpdb/plastid.genome.sequence.fasta

nohup perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/assembly/cp.assembly.pip.pl -i ass.cfg & # 组装核心流程

当bowtie超过指定大小立马组装

以下脚本为:
当bowtie超过50M或者3_Spades_Assembly已经创建,立马切断流程,进行组装
这样可以节约时间

#!/bin/bash

while true; do
    if  [ -d "analysis/assembly/"*"/pseudo/"*"/3_Spades_Assembly" ]; then
        echo "文件夹 '3_Spades_Assembly' 已创建"
        break  # 跳出循环
    fi
    if [ -f "analysis/assembly/"*"/pseudo/"*"/2_Bowtie_Mapping/map_pair_hits.1.fq" ]; then
        if ((`stat -c %s "analysis/assembly/"*"/pseudo/"*"/2_Bowtie_Mapping/map_pair_hits.1.fq"` > 50 * 1024 * 1024)); then
            echo "文件 'map_pair_hits.1.fq' 大小超过50M"
            break  # 跳出循环
        fi
    fi
    sleep 30  # 等待一段时间后再次进行检查
done

echo "条件满足，可以继续执行后续命令"

# 获取之前的进程ID
previous_pid=$(pgrep -f "/share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/assembly/cp.assembly.pip.pl")

# 终止之前的进程
if [ -n "$previous_pid" ]; then
    kill "$previous_pid"
fi

cd "analysis/assembly/"*"/pseudo/"*"/2_Bowtie_Mapping"
head -n 500000 map_pair_hits.1.fq > cut.1.fq
head -n 500000 map_pair_hits.2.fq > cut.2.fq
unicycler.py -1 cut.1.fq -2 cut.2.fq -o uni

三.组装调整

3.1批量运行

# 批量展示流程结果
for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;ll  analysis/assembly/$i/pseudo/$i/Final_Assembly/$i"_FULLCP.fsa"  ||  ll analysis/assembly/$i/pseudo/$i/Final_Assembly/$i"_afin_iter2.fa" || ll analysis/assembly/$i/pseudo/$i/4_Afin_Assembly/filtered_spades_contigs.fsa;done

## spa 批量移动3_Spades_Assembly的gfa
mkdir spa && for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;cp analysis/assembly/$i/pseudo/$i/3_Spades_Assembly/spades_iter1/assembly_graph.fastg spa/$i"_assembly_graph.fastg";done

## org 批量移动3_Spades_Assembly的gfa
mkdir org && for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;cp analysis/assembly/$i/pseudo/$i/3_Spades_Assembly/spades_iter1/extended_K121.assembly_graph.fastg org/$i"_K121.assembly.fastg" || cp analysis/assembly/$i/pseudo/$i/3_Spades_Assembly/spades_iter1/extended_spades/assembly_graph.fastg org/$i"_assembly_graph.fastg";done

# 批量查看bowtie的reads数
for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;cd analysis/assembly/$i/pseudo/$i/2_Bowtie_Mapping && showline- && cd -;done

# 批量2_Bowtie_Mapping cut
for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;cd analysis/assembly/$i/pseudo/$i/2_Bowtie_Mapping && head -n 300000 map_pair_hits.1.fq > cut.1.fq && head -n 300000 map_pair_hits.2.fq > cut.2.fq && cd -;done

# 批量unicycler
for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;cd analysis/assembly/$i/pseudo/$i/2_Bowtie_Mapping && echo `pwd`; unicycler.py -1 *.1.fq -2 *.2.fq -o uni &  cd -;done

# 批量移动unicycler的gfa
mkdir uni && for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;cp analysis/assembly/$i/pseudo/$i/2_Bowtie_Mapping/uni/002_overlaps_removed.gfa uni/$i"_002_overlaps_removed.gfa";done

# 批量移动finish
for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;cp analysis/assembly/$i/finish/$i"_FULLCP.fsa" fsa/$i"_FULLCP.fsa";done

# 批量把调整好的结果放回finish
for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;mkdir analysis/assembly/$i/finish -p ;cp ass/${i}_FULLCP.fsa analysis/assembly/$i/finish;done

1.批量bowtie
for i in `ls analysis/assembly/`;do echo $i; cd analysis/assembly/$i/pseudo/$i/ && cp_bowtie.align.pl -i  $i.fa*  -1 1_Trimmed_Reads/*P1.fq -2 1_Trimmed_Reads/*P2.fq  -o bowtie  & done

2.全部check
for i in $( ls /share/nas1/yuj/project/202201/GP-20211206-3816/analysis/assembly) ;do echo $i; cd /share/nas1/yuj/project/202201/GP-20211206-3816/analysis/assembly/$i/finish && check -1 ../pseudo/$i/4_Afin_Assembly/bowtie/*1.fq -2 ../pseudo/$i/4_Afin_Assembly/bowtie/*2.fq  -k 125 -m 3 -g *.fsa -o check;done

3.批量check
for i in {Mm_J3,Mm_L1,Mm_L2,Mm_L3,Mm_S1,Mm_S2,Mm_S3};do echo $i;cd /share/nas1/yuj/project/202201/GP-20211206-3816/analysis/assembly/$i/pseudo/$i/4_Afin_Assembly/ && check -1 bowtie/*1.fq -2 bowtie/*2.fq  -k 125 -m 3 -g fil*.fsa -o check ;done

3.2 100次迭代

cd analysis/assembly/*/pseudo/*/

1.组装
1)使用unicycler
mkdir -p bow1/uni/ # 创建bow1文件夹,里面有/uni/assembly.fasta
for i in {1..100} ;do echo $((i+1)) && cp_bowtie.align.pl -i bow$i/uni/assembly.fasta -1 1_Trimmed_Reads/*trimmed_P1.fq -2 1_Trimmed_Reads/*trimmed_P2.fq -o bow$((i+1)) && unicycler.py -1 bow$((i+1))/*1.fq -2 bow$((i+1))/*2.fq -o bow$((i+1))/uni --kmers 21,55,127 > log2 ;done & # 指定kmer

for i in {1..100} ;do echo $((i+1)) && cp_bowtie.align.pl -i bow$i/uni/assembly.fasta -1 1_Trimmed_Reads/*trimmed_P1.fq -2 1_Trimmed_Reads/*trimmed_P2.fq -o bow$((i+1)) && cd bow$((i+1)) && uni-  > log2 && cd - ;done & # 不指定kmer

2)使用spades
mkdir -p bow1/spades/
for i in {1..100} ;do echo $((i+1)) && cp_bowtie.align.pl -i bow$i/spades/scaffolds.fasta -1 1_Trimmed_Reads/*trimmed_P1.fq -2 1_Trimmed_Reads/*trimmed_P2.fq -o bow$((i+1)) && cd bow$((i+1)) && spa- && cd - ;done > 100tmp.log &

----监控文件夹数量
sh /share/nas1/yuj/script/mitochondrion/assembly/monitor.sh &

2.运行完毕
mkdir uni

for i in `ls analysis/assembly/`;do echo $i---------------------------------------------------------------------------------------------------------------------------------------------------------;cp analysis/assembly/$i/pseudo/$i/bow101/uni/002_overlaps_removed.gfa uni/$i"_002_overlaps_removed.gfa";done

3.3 使用延伸程序


序列延伸：
  长序列：perl /share/nas6/xul/program/mt2/assembly/src/extend.auto.pl -i1 target.fa(组装序列) -i2 all.fq/fa(原始数据r1、r2各运行一次) 
  短序列：perl /share/nas6/xul/program/mt2/assembly/src/extend.auto.pl -i1 target.fa -i2 all.fq/fa -sr  #优先运行这个

延伸后序列替换：
  perl /share/nas6/xul/program/mt2/assembly/src/extend.pl  -t target.fa(组装序列) -f _extend.map.fa -c _out.coords   -id  999(auto程序展示结果中的匹配reads)

3.4 手动步骤

3.4.1 bowtie

$ cp_bowtie.align.pl -i allseq.fa  -1 ../1_Trimmed_Reads/*P1.fq -2 ../1_Trimmed_Reads/*P2.fq  -o bowtie &
$ cp_bowtie.align.pl -i spa.fsa  -1 ../1_Trimmed_Reads/*P1.fq -2 ../1_Trimmed_Reads/*P2.fq  -o bowtie &
$ cp_bowtie.align.pl -i test.fa  -1 ../1_Trimmed_Reads/*P1.fq -2 ../1_Trimmed_Reads/*P2.fq  -o bowtie &
$ cp_bowtie.align.pl -i filtered_spades_contigs.fsa  -1 ../1_Trimmed_Reads/*P1.fq -2 ../1_Trimmed_Reads/*P2.fq  -o bowtie &

$ showline-
$ head -n 1000 map_pair_hits.1.fq > cut.1.fq && head -n 1000 map_pair_hits.2.fq > cut.2.fq
$ unicycler.py -1 cut.1.fq -2 cut.2.fq -o unicycler

3.4.2 check

$ check -1 bowtie/*1.fq -2 bowtie/*2.fq  -k 125 -m 3 -g test.fa -o check &
$ check -1 bowtie/*1.fq -2 bowtie/*2.fq  -k 125 -m 3 -g filtered_spades_contigs.fsa -o check &
$ check -1 ../2_Bowtie_Mapping/*1.fq -2 ../2_Bowtie_Mapping/*2.fq  -k 125 -m 3 -g test.fa -o check &
$ check -1 ../2_Bowtie_Mapping/*1.fq -2 ../2_Bowtie_Mapping/*2.fq  -k 125 -m 3 -g filtered_spades_contigs.fsa -o check &

3.4.3 correct

$ correct -i1 check/K125.scaffold.fa  -i2 check/tmp/K125.dump  -s 125 -p K125 -o correct

3.4.4 tie

4.tie
# 对原来的test.fa末尾tie
$ tie -i check/K125.scaffold.fa -d check/tmp/K125.dump  -k 125  -tp 146000 -p K125 -o tie
$ tie -i test.fa -d check/tmp/K125.dump  -k 125  -tp 15400 -p K125 -o tie
# result —— tie/K125.extend.fa

5.ssc
$ cp_ssc.pl -i test.fa  -o re_ssc.fa

四.组装软件

4.1unicycler

unicycler.py -1 *1.fq -2 *2.fq -o uni

4.2spades

spades.py --pe1-1 *1.fq --pe1-2 *2.fq --careful  -o spades               -k 35,45,55,75,97,107

spades.py --pe1-1 *1.fq.gz --pe1-2 *2.fq.gz -t 200 -k 97,107,117,127 -m 600 --careful --phred-offset 33 -o spades
# 全部数据组装  kmer可以设置的长一些

4.3Getorganelle

# 路径 /share/nas1/yuj/software/miniconda3/envs/getorganelle/bin/get_organelle_from_reads.py

# 其他参数  
-R 15 -k 21,45,65,85,105,115 -s 参考
1.常规使用
/share/nas1/yuj/software/miniconda3/envs/getorganelle/bin/get_organelle_from_reads.py -1 *1.fq -2 *2.fq -o org -F embplant_pt
get_organelle_from_reads.py -1 *1.fq -2 *2.fq -o org -F embplant_pt

# 更快的方法	
/share/nas1/yuj/software/miniconda3/envs/getorganelle/bin/get_organelle_from_reads.py -1 *1.fq -2 *2.fq -o org -F embplant_pt --fast                   -w 0.68(可不设置)

2.使用已有fastg组装
/share/nas1/yuj/software/miniconda3/envs/getorganelle/bin/get_organelle_from_reads.py -g assembly_graph.fastg -F embplant_pt -o org

4.4 串行运行

unicycler.py -1 *1.fq -2 *2.fq -o uni && get_organelle_from_reads.py -1 *1.fq -2 *2.fq -o org -F animal_mt && spades.py --pe1-1 *1.fq --pe1-2 *2.fq --careful  -o spades

五.下载参考(后接注释)

perl /share/nas6/xul/program/chloroplast/bin/cp_get_genbank_form_ncbi_with_ACCESSION.pl

叶绿体以前的记录

python3 /share/nas1/yuj/cp/get_ass_cfg.py -i ?????/2.clean_data #改物种名
nohup perl /share/nas6/pub/pipline/genome-assembly-seq/mitochondrial-genome-seq/v2.0/assembly/mt.assembly.pip.pl -i ass.cfg &
程序暂停  fg一直不好用  bg
cp_bowtie.align.pl -i test.fa  -1 ../1_Trimmed_Reads/*P1.fq -2 ../1_Trimmed_Reads/*P2.fq  -o bowtie &


----寻找组装用参考
cp   备份一下为test.fa,复制到ncbi找参考
perl /share/nas6/xul/program/chloroplast/bin/cp_get_genbank_form_ncbi_with_ACCESSION.pl 下载参考


一.重新延伸组装
1.先画图看看,大致了解
nuc xxxx.fasta  test.fa
mummerplot out.delta
return test.fa

2.生成bowtie文件夹(可选,建议生成)
cp_bowtie.align.pl -i test.fa  -1 ../1_Trimmed_Reads/xxxx.trimmed_P1.fq -2 ../1_Trimmed_Reads/xxxx.trimmed_P2.fq  -o bowtie &
耗时比较久,大约10多分钟
return bowtie/

3.开始检查,生成kmer库, 用bowtie/(优先建议),../1_Trimmed_Reads/(建议),../2_Bowtie_Mapping/,均可
check -1 ../1_Trimmed_Reads/xxxx.trimmed_P1.fq -2 ../1_Trimmed_Reads/xxxx.trimmed_P2.fq  -k 125 -m 3 -g test.fa -o check & #m为最小深度,选3 5 9 19等
return check/tmp/K125.dump

4.加载kmer库延伸
从ref截一段去kmer库找到kmer0TTTT
cp_extend_with_dump3.pl -i check/tmp/K125.dump  -s TTTT可以加入手动选择参数 -c
vim _extend3.fa删除125个
return _extend3.fa

5.画图比对删改调整
nucmer .../fasta/xxxx.fasta  _extend3.fa
mummerplot out.delta 
修改至成环
return cir_extend.fa

6.调整起点,用ir查看结构
ir cir_extend.fa 采用lsc做开头
return cir_extend.fa

7.检查结果
check -1 ../1_Trimmed_Reads/xxxx.trimmed_P1.fq -2 ../1_Trimmed_Reads/xxxx.trimmed_P2.fq  -k 125 -m 3 -g cir_extend.fa  -o check2 &
check -d check/tmp/K125.dump  -k 125 -m 3 -g cir_extend.fa -o check2 &第2步生成了kmer库,也可以直接用它
检查scaffold和覆盖度
return cir_extend.fa

8.结束,创建finish文件夹
xxxx_FULLCP.fsa 




二.使用已有结果修补
1.以物种9-Phalaenopsis_stobartiana为例
cp Phalaenopsis_stobartiana_afin_iter2.fa test.fa
nucmer ../../../../Phalaenopsis_stobartiana/pseudo/ref/fasta/MW531729.1.fasta  test.fa
mummerplot out.delta 
对序列进行操作,切一个环出来
return test.fa

2.生成bowtie文件夹(可选,建议)
cp_bowtie.align.pl -i test.fa  -1 ../1_Trimmed_Reads/*P1.fq -2 ../1_Trimmed_Reads/*P2.fq  -o bowtie &
耗时比较久,大约10多分钟
return bowtie/

3.利用bowtie里的俩文件进行检查,或者利用../2_Bowtie_Mapping/里的来check(注意:../2_Bowtie_Mapping/这里面是18年以前的物种,用的话有可能漏掉reads)
check -1 bowtie/*1.fq -2 bowtie/*2.fq  -k 125 -m 3 -g test.fa -o check &
找check里的scaffold,看看有没有N,本例中间末尾均有n
查看k125覆盖情况即kmer深度,若中间有N,则需要查看,末尾有,不看也行
大于70一般没问题,出现突降突升则有问题,寻找0-7的位置,/\t0\.[0-7]
覆盖深度低的地方全部替换为N,可以多换点,得到新的scaffold2
return check/K125.scaffold.fa

4.利用correct程序修补中间N
correct -i1 check/K125.scaffold.fa  -i2 check/tmp/K125.dump  -s 125 -p K125 -o correct
生成了correct,得到新补的gapfiller,中间没有gap了
找correct里的K125.iteration1.gapfilled.log查看
fastalength K125.gapfiller_finish.fa
return correct/K125.gapfiller_finish.fa

5.用tie延伸(针对末尾有N的情况)
tie -i check/K125.scaffold.fa -d check/tmp/K125.dump  -k 125  -tp 146000 -p K125 -o tie
进到tie文件夹,查看下
打开log文件找warn的位置,有可能出现分支,物种9没有分支
查看extend.fa的长度,有可能过长
return tie/extend.fa

6.tie延伸较长的解决办法
方法一:进行序列的反向互补,反向互补后的序列再tie延伸
方法二:
1.打开NCBI,粘贴进去,点击等待结果
2.比对完,下载参考序列放ref文件夹里,顺便挪到pseudo下
画参考图,找到最后面的正确位置  145895
找到145896处,后面全去掉,包括145896的T
3.用tie延伸一下,选择140000开始,本例延伸后依旧过长
得到新的extend与kmer库check一下得到check3文件夹
然后用最新的scaffold,此处为最终结果
return check3/K125.scaffold.fa

7.创建finish文件夹,改名,注意格式
mkdir ../../../finish
cp check3/K125.scaffold.fa  ../../../finish/Phalaenopsis_stobartiana_FULLMT.fsa

一.前期准备

二.组装

三.组装调整

3.1批量运行

3.2 100次迭代

3.3 使用延伸程序

3.4 手动步骤

3.4.1 bowtie

3.4.2 check

3.4.3 correct

3.4.4 tie

四.组装软件

4.1unicycler

4.2spades

4.3Getorganelle

4.4 串行运行

五.下载参考(后接注释)

叶绿体 以前的记录

叶绿体以前的记录