线粒体基因组分析¶

# 查看参考线粒体基因组信息
# 从NCBI下载参考序列
wget "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=NC_012920&rettype=fasta" -O human_mt_ref.fasta

# 查看序列长度
grep -v "^>" human_mt_ref.fasta | tr -d '\n' | wc -c
# 输出: 16569

# 安装NOVOPlasty
git clone https://github.com/ndierckx/NOVOPlasty.git
cd NOVOPlasty

# 准备配置文件 config.txt
cat > config.txt << 'EOF'
Project:
-----------------------
Project name          = MyMito
Type                  = mito
Genome Range          = 12000-22000
K-mer                 = 33
Max memory            = 8
Extended log          = 0
Save assembled reads  = no
Seed Input            = seed.fasta
Extend seed directly  = no
Reference sequence    = 
Variance detection    = no
Heteroplasmy          = 
HP excluded file      = 
Chloroplast sequence  = 

Dataset 1:
-----------------------
Read Length            = 150
Insert size           = 300
Platform              = illumina
Single/Paired         = PE
Combined reads        = 
Forward reads         = sample_R1.fastq.gz
Reverse reads         = sample_R2.fastq.gz
Store Hash            = yes

Optional:
-----------------------
Insert size auto      = yes
Use Quality Scores    = no
Output path           = ./output/
EOF

# 准备seed序列（通常用COX1或cytb基因片段）
# 可以从近缘物种的线粒体基因组中提取
echo ">COX1_seed" > seed.fasta
echo "ATGTTCGGGTGGTTTGGAAACTGACTCATC..." >> seed.fasta

# 运行NOVOPlasty
perl NOVOPlasty4.3.3.pl -c config.txt

# 查看输出
ls output/
# Circularized_assembly_1_MyMito.fasta  # 环化成功的组装结果
# log_MyMito.txt                         # 运行日志

# 安装GetOrganelle
pip install getorganelle
get_organelle_config.py --add animal_mt

# 运行组装
get_organelle_from_reads.py \
    -1 sample_R1.fastq.gz \
    -2 sample_R2.fastq.gz \
    -o getorganelle_output \
    -t 8 \
    -R 20 \
    -k 21,45,65,85,105 \
    -F animal_mt

# GetOrganelle优势：使用多种k-mer自动选择最优组装

# 将reads比对回组装结果
bwa index output/Circularized_assembly_1_MyMito.fasta
bwa mem -t 8 output/Circularized_assembly_1_MyMito.fasta \
    sample_R1.fastq.gz sample_R2.fastq.gz | \
    samtools sort -o mt_mapped.bam
samtools index mt_mapped.bam

# 计算覆盖深度
samtools depth mt_mapped.bam > coverage.txt

# 绘制覆盖曲线（R脚本）
Rscript -e '
cov <- read.table("coverage.txt", col.names=c("chr","pos","depth"))
pdf("mt_coverage.pdf", width=10, height=4)
plot(cov$pos, cov$depth, type="l", xlab="Position", ylab="Depth",
     main="Mitochondrial Genome Coverage")
abline(h=mean(cov$depth), col="red", lty=2)
dev.off()
cat("Mean coverage:", mean(cov$depth), "\n")
cat("Min coverage:", min(cov$depth), "\n")
'

# 检查基因组大小
grep -v "^>" output/Circularized_assembly_1_MyMito.fasta | tr -d '\n' | wc -c

# 方法1：使用MITOS2 web服务器
# 访问 http://mitos2.bioinf.uni-leipzig.de/
# 上传序列，选择遗传密码表，提交

# 方法2：命令行安装MITOS2
conda create -n mitos -c bioconda mitos
conda activate mitos

# 运行MITOS2注释
runmitos.py \
    -i output/Circularized_assembly_1_MyMito.fasta \
    -o mitos_output \
    -r /path/to/mitos2_refdata \
    -c 2 \               # 遗传密码表: 2=脊椎动物线粒体
    --best \             # 只保留最佳注释结果
    --rrna 1 \
    --trna 1 \
    --intron 0

# 查看注释结果
cat mitos_output/result.bed    # BED格式基因坐标
cat mitos_output/result.gff    # GFF格式注释
cat mitos_output/result.fas    # 各基因序列
cat mitos_output/result.seq    # 蛋白序列

# 方法3：使用MitoAnnotator（适合鱼类）
# 访问 http://mitofish.aori.u-tokyo.ac.jp/annotation/input/

# 安装
conda install -c bioconda trnascan-se

# 运行tRNA预测
tRNAscan-SE -O \
    -o trna_output.txt \
    output/Circularized_assembly_1_MyMito.fasta

# -O: 搜索线粒体tRNA（organellar mode）

# 从注释结果提取各基因序列
# 假设已有多个物种的COX1和CYTB序列

# 使用MAFFT进行多序列比对
mafft --auto --adjustdirection COX1_all_species.fasta > COX1_aligned.fasta
mafft --auto CYTB_all_species.fasta > CYTB_aligned.fasta

# 蛋白编码基因推荐使用密码子比对
# 先翻译为蛋白→比对→反推核酸比对
macse -prog alignSequences \
    -seq COX1_all_species.fasta \
    -gc_def 2 \
    -out_NT COX1_NT_aligned.fasta \
    -out_AA COX1_AA_aligned.fasta

# 使用trimAl修剪比对结果
trimal -in COX1_aligned.fasta -out COX1_trimmed.fasta -automated1

# 连接多基因比对（concatenation）
# 使用FASconCAT或AMAS
python3 AMAS.py concat \
    -i COX1_trimmed.fasta CYTB_trimmed.fasta ND1_trimmed.fasta \
    -f fasta \
    -d dna \
    -p partition.txt \
    --out-format fasta \
    --part-format raxml

# 查看partition文件
cat partition.txt
# DNA, COX1 = 1-1536
# DNA, CYTB = 1537-2691
# DNA, ND1 = 2692-3648

# 方法1：IQ-TREE（最大似然法，推荐）
# 自动模型选择 + 分区分析 + ultrafast bootstrap
iqtree2 -s concatenated.fasta \
    -p partition.txt \
    -m MFP \
    --merge rclusterf \
    -B 1000 \
    --alrt 1000 \
    -T AUTO \
    --prefix mt_phylo

# 输出文件：
# mt_phylo.treefile  - 最优树
# mt_phylo.iqtree    - 详细结果报告
# mt_phylo.contree   - 一致性树（含支持值）

# 方法2：RAxML-NG
raxml-ng --all \
    --msa concatenated.fasta \
    --model partition.txt \
    --bs-trees 200 \
    --threads auto \
    --prefix raxml_mt

# 方法3：MrBayes（贝叶斯法）
cat > mrbayes_cmd.nex << 'EOF'
begin mrbayes;
    set autoclose=yes nowarn=yes;
    execute concatenated.nex;
    charset COX1 = 1-1536;
    charset CYTB = 1537-2691;
    partition genes = 2: COX1, CYTB;
    set partition = genes;
    lset nst=6 rates=invgamma;
    unlink statefreq=(all) revmat=(all) shape=(all) pinvar=(all);
    prset applyto=(all) ratepr=variable;
    mcmc ngen=5000000 samplefreq=1000 printfreq=1000 
         nchains=4 nruns=2 savebrlens=yes;
    sump burnin=1250;
    sumt burnin=1250 contype=halfcompat;
end;
EOF

mb mrbayes_cmd.nex

# R脚本：使用ggtree绘制系统发育树
library(ggtree)
library(treeio)
library(ggplot2)

# 读取IQ-TREE结果
tree <- read.iqtree("mt_phylo.contree")

# 基本树形图
p <- ggtree(tree, layout="rectangular") +
    geom_tiplab(size=3) +
    geom_nodelab(aes(label=UFboot), size=2, hjust=-0.2) +
    geom_treescale() +
    theme_tree2()

ggsave("phylo_tree.pdf", p, width=10, height=12)

# 带Bootstrap支持度的彩色树
p2 <- ggtree(tree) +
    geom_tiplab(size=3) +
    geom_nodepoint(aes(color=UFboot), size=2) +
    scale_color_continuous(low="red", high="green",
                          name="Bootstrap") +
    theme_tree2()

ggsave("phylo_tree_bootstrap.pdf", p2, width=10, height=12)

# 使用OGDRAW绘制线粒体基因组环形图
# 在线工具: https://chlorobox.mpimp-golm.mpg.de/OGDraw.html
# 输入GenBank格式文件即可

# 或使用CGView
java -jar cgview.jar \
    -i mt_genome.xml \
    -o mt_genome_map.png \
    -f png

# === 完整流程一键脚本 ===

#!/bin/bash
# 线粒体基因组分析流程
set -e

SAMPLE="species_A"
THREADS=8
R1="raw_data/${SAMPLE}_R1.fastq.gz"
R2="raw_data/${SAMPLE}_R2.fastq.gz"
SEED="reference/COX1_seed.fasta"
CODE=2  # 脊椎动物线粒体遗传密码

# 1. 质控
fastp -i $R1 -o clean_R1.fq.gz \
      -I $R2 -O clean_R2.fq.gz \
      -w $THREADS -h fastp_report.html

# 2. 组装（GetOrganelle）
get_organelle_from_reads.py \
    -1 clean_R1.fq.gz -2 clean_R2.fq.gz \
    -o assembly_output -t $THREADS \
    -R 20 -k 21,45,65,85,105 -F animal_mt

# 3. 注释
runmitos.py \
    -i assembly_output/*.fasta \
    -o annotation_output \
    -r /db/mitos2_refdata -c $CODE --best

# 4. 提取蛋白编码基因
python3 extract_genes.py annotation_output/result.gff \
    assembly_output/*.fasta genes/

# 5. 多序列比对（每个基因单独比对）
for gene in genes/*.fasta; do
    mafft --auto $gene > aligned/$(basename $gene)
done

# 6. 修剪
for aln in aligned/*.fasta; do
    trimal -in $aln -out trimmed/$(basename $aln) -automated1
done

# 7. 连接比对
python3 AMAS.py concat -i trimmed/*.fasta \
    -f fasta -d dna --out-format fasta \
    -p partition.txt

# 8. 建树
iqtree2 -s concatenated.fasta -p partition.txt \
    -m MFP -B 1000 --alrt 1000 -T $THREADS

echo "Pipeline completed!"

GTR+I+G: 最通用，参数最多
HKY+G: 较简单，转换/颠换区分
K2P: 最简单，等碱基频率假设
对于线粒体蛋白编码基因：
- 建议按密码子位点分区
- 第3位点通常需要GTR+G
- 第1、2位点用较简单模型即可
使用ModelFinder自动选择: -m MFP