三代测序数据处理（Nanopore/PacBio）¶

# === Dorado basecalling ===
# 新版Dorado（v0.7+）支持简化的"模型复合体"语法，自动根据数据选择匹配模型
# 推荐使用简化语法：fast / hac / sup，Dorado自动检测化学版本

# 基本basecalling（推荐方式：简化模型选择）
dorado basecaller sup \
  /data/pod5_files/ \
  --device cuda:all \
  > basecalled.bam

# 输出FASTQ格式
dorado basecaller sup \
  /data/pod5_files/ \
  --device cuda:all \
  --emit-fastq \
  > basecalled.fastq

# 输出BAM + 同时比对到参考基因组
dorado basecaller sup \
  /data/pod5_files/ \
  --device cuda:0 \
  --reference /ref/GRCh38.fa \
  > basecalled_aligned.bam

# 带修饰碱基检测（在模型复合体中用逗号添加修饰类型）
dorado basecaller sup,5mCG_5hmCG \
  /data/pod5_files/ \
  --device cuda:all \
  > basecalled_mods.bam

# Duplex calling（超高准确度）
dorado duplex sup \
  /data/pod5_files/ \
  --device cuda:all \
  > duplex.bam

# 也可使用完整模型名称（适合指定特定版本）
# dorado basecaller dna_r10.4.1_e8.2_400bps_sup@v5.2.0 /data/pod5_files/

# 从BAM提取FASTQ
samtools fastq basecalled.bam > basecalled.fastq

# === 旧版Guppy（已停止更新，仅供参考）===
# 注意：Guppy已被Dorado完全替代，R9.4.1和RNA002模型在Dorado v1.0.0中已弃用
# 新项目请统一使用Dorado
# guppy_basecaller \
#   -i /data/fast5_files/ \
#   -s /output/basecalled/ \
#   -c dna_r10.4.1_e8.2_400bps_sup.cfg \
#   --device cuda:0 \
#   --records_per_fastq 0 \
#   --compress_fastq

# === 生成CCS/HiFi reads ===
# 从subreads.bam生成HiFi reads
ccs subreads.bam hifi_reads.bam \
  --min-rq 0.99 \          # 最低准确度（Q20）
  --min-passes 3 \         # 最少通过次数
  --num-threads 32 \
  --log-level INFO \
  --report-file ccs_report.txt

# 查看CCS报告
cat ccs_report.txt

# 提取FASTQ
samtools fastq hifi_reads.bam > hifi_reads.fastq.gz

# 统计
samtools stats hifi_reads.bam > hifi_stats.txt

# === 使用pbmm2比对（PacBio推荐）===
pbmm2 align /ref/GRCh38.fa hifi_reads.bam aligned.bam \
  --sort \
  --preset HiFi \
  --sample sample1 \
  -j 16

# === 使用actc生成HiFi（新版本仪器）===
# Revio/Sequel IIe直接输出HiFi，无需手动CCS
# 输出：*.hifi_reads.bam

# === NanoPlot ===
# 从FASTQ质控
NanoPlot --fastq basecalled.fastq.gz \
  -o nanoplot_output/ \
  -t 8 \
  --title "Sample1 Nanopore QC" \
  --plots hex dot \
  --N50 \
  --loglength

# 从BAM质控
NanoPlot --bam basecalled_aligned.bam \
  -o nanoplot_bam/ \
  -t 8 \
  --title "Sample1 Aligned" \
  --plots hex dot

# 从sequencing_summary质控（更快）
NanoPlot --summary sequencing_summary.txt \
  -o nanoplot_summary/ \
  -t 8

# === NanoQC（位置级质量）===
nanoQC -o nanoqc_output/ basecalled.fastq.gz

# === NanoStat（纯数字统计）===
NanoStat --fastq basecalled.fastq.gz \
  -t 8 \
  -n nanostats_results.txt

# === chopper（质量过滤+长度过滤）===
# 替代NanoFilt（更快）
gunzip -c basecalled.fastq.gz | \
  chopper --quality 10 \
          --minlength 1000 \
          --maxlength 100000 \
          --threads 8 | \
  gzip > filtered.fastq.gz

# === PacBio HiFi QC ===
# 使用samtools stats
samtools stats hifi_reads.bam > hifi_stats.txt

# 或使用NanoPlot（也支持PacBio）
NanoPlot --fastq hifi_reads.fastq.gz \
  -o nanoplot_hifi/ \
  -t 8 \
  --title "PacBio HiFi QC"

# 基本统计
echo "Total reads: $(samtools view -c hifi_reads.bam)"
echo "Total bases: $(samtools stats hifi_reads.bam | grep 'bases mapped' | cut -f3)"

# === Nanopore数据比对 ===
minimap2 -ax map-ont \
  -t 16 \
  --secondary=no \
  --MD \
  /ref/GRCh38.fa \
  filtered.fastq.gz | \
  samtools sort -@ 8 -o sample_ont_sorted.bam -

samtools index sample_ont_sorted.bam

# === PacBio HiFi比对 ===
minimap2 -ax map-hifi \
  -t 16 \
  --secondary=no \
  --MD \
  /ref/GRCh38.fa \
  hifi_reads.fastq.gz | \
  samtools sort -@ 8 -o sample_hifi_sorted.bam -

samtools index sample_hifi_sorted.bam

# === PacBio CLR比对 ===
minimap2 -ax map-pb \
  -t 16 \
  /ref/GRCh38.fa \
  clr_reads.fastq.gz | \
  samtools sort -@ 8 -o sample_clr_sorted.bam -

# === 长读长RNA比对（如Nanopore direct RNA）===
minimap2 -ax splice \
  -t 16 \
  --junc-bed /ref/known_junctions.bed \
  /ref/GRCh38.fa \
  direct_rna.fastq.gz | \
  samtools sort -@ 8 -o rna_sorted.bam -

# === 比对质量统计 ===
samtools flagstat sample_ont_sorted.bam
samtools coverage sample_ont_sorted.bam > coverage.txt

# 查看比对率和平均覆盖度
echo "Mapping rate:"
samtools flagstat sample_ont_sorted.bam | grep "mapped ("

# 深度统计
mosdepth -t 8 --by 1000 sample_depth sample_ont_sorted.bam

# === Flye（Nanopore/PacBio通用）===
flye --nano-hq filtered.fastq.gz \    # Nanopore Q20+
  --out-dir flye_assembly/ \
  --threads 32 \
  --genome-size 5m                     # 预估基因组大小

# PacBio HiFi
flye --pacbio-hifi hifi_reads.fastq.gz \
  --out-dir flye_hifi/ \
  --threads 32 \
  --genome-size 3g

# === hifiasm（PacBio HiFi最佳选择，v0.25+支持ONT组装）===
hifiasm -o sample_asm \
  -t 32 \
  hifi_reads.fastq.gz

# 转换GFA到FASTA
awk '/^S/{print ">"$2"\n"$3}' sample_asm.bp.p_ctg.gfa > assembly.fa

# HiFi + ultralong ONT联合组装（T2T级）
# hifiasm -o asm_hybrid -t 32 --ul ont_ultralong.fastq.gz hifi_reads.fastq.gz

# Trio模式（需父母短reads，先用yak计k-mer）
# yak count -k31 -b37 -t16 -o pat.yak paternal.fq.gz
# yak count -k31 -b37 -t16 -o mat.yak maternal.fq.gz
# hifiasm -o asm_trio -t 32 -1 pat.yak -2 mat.yak hifi_reads.fastq.gz

# === 组装校正（Nanopore数据可能需要）===
# 使用Medaka（Nanopore自校正）
medaka_consensus \
  -i filtered.fastq.gz \
  -d flye_assembly/assembly.fasta \
  -o medaka_output/ \
  -t 16 \
  -m r1041_e82_400bps_sup_v4.3.0

# === 组装质量评估 ===
# QUAST
quast assembly.fa \
  -r /ref/reference.fa \
  -o quast_results/ \
  -t 16

# BUSCO完整性
busco -i assembly.fa \
  -m genome \
  -l bacteria_odb10 \
  -c 16 \
  -o busco_results

# === Sniffles2（Nanopore/PacBio SV检测）===
sniffles --input sample_ont_sorted.bam \
  --vcf sv_calls.vcf \
  --reference /ref/GRCh38.fa \
  --threads 16 \
  --minsvlen 50 \
  --mapq 20

# === cuteSV（另一个SV caller）===
cuteSV sample_ont_sorted.bam \
  /ref/GRCh38.fa \
  cutesv_output.vcf \
  ./cutesv_work_dir/ \
  --threads 16 \
  --sample sample1 \
  --genotype \
  --min_size 50 \
  --max_size 100000

# === PBSV（PacBio官方SV工具）===
pbsv discover sample_hifi_sorted.bam sample.svsig.gz
pbsv call /ref/GRCh38.fa sample.svsig.gz sv_output.vcf \
  --ccs --num-threads 16

# ===== 完整Nanopore数据处理流水线 =====

SAMPLE="sample1"
REF="/ref/GRCh38.fa"
THREADS=32

# 1. Basecalling（使用简化模型选择，Dorado自动检测化学版本）
dorado basecaller sup \
  /data/pod5/${SAMPLE}/ \
  --device cuda:all \
  > ${SAMPLE}_basecalled.bam

# 转FASTQ
samtools fastq ${SAMPLE}_basecalled.bam | gzip > ${SAMPLE}.fastq.gz

# 2. 质控
NanoPlot --fastq ${SAMPLE}.fastq.gz -o ${SAMPLE}_qc/ -t ${THREADS} --N50

# 3. 过滤
gunzip -c ${SAMPLE}.fastq.gz | \
  chopper -q 10 --minlength 1000 -t ${THREADS} | \
  gzip > ${SAMPLE}_filtered.fastq.gz

# 4. 比对
minimap2 -ax map-ont -t ${THREADS} --MD ${REF} ${SAMPLE}_filtered.fastq.gz | \
  samtools sort -@ 8 -o ${SAMPLE}_sorted.bam -
samtools index ${SAMPLE}_sorted.bam

# 5. 比对统计
samtools flagstat ${SAMPLE}_sorted.bam
mosdepth -t 8 ${SAMPLE}_depth ${SAMPLE}_sorted.bam

# 6. 组装（可选）
flye --nano-hq ${SAMPLE}_filtered.fastq.gz \
  --out-dir ${SAMPLE}_assembly/ \
  --threads ${THREADS} \
  --genome-size 5m

# 7. SV检测（可选）
sniffles --input ${SAMPLE}_sorted.bam --vcf ${SAMPLE}_sv.vcf \
  --reference ${REF} --threads ${THREADS}

# ===== 完整PacBio HiFi流水线 =====
# 1. CCS（如果仪器未自动生成HiFi）
ccs subreads.bam hifi.bam --min-rq 0.99 --num-threads ${THREADS}

# 2. 比对
minimap2 -ax map-hifi -t ${THREADS} ${REF} hifi.fastq.gz | \
  samtools sort -@ 8 -o hifi_sorted.bam -
samtools index hifi_sorted.bam

# 3. 组装
hifiasm -o asm -t ${THREADS} hifi.fastq.gz