1️.MITOS注释

## for循环步骤
for i in $(ls analysis/annotation);do echo $i;cd analysis/annotation/$i &&  mkdir mitos2 && cd -;done
for i in $(ls analysis/annotation);do echo $i;cd analysis/annotation/$i && micromamba run -n mitos2 runmitos.py -i $i"_FULLMT.fsa"  -o mitos2/ -c 4 -R /share/nas1/yuj/software/mitos2_refdata/0.3_ref -r refseq63f --maxtrnaovl 100 --maxrrnaovl 100 --rrna 2 & done # 注意，这里是refseq63f

## 步骤
cd analysis/annotation/xxx  
mkdir mitos2
micromamba run -n mitos2 runmitos.py -i *.fsa -o mitos2/ -c 4 -R /share/nas1/yuj/software/mitos2_refdata/0.3_ref -r refseq63f --maxtrnaovl 100 --maxrrnaovl 100 --best --rrna 2

二.MITOS注释修改

2.1 解析注释结果

## for循环步骤
1）转换结果
for i in $(ls analysis/annotation);do echo $i;cd analysis/annotation/$i && python3 /share/nas1/yuj/script/mitochondrion/annotation/mt_mitos2info_table.py -i mitos2/result.mitos -fa $i"_FULLMT.fsa"  -o gene.info &&  python3 /share/nas1/yuj/script/fungi/annotation/fungi_parse_info_table.py -i gene.info -o gene.annotation.info2 -n 4 &&  cd - & done

2）看一下基因数量
for i in $(ls analysis/annotation);do echo $i;grep "Total" analysis/annotation/$i/gene.annotation.info2;done

2.2 修改注释

2.2.1 共线性好，用一个ref（不推荐）

1.总览
count=0 && for i in $(grep 'CDS' gene.annotation.info2 | awk '{print $3 "_:_" $2}'); do count=$((count+1)) &&  echo "" && echo "###############" $count ":"  $(echo "$i" | awk -F "_:_" '{print $1}') "###############" && fu_add.py  -i *.fsa -p $(echo "$i" | awk -F "_:_" '{print $2}') -f ;done  > sample_ann.log

count=0 && for i in $(grep 'CDS' *.gbk.ann | awk '{print $3 "_:_" $2}'); do count=$((count+1)) &&  echo "" && echo "###############" $count ":"  $(echo "$i" | awk -F "_:_" '{print $1}') "###############" && fu_add.py  -i *.fasta  -p $(echo "$i" | awk -F "_:_" '{print $2}') -f ;done  > ref_ann.log

2.提取cds序列，挑选一个作为参考
cd ref_ann
python3 /share/nas1/yuj/script/fungi/annotation/from_gbk_get_cds.py -i gbk/ -o out  -d
ir.py -s7 -n 4 -i out/cds/cds_OP453745.1.fasta -o cds_OP453745.1.pep && cp out/cds/cds_OP453745.1.fasta ./

tblastn -query cds_NC_022933.1.pep -subject  ../*.fsa -outfmt 6 | cut -f 1,7,8,9,10 > cdsNC_022933.1.pep.log
blastn -query cds_NC_022933.1.fasta -subject  ../*.fsa -outfmt 6 | cut -f 1,7,8,9,10 > cdsNC_022933.1.fa.log
cd ../

3.进一步处理这些比对出的位置
## for循环直接全部处理了
cd ref_ann && for i in cds*log; do for j in {ATP6,ATP8,ATP9,COX1,COX2,COX3,CYTB,ND1,ND2,ND3,ND4,ND4L,ND5,ND6,RPS3}; do  echo ""; echo $j"-----------------------------------------"; grep -i  $j $i ; done > ../${i%.log}.txt ; done && cd ..
1	1481	6276	7756
1	1481	6276	7756
从结果里挑选最符合的

2.2.2 共线性有点差，把所有参考都用上（推荐）

mkdir ann_tmp
1.提取出cds及对应蛋白序列
python3 /share/nas1/yuj/script/fungi/annotation/from_gbk_get_cds.py -i ref_ann/gbk/ -o ref_ann/out  -d
cp ref_ann/out/cds/*.fasta ann_tmp/ && python3 /share/nas1/yuj/script/mitochondrion/phytree/mt_extract2mafft_V1.5.py -i ann_tmp/ -o1 gene/extract -o2 gene/mafft -c2
停顿几秒等比对完
cp *.fsa ann_tmp && cd ann_tmp
for i in cds*.fasta;do ir.py -s7 -n 4 -i $i -o ${i%.fasta}.pep;done

2.利用核酸和蛋白序列blast比对
for i in cds*.fasta;do blastn -query $i  -subject *.fsa -outfmt 6 -max_hsps 1 > $i.log;done
for i in cds*.fasta;do tblastn -query ${i%.fasta}.pep -subject *.fsa -outfmt 6 -max_hsps 1 > ${i%.fasta}.pep.log;done
for i in cds_*.log;do echo $i;cut -f 1,7,8,9,10 $i > ${i%.log}.txt;done

3.进一步处理这些比对出的位置
## for循环直接全部处理了
# for gene in {ATP6,ATP8,COX1,COX2,COX3,CYTB,ND1,ND2,ND3,ND4L,ND4,ND5,ND6}; do echo $gene;for i in cds*txt; do grep -i "\b$gene\b" $i ; done | cut -f 2,3,4,5 | sort -k 3 -n ;done > ../allcds.blast.txt && cd .. # 精准匹配

for gene in {ATP6,ATP8,COX1,COX2,COX3,CYTB,ND1,ND2,ND3,ND4L,ND4,ND5,ND6}; do echo $gene;for i in cds*txt; do grep -i  $gene $i ; done | cut -f 2,3,4,5 | sort -k 3 -n ;done > ../allcds.blast.txt && cd .. # 真菌提取的基因名字带有其他信息，不能用精准匹配
1	1481	6276	7756
1	1481	6276	7756
从结果里挑选最符合的

## 比如这个是查找nd4基因
for i in cds*txt;do grep -i -v nd4l $i| grep -i nd4;done | cut -f 2,3,4,5 | sort -k 3

2.3 回到线粒体注释的2.3章节往下

count=0 && for i in $(grep 'CDS' gene.annotation.info | awk '{print $3 "_:_" $2}'); do count=$((count+1)) && echo "" && echo "###############" $count ":" $(echo "$i" | awk -F "_:_" '{print $1}') "###############" && fu_add.py -i *.fsa -p $(echo "$i" | awk -F "_:_" '{print $2}') -f ;done > sample_ann.log

三.高级分析

3.1 pi（叶绿体程序）

## 完整程序
perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/script/pi.analysis/pi.no.ir.analysis.pl -i gbk/  -p Huperzia_crispata -o ./

## 仅重新绘图
python3 /share/nas1/yuj/script/fungi/advance/pi/fu_plot.pi.no.ir.py -i pi.stat.xls -o pi.gene.png