620_多序列比对工具对比¶

# === 三个工具一起安装 ===
conda install -c bioconda mafft muscle clustalo  # 一键安装

# 验证
mafft --version    # MAFFT版本
muscle --version   # MUSCLE版本（v5.x）
clustalo --version # Clustal Omega版本

# === 自动选择策略（推荐）===
mafft --auto input.fasta > aligned.fasta  # 自动选最佳策略

# === 最精确模式（<200条序列）===
mafft --maxiterate 1000 \     # 最多1000次迭代
    --localpair \              # 局部成对比对（L-INS-i模式）
    input.fasta > aligned_linsi.fasta

# 简写
linsi input.fasta > aligned_linsi.fasta  # L-INS-i简写命令

# === 快速模式（>200条序列）===
mafft --retree 2 \             # FFT-NS-2模式
    input.fasta > aligned_fast.fasta

# === 超大数据集（>30000条序列）===
mafft --thread 16 \            # 16线程
    --parttree \               # PartTree算法（最快）
    input.fasta > aligned_huge.fasta

# === 添加新序列到已有比对 ===
mafft --add new_seqs.fasta \   # 新增序列
    --reorder \                 # 重新排序
    existing_alignment.fasta > updated.fasta

# === MUSCLE v5基本用法 ===
muscle -align input.fasta \    # 输入文件
    -output aligned.fasta \    # 输出文件
    -threads 8                 # 线程数

# === 超大数据集用Super5算法 ===
muscle -super5 input.fasta \   # Super5模式（更快）
    -output aligned_super5.fasta \
    -threads 16

# === MUSCLE v3（旧版，部分环境仍在用）===
muscle -in input.fasta \       # 输入
    -out aligned.fasta \       # 输出
    -maxiters 16               # 最大迭代次数

# === 基本用法 ===
clustalo \
    -i input.fasta \           # 输入文件
    -o aligned.fasta \         # 输出文件
    --outfmt=fasta \           # 输出格式
    --threads=8 \              # 线程数
    --force                    # 覆盖已有文件

# === 蛋白质大数据集 ===
clustalo \
    -i proteins.fasta \
    -o aligned_proteins.fasta \
    --full \                   # 完整距离矩阵（更准但更慢）
    --threads=16

# === 添加序列到已有比对 ===
clustalo \
    -i existing_alignment.fasta \
    --profile1=existing_alignment.fasta \
    --profile2=new_sequences.fasta \
    -o updated_alignment.fasta

# === 用trimAl修剪比对 ===
conda install -c bioconda trimal  # 安装trimAl

# 自动修剪
trimal -in aligned.fasta \
    -out trimmed.fasta \
    -automated1  # 自动选择最佳修剪策略

# 按gap比例修剪
trimal -in aligned.fasta \
    -out trimmed.fasta \
    -gt 0.8  # 去除80%以上样本有gap的列

# 查看比对统计
trimal -in aligned.fasta -sgt  # 显示每列的gap比例

# === R语言评估比对质量 ===
library(Biostrings)  # 加载Biostrings
library(ggplot2)     # 加载ggplot2

# 读取比对
aln <- readDNAMultipleAlignment("aligned.fasta", format = "fasta")

# 计算一致性
cons <- consensusMatrix(aln)  # 一致性矩阵
gap_ratio <- cons["-", ] / nrow(aln)  # 每列的gap比例

# 画gap分布图
df <- data.frame(position = 1:length(gap_ratio), gap = gap_ratio)
ggplot(df, aes(x = position, y = gap)) +
    geom_line(color = "steelblue") +
    labs(x = "Position", y = "Gap Ratio", title = "Gap Distribution") +
    theme_minimal()

ggsave("gap_distribution.pdf", width = 12, height = 4)

#!/bin/bash
# 用三个工具分别比对，然后比较结果

INPUT="sequences.fasta"
THREADS=8

# === MAFFT ===
time mafft --auto --thread ${THREADS} ${INPUT} > mafft_aligned.fasta
echo "MAFFT完成"

# === MUSCLE ===
time muscle -align ${INPUT} -output muscle_aligned.fasta -threads ${THREADS}
echo "MUSCLE完成"

# === Clustal Omega ===
time clustalo -i ${INPUT} -o clustalo_aligned.fasta --threads=${THREADS} --force
echo "Clustal Omega完成"

# === 比较比对长度 ===
echo "=== 比对长度比较 ==="
for f in *_aligned.fasta; do
    len=$(head -2 "$f" | tail -1 | wc -c)
    echo "$f: ${len} bp"
done

# === 分别建树比较 ===
for f in *_aligned.fasta; do
    name=$(basename "$f" _aligned.fasta)
    iqtree2 -s "$f" -m MFP -bb 1000 -nt AUTO -pre "tree_${name}" -redo
done