05植物线粒体流程1-组装
二代（主要）
0.创建项目目录，然后创建data目录，把二代和三代数据复制过来（软连接也可，三代要fasta或者fastq格式）
mkdir data/2dai -p
cd data/2dai
ln -s xxx/2.clean_data/*/*_R1_001.fastq.gz
ln -s xxx/2.clean_data/*/*_R2_001.fastq.gz
1.全部数据组装
nohup spades.py --pe1-1 *1.fastq.gz --pe1-2 *2.fastq.gz -o spades -k 55 -m 2000 --only-assembler &
cd spades/
# export ref=/share/nas6/xul/program/mt2/assembly/database/ref.gene.fa
2.根据核心基因比对上的edge，获取相关联的node
nucmer $ref assembly_graph.fastg
show-coords out.delta
vim id # 填入edge的ID
d=50 && perl /share/nas6/xul/program/mt2/assembly/src/filter_node_from_fastg.pl -i id -g assembly_graph.fastg -o $d.target.fastg -d $d & # -d 深度阈值，大于该阈值被保留
随后用bandage打开target.gfa，差不多筛选出包含整个线粒体的序列basic.gfa，不用精修，上传到2代目录下
3.重新组装
data/2dai]$
gfa2fa basic.gfa > basic.gfa.fa && cp_bowtie.align.pl -i basic.gfa.fa -1 *R1*.fastq.gz -2 *R2*.fastq.gz -o bowtie && cd bowtie && nohup spades.py --pe1-1 *1.fq --pe1-2 *2.fq --careful -o spades > logs &
# 装不出来再用 -k 55
nohup unicycler.py -1 *1.fq -2 *2.fq -o uni --kmers 55 > logu & # 先不用吧
简单修剪一下
把graph.gfa上传到data目录
## 删除overlap
/share/nas6/xul/program/bacterial/assembly/Unicycler/remove_spades_gfa_overlap.py -i graph.gfa -o graph.rmlap.gfa
## 改名
perl /share/nas6/xul/program/mt2/assembly/src/change_gfa_name.pl graph.rmlap.gfa > graph.rmlap.rename.gfa
Tips
nad1：5个
nad2：5个
nad5：5个
4.根据三代数据解环，
回到上面三代部分做到第6步，略过第5步
也就是完成校正
cd data
## 4.1 第一个是比对后查看所有的contig覆盖情况
## 找到所有的连接关系：
gfa2fa graph.rmlap.rename.gfa > contig.fa
# fq=3dai/map_basic.fa
# fq=3dai/map_gene.fa
fq=3dai/__extend.map.fa && file=contig.fa && nu=1000000 && minimap2 -t 20 -ax map-pb $file $fq > 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+)[MID]/g;if(($l-$n)/$l>0.7 and (($l-$n)>900)){print}}'|sort -k 4 -n| tail -n $nu > tmp.sam && cat tmp.sam |cut -f 10 |perl -lane 'print ">",++$i;print $F[0]' > extend.map.fa && nucmer $file extend.map.fa -c 1 --maxmatch && show-coords out.delta > out.coords && mummerplot out.delta
感觉这一步看不看都行
4.2需要先完成2代校正3代
## 4.2 第二个是比对后，直接把确定的分支节点拆开，拆分形成新的gfa文件
perl /share/nas6/xul/program/mt2/assembly/src/split_node_with_long_reads2.pl graph.rmlap.gfa ../sample_corrected-reads.fa
或
perl /share/nas6/xul/program/mt2/assembly/src/split_node_with_long_reads2.pl graph.rmlap.rename.gfa ../sample_corrected-reads.fa
# fastalength extend.map.fa > reads_len.txt && perl /share/nas6/xul/program/mt2/assembly/src/split_definite_node_use_coords_in_gfa_file_v2.pl -i1 contig.fa -i2 out.coords -i4 reads_len.txt -i graph.gfa
-------------------------------------------
！！！开始把4.2生成的_final.gfa放到bandage里处理！！！
-------------------------------------------
## 4.3 第三个是针对某个节点，比对后查看它的连接关系，进一步确认用的
## 通过mummer可视化对单独的节点进行判断
grep -A 1 ">12$" __final.gfa.fa > 1.fa && file=1.fa && nu=1000000 && minimap2 -t 20 -ax map-pb $file __extend.map.fa > 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+)[MID]/g;if(($l-$n)/$l>0.7 and (($l-$n)>900)){print}}'|sort -k 4 -n| tail -n $nu > tmp.sam && cat tmp.sam |cut -f 10 |perl -lane 'print ">",++$i;print $F[0]' > extend.map.fa && nucmer __final.gfa.fa extend.map.fa -c 1 --maxmatch && mummerplot out.delta
或
grep -A 1 ">9$" contig.fa > 1.fa && file=1.fa && nu=1000000 && minimap2 -t 20 -ax map-pb $file ../../_extend.map.fa > 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+)[MID]/g;if(($l-$n)/$l>0.7 and (($l-$n)>900)){print}}'|sort -k 4 -n| tail -n $nu > tmp.sam && cat tmp.sam |cut -f 10 |perl -lane 'print ">",++$i;print $F[0]' > extend.map.fa && nucmer contig.fa extend.map.fa -c 1 --maxmatch && mummerplot out.delta
看着也没什么用
## 4.4 第四个是查看某个分支节点的左右两面连接情况，不直接拆分，第二个有时候没法把确定的节点完全拆分，这个时候会用到这个。
## 通过mummer统计连接个数
fastalength extend.map.fa > reads_len.txt && perl /share/nas6/xul/program/mt2/assembly/src/find_contig_link_use_mummer.pl -i1 contig.fa -i2 out.coords -i4 reads_len.txt -p |grep "\-77\-"
或
perl /share/nas6/xul/program/mt2/assembly/src/find_contig_link_use_mummer.pl -i1 __final.gfa.fa -i2 __out.coords -i4 __reads_len.txt -p |grep "\-77\-"
继续调整获得final.gfa
上传至data目录
-------------------------------------------------------------------------------------------------------------------------
## gfa改名，可以用
perl /share/nas6/xul/program/mt2/assembly/src/change_gfa_name.pl final.gfa > final.rename.gfa
## 删除overlap，用
/share/nas6/xul/program/bacterial/assembly/Unicycler/remove_spades_gfa_overlap.py -i final.rename.gfa -o final.rmlap.gfa
## 简化序列，先不用
perl -lane 'if(/^S/){$len=length $F[2];print "$F[0]\t$F[1]\t*\tLN:i:$len"}else{print}' assembly_graph.fastg > target.gfa
调整
5.对gfa序列进行校正：
cd data
perl /share/nas6/xul/program/mt2/assembly/src/pilon_gfa.pl -g final.rename.gfa -i2 2dai/1bowtie/map_pair_hits.* -c cp.fasta -o polish
多环重命名
cat final.raw.fa | grep -v ">" | perl -lane 'BEGIN{$/=undef};@seq=split/\n/,$_; for(sort{length $b <=> length $a}@seq){print ">",++$i," circular\n$_"}' > final.fa
6.创建文件夹
cd ../
prefix=Eremochloa_ophiuroides_Munro_Hack.
mkdir analysis/assembly/$prefix -p
cp data/polish/correctted.gfa.fa analysis/assembly/$prefix/$prefix.fasta
cp data/polish/correctted.gfa analysis/assembly/$prefix/$prefix.gfa
三代（辅助）
1.创建项目目录，然后创建data目录，把二代和三代数据复制过来（软连接也可，三代要fasta或者fastq格式）
mkdir data/3dai -p
cd data/3dai
ln -s xxx/data/*/pass.fq.gz
# export ref=/share/nas6/xul/program/mt2/assembly/database/ref.gene.fa
2.3代数据mapping核心基因
file=pass.fq* && nu=1000000 && minimap2 -t 20 -ax map-pb $ref $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 > 50){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 $ref map_gene.fa -p map_gene && show-coords map_gene.delta > map_gene.coords && mummerplot map_gene.delta &
## 或者，3代数据mapping2代的初步基因组序列，然后map.fa比对核心基因看共线性
file=pass.fq* && nu=1000000 && minimap2 -t 20 -ax map-pb basic.gfa.fa $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 > 50){print}}'|sort -k 4 -n| tail -n $nu > tmp.sam && cat tmp.sam |cut -f 10 |perl -lane 'print ">",++$i;print $F[0]' > map_basic.fa
nucmer $ref map_basic.fa -p map_basic && show-coords map_basic.delta > map_basic.coords && mummerplot map_basic.delta
3.挑选有完整覆盖某个基因的reads进行延伸
grep -A 1 ">60$" map*.fa > test.fa # 把id为60的序列取出来作为延伸基础
4.自动延伸，往往不太🌟
perl /share/nas6/xul/program/mt2/assembly/src/extend.auto.pl -i1 test.fa -i2 map*.fa -polish -re # re:延伸失败时,选择另一条reads来延伸
4.手动延伸
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）。这一步会产生_extend.map.fa文件
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文件
重复循环操作
5.查看组装结果和参考的共线性
nucmer test.fa /share/nas6/xul/program/mt2/assembly/database/ref.gene.fa -c 1 --maxmatch && mummerplot out.delta
6.校正比对上的reads
## 三代自我校正
nohup /share/nas6/zhangxq/biosoft/canu-master/Linux-amd64/bin/canu -correct -p correct -d correct genomeSize=400k useGrid=false -nanopore-raw _extend.map.fa &
nohup /share/nas6/zhangxq/biosoft/canu-master/Linux-amd64/bin/canu -correct -p correct -d correct genomeSize=400k useGrid=false -nanopore-raw map_gene.fa &
nohup /share/nas6/zhangxq/biosoft/canu-master/Linux-amd64/bin/canu -correct -p correct -d correct genomeSize=400k useGrid=false -nanopore-raw map_basic.fa &
## 二代校正三代
需要先完成 二代数据的 3.重新组装
cd correct
perl /share/nas6/xul/program/mt2/assembly/src/TGS_correct_by_NGS.pl -s correct.correctedReads.fasta* -o ./ -i ../../2dai/1bowtie/map_pair_hits.*
或
perl /share/nas6/xul/program/mt2/assembly/src/TGS_correct_by_NGS.pl -s correct.correctedReads.fasta* -o ./ -i ../bowtie/map_pair_hits.*
7.校正组装结果
三代校正
设置初始值
genome=test.fa && cp $genome genome.lgspolish.fa && lgsreads=_extend.map.fa
genome=all.test.fa && cp $genome genome.lgspolish.fa && lgsreads=_extend.map.fa
迭代校正：
for i in {1..6} ;do mv genome.lgspolish.fa $i.fa && genome=$i.fa; minimap2 -ax map-pb -t 20 ${genome} ${lgsreads} |samtools sort -m 2g --threads 20 -o genome.lgs.bam && samtools index genome.lgs.bam && ls `pwd`/genome.lgs.bam > pb.map.bam.fofn && python /share/nas6/xul/soft/NextPolish/NextPolish/lib/nextpolish2.py -g ${genome} -l pb.map.bam.fofn -r clr -p 20 -a -s -o genome.lgspolish.fa ;done
二代校正：
perl /share/nas6/xul/program/fungi/assembly/pipline/src/polish/pilon.pl -i1 genome -i2 map_pair_hits.* -o polish
其他
二代组装
然后尝试下GSAT或者其他软件试试。
先使用getorganelle，数据库用embplant_mt，看组装出来的结构是否完整。
mamba activate getorganelle
/share/nas1/yuj/software/miniconda3/envs/getorganelle/bin/get_organelle_from_reads.py -1 *1.fq -2 *2.fq -o org -F embplant_mt
mamba deactivate
## 路径 /share/nas1/yuj/software/miniconda3/envs/getorganelle/bin/get_organelle_from_reads.py
## 其他参数
-R 15 -k 21,45,65,85,105,115 -s 参考
## 更快的方法
/share/nas1/yuj/software/miniconda3/envs/getorganelle/bin/get_organelle_from_reads.py -1 *1.fq -2 *2.fq -o org -F embplant_mt --fast -w 0.68(可不设置)
## 使用已有fastg组装
/share/nas1/yuj/software/miniconda3/envs/getorganelle/bin/get_organelle_from_reads.py -g assembly_graph.fastg -F embplant_mt -o org
spades -k 55 的纯二代组装
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
三代组装
ls /share/nas2/seqloader/xianliti/GP-20231113-7280_shanxizhongyiyaodaxue_2samples_xianliti2/3dai/Polygala_tenuifolia_Willd/pass.fq.gz > input.fofn
cp /share/nas1/yuj/software/NextDenovo/test_data/run.cfg ./nextdenovo.cfg
/share/nas1/yuj/software/PMAT/bin/PMAT autoMito -i /share/nas2/seqloader/xianliti/GP-20231113-7280_shanxizhongyiyaodaxue_2samples_xianliti2/3dai/Polygala_tenuifolia_Willd/pass.fq.gz -o ./pmat1 -st ONT -g 800m -cs nextDenovo -np /share/nas1/yuj/software/NextDenovo/nextDenovo -cfg nextdenovo.cfg -cp /share/nas6/zhangxq/biosoft/canu-master/Linux-amd64/bin/canu -m