05线粒体流程1-组装

#线粒体 #组装 #SOP

一. 3代 —— Nanopore/PacBio(Hifi数据)

1.0 使用延伸程序

ass=sample.fa && file=pass.fq* && nu=1000000 && minimap2 -t 20 -ax map-ont $ass $file > extend.sam && samtools view -F 4 extend.sam  | perl -lane 'print unless($F[9] eq "*")' |  perl -ane 'print if(/^@/);if(/NM:i:(\d+)/){$n=$1;$l=0;$l+=$1 while $F[5]=~ /(\d+)[M]/g;if($l > 5000){print}}'|sort -k 4 -n| tail -n $nu > tmp.sam && cat tmp.sam  |cut -f 10 |perl -lane 'print ">",++$i;print $F[0]'  > map_gene.fa  && nucmer $ass map_gene.fa -p map_gene && show-coords  map_gene.delta  > map_gene.coords && mummerplot  map_gene.delta & # 包含下面1.1和1.2

grep -A 1 ">201$" map*.fa > test.fa

perl /share/nas6/xul/program/mt2/assembly/src/extend.auto.pl -i1 test.fa -i2 map*.fa -polish -re

perl /share/nas6/xul/program/mt2/assembly/src/extend.auto.pl -i1 test.fa -i2 map_gene.fa -tail 1000000 -n 0 
# -tail：显示的reads条数，默认300。-n为0时可以查看共线性（mummerplot _out.delta） 

perl /share/nas6/xul/program/mt2/assembly/src/extend.pl -t test.fa -f _extend.map.fa -c _out.coords -i -id 10086 
# 使用10086号reads进行延伸并替换test.fa文件

1.1 minimap2 —— 映射

0.PacBio(Hifi数据)需要用bedtools-bamtofastq转换格式 
bamToFastq -i *.bam -fq unmapped.fastq # PacBio

1.
minimap2 -ax map-ont ref.fa 测序数据.fq.gz > aln.sam  # Nanopore/ont
minimap2 -ax map-pb  ref.fa unmapped.fastq > aln.sam  # PacBio

2.使用igv查看一下覆盖深度
samtools view -h -F 4 aln.sam | samtools sort > 01_map.sort.bam && samtools index 01_map.sort.bam
igv 加载基因组和排序后的bam/sam

2.使用samtools统计覆盖深度
samtools depth -a 01_map.sort.bam > depth.txt 
# perl /share/nas6/xul/program/chloroplast/assembly/draw_line_depth.pl depth.txt && convert test.svg depth.png
perl /share/nas1/yuj/program/chloroplast/assembly/draw_line_depth.pl depth.txt depth.svg && convert depth.svg depth.png

samtools depth  01_map.sort.bam |awk '{print "mt1",$2,$2,$3}' > sample.depth.txt
python3 /share/nas1/yuj/script/chloroplast/annotation/draw_line_depth_v2.py sample.depth.txt line_depth.png

1）使用minimap2将3代数据回比基因组
2）使用samtools提取比对到基因组上的数据
3）对sam文件进行排序
4）使用samtools统计覆盖情况
5）可视化覆盖度

1.2 samtools —— 筛选

samtools view -F 4 aln.sam  | perl -lane 'print unless($F[9] eq "*")' |  perl -ane 'print if(/^@/);if(/NM:i:(\d+)/){$n=$1;$l=0;$l+=$1 while $F[5]=~ /(\d+)[M]/g;if($l > 5000){print}}'|sort -k 4 -n> tmp.sam && cat tmp.sam  |cut -f 10 |perl -lane 'print ">",++$i;print $F[0]'    > 5000map_gene.fa

samtools view -F 4 aln.sam | perl -lane 'print unless($F[9] eq "*")' | perl -ane 'print if(/^@/);if(/NM:i:(\d+)/){$n=$1;$l=0;$l+=$1 while $F[5]=~ /(\d+)[M]/g;if($l > 5000 && $l < 20000){print}}' | sort -k 4 -n > tmp.sam && cat tmp.sam | cut -f 10 | perl -lane 'print ">",++$i;print $F[0]' > 5000_20000map_gene.fa

# if($l > 1000) 筛选大于1000bp / if($l > 5000) 筛选大于5000bp
# 生成 map_gene.fa

1.3 awk —— 查看筛选后的平均长度

# 平均长度
$ ir.py -i 5000map_gene.fa -l | awk 'BEGIN {n=0;count=0} {count=count+$1;n=n+1} END{print "count/n is ",count,n;print "adv is",count/n;}'
$ fl map_gene.fa | awk 'BEGIN {n=0;count=0} {count=count+$1;n=n+1} END{print "count/n is ",count,n;print "adv is",count/n;}'

# 最大长度
$ ir.py -i 5000map_gene.fa -l | awk 'BEGIN {max = 0} {if ($1+0 > max+0) max=$1} END {print "Max=", max,"bp"}'
$ fl map_gene.fa > len.txt && awk  'BEGIN {max = 0} {if ($1+0 > max+0) max=$1} END {print "Max=", max,"bp"}' len.txt

# 最小长度
$ ir.py -i 5000map_gene.fa -l | awk 'BEGIN {min = 1000000} {if ($1+0 < min+0) min=$1} END {print "Min=", min,"bp"}'
$ fl map_gene.fa > len.txt && awk 'BEGIN {min = 1000000} {if ($1+0 < min+0) min=$1} END {print "Min=", min,"bp"}' len.txt

1.4 canu —— 校正三代数据

nohup /share/nas6/zhangxq/biosoft/canu-master/Linux-amd64/bin/canu -correct -p correct -d correct genomeSize=10k useGrid=false -nanopore-raw map_gene.fa & # Nanopore/ont

nohup /share/nas6/zhangxq/biosoft/canu-master/Linux-amd64/bin/canu -correct -p correct -d correct genomeSize=11k useGrid=false -pacbio-raw map_gene.fa & # pacbio

# genomeSize=10k 基因组预估大小
# 生成 correct/correct.correctedReads.fasta.gz

1.5 gz -d [file] —— 解压缩

gunzip –c correct.correctedReads.fasta.gz > correct.correctedReads.fasta #保留源文件
gunzip correct.correctedReads.fasta.gz #不保留源文件

# 压缩命令 gzip [file]
# 解压文件 gunzip [file]/gzip -d [file]

# 压缩保留源文件 gzip –c filename > filename.gz
# 解压保留源文件 gunzip –c filename.gz > filename

二. 2代数据

2.0 示例

# 如3816从mmj3开始
cp_bowtie.align.pl -i filtered_spades_contigs.fsa  -1 ../1_Trimmed_Reads/*P1.fq -2 ../1_Trimmed_Reads/*P2.fq  -o bowtie 
# 查看最开始的check文件夹
ir.py -i filtered_spades_contigs.fsa -o test.fa
check -1 bowtie/*1.fq -2 bowtie/*2.fq  -k 125 -m 3 -g test.fa -o check &
# 再次查看
tie -i test.fa -d check/tmp/K125.dump  -k 125  -tp 19100 -p K125 -o tie
nuc GQ903339.1.fa tie/*.fa && mum  
# 改图改名
nuc GQ903339.1.fa tie/*.fa && mum
check -1 bowtie/*1.fq -2 bowtie/*2.fq  -k 125 -m 3 -g tie/*.fa -o check2 &
# 查看2
vim check2/K125.coverage.txt  
# 没有问题
# check2对应的序列为结果,改名

2.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

2.2 组装

数据库位置
/share/nas6/pub/pipline/genome-assembly-seq/mitochondrial-genome-seq/v2.0/assembly/dat/mtGenome_db_20201120/allmtseq_20201120.fa

## 组装流程
nohup perl /share/nas6/pub/pipline/genome-assembly-seq/mitochondrial-genome-seq/v2.0/assembly/mt.assembly.pip.pl -i ass.cfg &
# 样本多出错的话,去掉nohup 去掉&

## 手动运行
perl /share/nas6/pub/pipline/genome-assembly-seq/mitochondrial-genome-seq/v2.0/assembly/mtDNA_assembly.pl -1 1_Trimmed_Reads/Neuroctenus_yunnanensis_Hsiao.trimmed_P1.fq -2 1_Trimmed_Reads/Neuroctenus_yunnanensis_Hsiao.trimmed_P2.fq -p Neuroctenus_yunnanensis_Hsiao -o analysis/assembly/Neuroctenus_yunnanensis_Hsiao

## 依次和参考看一下起点
for i in ref_ass/fasta/*;do echo $i; nuc $i *.fsa -l 1 && awk 'NR==4 {print $3}' out.delta;done

2.3 手动修改

2.3.1 串行

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 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

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 bowtie/*1.fq -2 bowtie/*2.fq  -k 125 -m 3 -g tie/K125.extend.fa -o check2 &
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.correct
correct -i1 check/K125.scaffold.fa  -i2 check/tmp/K125.dump  -s 125 -p K125 -o correct

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

2.3.2 使用延伸程序

序列延伸：
  长序列：perl /share/nas6/xul/program/mt2/assembly/src/extend.auto.pl -i1 target.fa(组装序列) -i2 all.fq(原始数据r1、r2各运行一次) 
  短序列：perl /share/nas6/xul/program/mt2/assembly/src/extend.auto.pl -i1 target.fa -i2 all.fq -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)

2.3.3 组装迭代100次

cd analysis/assembly/*/pseudo/*/

----组装
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 &
monitor.sh脚本内容如下所示

使用以下的Shell脚本来实现每时每刻只保留最新的5个bow文件，并且监控间隔为5分钟：

#!/bin/bash

# 设置监控间隔（单位：秒）
interval=300

while true; do
  # 获取当前目录下的所有bow文件夹，并按修改时间排序
  folders=$(ls -td bow[0-9]*/)

  # 计算当前文件夹数量
  count=$(echo "$folders" | wc -l)

  # 如果文件夹数量超过5个，删除多余的文件夹
  if [ $count -gt 5 ]; then
    delete_count=$((count - 5))
    delete_folders=$(echo "$folders" | tail -n $delete_count)
    rm -rf $delete_folders
    echo "Deleted $delete_count folders."
  fi

  sleep $interval
done

将上述脚本保存为一个名为monitor.sh的文件，并确保该文件具有可执行权限（可以使用chmod +x monitor.sh命令赋予执行权限）。
这个脚本使用一个无限循环来持续监控当前目录下的bow文件夹。它首先获取所有bow文件夹并按修改时间进行排序。然后，它检查当前文件夹数量是否超过5个，如果超过，则删除最旧的文件夹，以保持最新的5个文件夹。
脚本中的interval变量设置了监控间隔，单位是秒。在每次循环结束后，脚本会暂停指定的时间间隔，然后再次进行检查和删除操作。
你可以在终端中运行这个脚本，让它在后台持续监控并执行删除操作。执行命令./monitor.sh &即可将脚本放入后台运行。如果需要停止脚本的执行，可以使用killall monitor.sh命令来终止它。
请注意，这个脚本假设bow文件夹的命名是按照生成顺序进行编号的（例如：bow1、bow2、bow3...）。如果文件夹的命名方式不同，请根据实际情况修改脚本中的文件夹名称匹配规则。

2.3.4 通过节点和深度筛选gfa

cd 1_Trimmed_Reads/spades

blastn -query ref.fa -subject assembly_graph.fastg -outfmt 6 -out blast.log

获得node序号,写入id.txt

perl /share/nas6/xul/program/mt2/assembly/src/filter_node_from_fastg.pl -g assembly_graph.fastg  -i id.txt -d 99(深度)


将pe150 reads连接成更长的reads
/share/nas6/zhouxy/biosoft/pandaseq/current/bin/pandaseq  -f map_pair_hits.1.fq  -r map_pair_hits.2.fq  > merge.fa  2> log

2.3.5 批量运行

(1)全部修改起点
for i in *.fsa; do nuc ref.fa $i -l 1 && mt_move_pos.py -fa $i -s `awk 'NR==4 {print $3}' out.delta` ;done
---- 依次和参考看一下起点
for i in ref_ann/fasta/*;do echo $i; nuc $i M1_FULLMT.fsa -l 1 && awk 'NR==4 {print $3}' out.delta;done

(2)全部bowtie
for i in $( ls analysis/assembly) ;do echo $i; cd analysis/assembly/$i/pseudo/$i/4_Afin_Assembly && cp_bowtie.align.pl -i filtered_spades_contigs.fsa  -1 ../1_Trimmed_Reads/*P1.fq -2 ../1_Trimmed_Reads/*P2.fq -o bowtie & done
## 把参考置于项目目录下
for i in $( ls analysis/assembly) ;do echo $i; cd analysis/assembly/$i/pseudo/$i/ && cp_bowtie.align.pl -i ../../../../../allass.fa  -1 1_Trimmed_Reads/*P1.fq -2 1_Trimmed_Reads/*P2.fq -o bowtie-ass & done
for i in M{1..30};do cd  */analysis/assembly/$i/pseudo/$i/ && cp_bowtie.align.pl -i  ../../../../../../M12_FULLMT.fsa -1 1_Trimmed_Reads/*P1.fq -2 1_Trimmed_Reads/*P2.fq -o bowtie-ass && cd - & done

(3)全部check
for i in $( ls analysis/assembly) ;do echo $i; cd 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
for i in {Mm_J3,Mm_L1,Mm_L2,Mm_L3,Mm_S1,Mm_S2,Mm_S3};do echo $i;cd 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

(4)批量spa/uni组装
---- 1_Trimmed_Reads
for i in `ls analysis/assembly/`;do echo $i;cd analysis/assembly/$i/pseudo/$i/1_Trimmed_Reads && spa- > log1 & done
---- 2_Bowtie_Mapping
for i in `ls analysis/assembly/`;do echo $i;cd analysis/assembly/$i/pseudo/$i/2_Bowtie_Mapping && uni- > log2 & done
for i in M{1..30};do cd  */analysis/assembly/$i/pseudo/$i/2_Bowtie_Mapping && uni- > log2 & done
--- bowtie-ass
for i in `ls analysis/assembly/`;do echo $i;cd analysis/assembly/$i/pseudo/$i/bowtie-ass && spa- > log1 & done
for i in `ls analysis/assembly/`;do echo $i;cd analysis/assembly/$i/pseudo/$i/bowtie-ass && uni- > log2 & done
for i in M{1..30};do cd  */analysis/assembly/$i/pseudo/$i/bowtie-ass && uni- > log2 & done

(5)批量复制gfa
GPXXX]$
--- 3_Spades_Assembly
mkdir pip-spa && for i in `ls analysis/assembly/`;do echo $i;cp analysis/assembly/$i/pseudo/$i/3_Spades_Assembly/spades_iter1/assembly_graph.fastg pip-spa/$i".fastg";done
---- 2_Bowtie_Mapping
mkdir pip-uni && for i in `ls analysis/assembly/`;do echo $i;cp analysis/assembly/$i/pseudo/$i/2_Bowtie_Mapping/uni/002* pip-uni/$i".gfa";done
mkdir pip-uni && for i in M{1..30};do echo $i;cp */analysis/assembly/$i/pseudo/$i/2_Bowtie_Mapping/uni/002* pip-uni/$i".gfa";done
---- bowtie-ass
mkdir spa-ass && for i in `ls analysis/assembly/`;do echo $i;cp analysis/assembly/$i/pseudo/$i/bowtie-ass/spa*/assembly_graph.fastg spa-ass/$i".gfa";done
mkdir uni-ass && for i in `ls analysis/assembly/`;do echo $i;cp analysis/assembly/$i/pseudo/$i/bowtie-ass/uni/002* uni-ass/$i".gfa";done
mkdir uni-ass && for i in M{1..30};do echo $i;cp */analysis/assembly/$i/pseudo/$i/bowtie-ass/uni/002* uni-ass/$i".gfa";done

(6)批量构建bow1文件夹  (上面接(4)批量spa/uni组装)
for i in `ls analysis/assembly/`;do echo $i;mkdir -p analysis/assembly/$i/pseudo/$i/bow1/uni/ ;cp  analysis/assembly/$i/pseudo/$i/2_Bowtie_Mapping/uni/assembly.fasta analysis/assembly/$i/pseudo/$i/bow1/uni/; done

(7)把要循环的命令写进文件 单引号可以正确写入
for i in `ls analysis/assembly/`;do echo $i;echo '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 > log2  ;done &' > analysis/assembly/$i/pseudo/$i/uni100.sh;done

(8)进到目录里执行循环的命令
for i in `ls analysis/assembly/`;do echo $i;cd analysis/assembly/$i/pseudo/$i/ && sh uni100.sh && sh /share/nas1/yuj/script/mitochondrion/assembly/monitor.sh && cd - & done

(9)删除文件夹
for i in `ls analysis/assembly/`;do echo $i;echo 'for i in {1..100} ;do echo $i && rm bow$((i)) -r & done' > analysis/assembly/$i/pseudo/$i/rm100.sh;done
for i in `ls analysis/assembly/`;do echo $i;cd analysis/assembly/$i/pseudo/$i/ && sh rm100.sh && cd -;done

2.4 组装软件

2.4.1 Unicycler混合组装

# ----------------------------2+3----------------------------
unicycler.py -1 *1.fq -2 *2.fq -l *.fasta -o 2+3_correct_unicyc # 校正后
unicycler.py -1 *1.fq -2 *2.fq -l *.fa -o 2+3_unicyc # 前

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

2.4.2 Spades混合组装

# ----------------------------2+3----------------------------
# illumina混合pacbio组装
spades.py -1 *1.fq -2 *2.fq --pacbio *.fa -t 52 -m 400 -o 2+3_hybrid_pacbio # 前
spades.py -1 *1.fq -2 *2.fq --pacbio *.fasta -t 52 -m 400 -o 2+3_correct_hybrid_pacbio # 后

# illumina混合nanopore组装
spades.py -1 *1.fq -2 *2.fq --nanopore *.fasta -t 52 -m 400 -o 2+3_correct_hybrid_nanopore # 后

# ----------------------------2----------------------------
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

2.4.3 Getorganelle

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

# 其他参数   
-R 15 -k 21,45,65,85,105,115 -s 参考

一.动物线粒体
1.常规使用
get_organelle_from_reads.py -1 *1.fq -2 *2.fq -o org -F animal_mt

# 更快的方法
get_organelle_from_reads.py -1 *1.fq -2 *2.fq -o org-fast -F animal_mt --fast            
-w 0.68(可不设置)

2.使用已有fastg组装
get_organelle_from_reads.py -g assembly_graph.fastg -F animal_mt -o org

二.真菌线粒体
get_organelle_from_reads.py -1 *1.fq -2 *2.fq  -R 10 -k 21,45,65,85,105 -F fungus_mt -o fungus_mt_out

2.4.4 Mitoz组装

python3 ~/MitoZ/version_2.4-alpha/release_MitoZ_v2.4-alpha/MitoZ.py all --genetic_code 5 --clade Arthropoda --outprefix ZZZ --thread_number 12 --fastq1 raw.1.fq.gz --fastq2 raw.2.fq.gz --fastq_read_length 150 --insert_size 250 --run_mode 2 --filter_taxa_method 1 --requiring_taxa 'Arthropoda'

python3 MitoZ.py assemble --genetic_code 5 --clade Arthropoda --outprefix test \
--thread_number 8 \
--fastq1 clean.1.fq.gz \
--fastq2 clean.2.fq.gz \
--fastq_read_length 150 \
--insert_size 250 \
--run_mode 2 \
--filter_taxa_method 1 \
--requiring_taxa 'Arthropoda'

# 参数
--genetic_code 为MitoZ设置正确的遗传code，节肢动物使用无脊椎动物的mt_code 5；哺乳动物使用脊椎动物的mt_code 2。
--clade 节肢动物使用'Arthropoda'，脊椎动物使用'Chordata'。
--outprefix 输出文件的前缀
--thread_number 线程数
--fastq1,2 双端测序的原始下机数据
--fastq_read_length 一端测的碱基数目
--insert_size 插入片段的大小
--run_mode 2快速比对；3multi-Kmer方法进行比对，先传参数2，如果结果中蛋白编码基因缺失，再传参3运行
--filter_taxa_method 1对测序数据进行过滤，比如该物种为节肢动物，不是节肢动物的测序reads会被过滤掉；3不进行过滤
--requiring_taxa 确认物种属于哪个类群，来对数据进行过滤

2.4.5 串行运行

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

叶绿体

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