611_宏基因组CRISPR阵列分析¶

# === 方法1：conda安装（推荐）===
conda create -n crispr python=3.8  # 创建环境
conda activate crispr  # 激活环境
conda install -c bioconda -c conda-forge crisprcasfinder  # 安装

# === 方法2：Docker安装（避免依赖问题）===
docker pull unlhcc/crisprcasfinder  # 拉取Docker镜像

# === 验证安装 ===
perl CRISPRCasFinder.pl --help  # 查看帮助

# conda安装
conda install -c bioconda minced  # 从bioconda安装

# 验证安装
minced --help  # 查看帮助，当前版本0.4.2

# === 基本用法：扫描单个基因组 ===
minced genome.fasta output.txt  # 输入基因组，输出到文本文件

# === 输出GFF格式（方便下游分析）===
minced -gff genome.fasta output.txt output.gff  # 同时输出GFF

# === 自定义参数 ===
minced \
    -minNR 3 \          # 最少repeat数（默认3）
    -minRL 23 \         # repeat最短长度（默认23bp）
    -maxRL 47 \         # repeat最长长度（默认47bp）
    -minSL 26 \         # spacer最短长度（默认26bp）
    -maxSL 50 \         # spacer最长长度（默认50bp）
    genome.fasta \      # 输入文件
    output.txt          # 输出文件

# === 批量处理多个基因组 ===
for genome in genomes/*.fasta; do
    name=$(basename "$genome" .fasta)  # 提取文件名
    minced "$genome" "results/${name}_crispr.txt" \
        "results/${name}_crispr.gff"  # 输出到results目录
    echo "完成: $name"  # 打印进度
done

# MinCED输出示例：
Sequence 'contig_1' (125678 bp)
CRISPR 1   Range: 45230 - 45890
POSITION        REPEAT                          SPACER
--------        ------                          ------
45230           GTTTTAGCCCCGATCAGGGGTTCGCCC     AATCGTTCAGCTCCGATCAAGCTGG
45287           GTTTTAGCCCCGATCAGGGGTTCGCCC     TGCAATCGGCTAGCTAACGTTCAAG
45344           GTTTTAGCCCCGATCAGGGGTTCGCCC     CCGATCAAGCTGGAATCGTTCAGCT
45401           GTTTTAGCCCCGATCAGGGGTTCGCCC

# === 标准分析 ===
perl CRISPRCasFinder.pl \
    -in genome.fasta \         # 输入基因组
    -out results_dir \         # 输出目录
    -cas \                     # 同时检测Cas蛋白
    -so /path/to/sel392v2.so \ # CRISPR repeat识别模型
    -cf /path/to/CasFinder-2.0.3 \  # Cas蛋白检测程序路径
    -cpuM 8                    # 线程数

# === 宏基因组模式（重要！）===
perl CRISPRCasFinder.pl \
    -in metagenome_contigs.fasta \  # 宏基因组组装的contig
    -meta \                         # 启用宏基因组模式（关键参数！）
    -mr 23 \                        # repeat最短长度
    -xr 50 \                        # repeat最长长度
    -ms 23 \                        # spacer最短长度
    -xs 50 \                        # spacer最长长度
    -lMin 3 \                       # 最少repeat数量
    -out metaquast_results \        # 输出目录
    -cas \                          # 检测Cas蛋白
    -cpuM 16                        # 线程数

# CRISPRCasFinder能自动分型，结果在JSON文件中
# 查看分型结果
cat results_dir/result.json | python -m json.tool | head -50

# 主要类型：
# Type I（最常见）：含Cas3，切割外来DNA
# Type II（CRISPR-Cas9来源）：含Cas9
# Type III：能切RNA
# Type V：含Cas12
# Type VI：含Cas13，切RNA

# === 提取spacer序列 ===
# 从MinCED结果中提取spacer
grep -v "^POSITION\|^-\|^Sequence\|^CRISPR\|^Repeat" \
    output.txt | awk '{print $3}' | \
    awk 'NF{print ">spacer_"NR"\n"$0}' > spacers.fasta  # 格式化为FASTA

# === 将spacer比对到噬菌体数据库 ===
blastn \
    -query spacers.fasta \           # spacer序列
    -db phage_database \             # 噬菌体数据库
    -out spacer_blast_results.txt \  # 输出结果
    -outfmt 6 \                      # 表格格式
    -evalue 1e-5 \                   # E值阈值
    -max_target_seqs 5               # 每个spacer最多5个匹配

# === 统计spacer匹配到哪些噬菌体 ===
awk '{print $2}' spacer_blast_results.txt | \
    sort | uniq -c | sort -rn | head -20  # 统计匹配频率

#!/bin/bash
# 宏基因组CRISPR阵列完整分析流程

INPUT="metagenome_contigs.fasta"   # 输入文件
OUTDIR="crispr_analysis"           # 输出目录
THREADS=16                          # 线程数

mkdir -p ${OUTDIR}  # 创建输出目录

# === 第1步：用MinCED快速扫描 ===
echo "Step 1: MinCED快速扫描..."
minced -minNR 3 \
    ${INPUT} \
    ${OUTDIR}/minced_results.txt \
    ${OUTDIR}/minced_results.gff

# 统计发现的CRISPR阵列数
n_crispr=$(grep -c "^CRISPR" ${OUTDIR}/minced_results.txt)
echo "MinCED发现 ${n_crispr} 个CRISPR阵列"

# === 第2步：用CRISPRCasFinder深入分析 ===
echo "Step 2: CRISPRCasFinder详细分析..."
perl CRISPRCasFinder.pl \
    -in ${INPUT} \
    -meta \
    -cas \
    -out ${OUTDIR}/ccf_results \
    -cpuM ${THREADS}

# === 第3步：提取和分析spacer ===
echo "Step 3: 提取spacer序列..."
# 从CRISPRCasFinder结果中提取spacer
python3 << 'EOF'
import json  # 导入json模块
import os    # 导入os模块

# 读取CRISPRCasFinder的JSON结果
result_dir = "crispr_analysis/ccf_results"
json_file = os.path.join(result_dir, "result.json")

with open(json_file, 'r') as f:
    data = json.load(f)  # 解析JSON

# 提取所有spacer
spacer_count = 0
with open("crispr_analysis/all_spacers.fasta", 'w') as out:
    for seq_name, seq_data in data.get("Sequences", {}).items():
        for crispr in seq_data.get("Crisprs", []):
            for spacer in crispr.get("Spacers", []):
                spacer_count += 1
                out.write(f">spacer_{spacer_count}_{seq_name}\n")
                out.write(f"{spacer['Sequence']}\n")

print(f"共提取 {spacer_count} 个spacer序列")
EOF

echo "分析完成！结果在 ${OUTDIR}/ 目录"