基因融合检测（STAR-Fusion/Arriba/FusionCatcher）¶

# ===== 数据质控 =====
# FastQC质量检查
fastqc -t 8 -o fastqc_output/ sample_R1.fastq.gz sample_R2.fastq.gz

# MultiQC汇总
multiqc fastqc_output/ -o multiqc_report/

# Fastp质量过滤（推荐用于融合检测前）
fastp \
  --in1 sample_R1.fastq.gz \
  --in2 sample_R2.fastq.gz \
  --out1 clean_R1.fastq.gz \
  --out2 clean_R2.fastq.gz \
  --thread 8 \
  --qualified_quality_phred 20 \
  --length_required 50 \
  --detect_adapter_for_pe \
  --html fastp_report.html

# ===== 参考基因组准备 =====
# 下载参考基因组和注释
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/GRCh38.primary_assembly.genome.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.annotation.gtf.gz
gunzip *.gz

# STAR索引构建（STAR-Fusion和Arriba都需要STAR比对）
STAR --runMode genomeGenerate \
  --runThreadN 16 \
  --genomeDir star_index/ \
  --genomeFastaFiles GRCh38.primary_assembly.genome.fa \
  --sjdbGTFfile gencode.v44.annotation.gtf \
  --sjdbOverhang 150    # read长度 - 1

# ===== STAR-Fusion安装 =====
# Docker方式（推荐）
docker pull trinityctat/starfusion:latest

# Conda方式
conda install -c bioconda star-fusion

# ===== 下载CTAT资源库 =====
# STAR-Fusion需要专门的参考资源库
wget https://data.broadinstitute.org/Trinity/CTAT_RESOURCE_LIB/GRCh38_gencode_v44_CTAT_lib_Oct292023.plug-n-play.tar.gz
tar -xzf GRCh38_gencode_v44_CTAT_lib_Oct292023.plug-n-play.tar.gz

# ===== 运行STAR-Fusion =====
# 方法1：从FASTQ直接运行（推荐）
STAR-Fusion \
  --left_fq clean_R1.fastq.gz \
  --right_fq clean_R2.fastq.gz \
  --genome_lib_dir ctat_genome_lib_build_dir/ \
  --output_dir star_fusion_output/ \
  --CPU 16 \
  --FusionInspector validate \
  --examine_coding_effect \
  --denovo_reconstruct

# 方法2：从已有STAR比对结果运行
# 先用STAR比对（需要特殊参数）
STAR \
  --runThreadN 16 \
  --genomeDir star_index/ \
  --readFilesIn clean_R1.fastq.gz clean_R2.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outReadsUnmapped None \
  --twopassMode Basic \
  --chimSegmentMin 12 \
  --chimJunctionOverhangMin 8 \
  --chimOutJunctionFormat 1 \
  --alignSJDBoverhangMin 10 \
  --alignMatesGapMax 100000 \
  --alignIntronMax 100000 \
  --alignSJstitchMismatchNmax 5 -1 5 5 \
  --chimMultimapScoreRange 3 \
  --chimScoreJunctionNonGTAG -4 \
  --chimMultimapNmax 20 \
  --chimNonchimScoreDropMin 10 \
  --peOverlapNbasesMin 12 \
  --peOverlapMMp 0.1 \
  --alignInsertionFlush Right \
  --alignSplicedMateMapLmin 30 \
  --alignSplicedMateMapLminOverLmate 0 \
  --chimOutType Junctions WithinBAM SoftClip

# 然后运行STAR-Fusion
STAR-Fusion \
  --genome_lib_dir ctat_genome_lib_build_dir/ \
  -J Chimeric.out.junction \
  --output_dir star_fusion_output/

# ===== 查看结果 =====
# 主要输出文件：star-fusion.fusion_predictions.tsv
head star_fusion_output/star-fusion.fusion_predictions.tsv
# 列：FusionName, JunctionReadCount, SpanningFragCount,
#     est_J, est_S, SpliceType, LeftGene, LeftBreakpoint,
#     RightGene, RightBreakpoint, JunctionReads, SpanningFrags, ...

# ===== Arriba安装 =====
conda install -c bioconda arriba

# 或下载预编译版本
wget https://github.com/suhrig/arriba/releases/download/v2.5.1/arriba_v2.5.1.tar.gz
tar -xzf arriba_v2.5.1.tar.gz

# ===== 运行Arriba（与STAR联合） =====
# Arriba提供了一个封装脚本，同时运行STAR和Arriba
arriba/run_arriba.sh \
  star_index/ \
  gencode.v44.annotation.gtf \
  GRCh38.primary_assembly.genome.fa \
  arriba/database/blacklist_hg38_GRCh38_v2.5.0.tsv.gz \
  arriba/database/known_fusions_hg38_GRCh38_v2.5.0.tsv.gz \
  arriba/database/protein_domains_hg38_GRCh38_v2.5.0.gff3 \
  16 \
  clean_R1.fastq.gz clean_R2.fastq.gz

# 或手动分步运行
# 步骤1：STAR比对
STAR \
  --runThreadN 16 \
  --genomeDir star_index/ \
  --readFilesIn clean_R1.fastq.gz clean_R2.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM Unsorted \
  --outSAMunmapped Within \
  --outBAMcompression 0 \
  --chimSegmentMin 10 \
  --chimOutType WithinBAM SoftClip \
  --chimJunctionOverhangMin 10 \
  --chimScoreDropMax 30 \
  --chimScoreJunctionNonGTAG 0 \
  --chimScoreSeparation 1 \
  --alignSJstitchMismatchNmax 5 -1 5 5 \
  --chimSegmentReadGapMax 3 \
  --outStd BAM_Unsorted | \
# 步骤2：Arriba检测融合
arriba \
  -x /dev/stdin \
  -o fusions.tsv \
  -O fusions.discarded.tsv \
  -a GRCh38.primary_assembly.genome.fa \
  -g gencode.v44.annotation.gtf \
  -b arriba/database/blacklist_hg38_GRCh38_v2.5.0.tsv.gz \
  -k arriba/database/known_fusions_hg38_GRCh38_v2.5.0.tsv.gz \
  -t arriba/database/known_fusions_hg38_GRCh38_v2.5.0.tsv.gz \
  -p arriba/database/protein_domains_hg38_GRCh38_v2.5.0.gff3

# ===== Arriba可视化 =====
# 先排序BAM
samtools sort -@ 8 -o Aligned.sortedByCoord.out.bam Aligned.out.bam
samtools index Aligned.sortedByCoord.out.bam

# 绘制融合基因图
draw_fusions.R \
  --fusions=fusions.tsv \
  --output=fusions.pdf \
  --alignments=Aligned.sortedByCoord.out.bam \
  --annotation=gencode.v44.annotation.gtf \
  --cytobands=arriba/database/cytobands_hg38_GRCh38_v2.5.0.tsv \
  --proteinDomains=arriba/database/protein_domains_hg38_GRCh38_v2.5.0.gff3

# ===== 查看结果 =====
# fusions.tsv包含高可信融合
# fusions.discarded.tsv包含被过滤的候选
head -n5 fusions.tsv
# 列：gene1, gene2, strand1, strand2, breakpoint1, breakpoint2,
#     site1, site2, type, direction1, direction2, split_reads1,
#     split_reads2, discordant_mates, coverage1, coverage2,
#     confidence, reading_frame, tags, ...

# ===== FusionCatcher安装 =====
conda install -c bioconda fusioncatcher

# 下载参考数据
fusioncatcher-build -g homo_sapiens -o fusioncatcher_data/

# ===== 运行FusionCatcher =====
fusioncatcher \
  -d fusioncatcher_data/ \
  -i clean_R1.fastq.gz,clean_R2.fastq.gz \
  -o fusioncatcher_output/ \
  -p 16 \
  --skip-blat  # 可选，跳过BLAT加速

# ===== 查看结果 =====
# 主要结果文件
cat fusioncatcher_output/final-list_candidate-fusion-genes.txt
# 按可信度从高到低排列

# GFF格式输出（可用于IGV可视化）
cat fusioncatcher_output/final-list_candidate-fusion-genes.GFF

# ===== R中整合多工具融合检测结果 =====

# 读取各工具结果
# STAR-Fusion
sf <- read.table("star_fusion_output/star-fusion.fusion_predictions.tsv",
                  header = TRUE, sep = "\t", comment.char = "")
sf$fusion_id <- sf$X.FusionName
sf$tool <- "STAR-Fusion"
sf$junction_reads <- sf$JunctionReadCount
sf$spanning_reads <- sf$SpanningFragCount

# Arriba
arriba <- read.table("fusions.tsv", header = TRUE, sep = "\t", comment.char = "")
arriba$fusion_id <- paste0(arriba$X.gene1, "--", arriba$gene2)
arriba$tool <- "Arriba"
arriba$junction_reads <- arriba$split_reads1 + arriba$split_reads2
arriba$spanning_reads <- arriba$discordant_mates

# FusionCatcher
fc <- read.table("fusioncatcher_output/final-list_candidate-fusion-genes.txt",
                  header = TRUE, sep = "\t")
fc$fusion_id <- paste0(fc$Gene_1_symbol.5end_fusion_partner.,
                        "--", fc$Gene_2_symbol.3end_fusion_partner.)
fc$tool <- "FusionCatcher"
fc$junction_reads <- fc$Spanning_unique_reads
fc$spanning_reads <- fc$Spanning_pairs

# 统一格式
normalize_fusion_name <- function(name) {
  genes <- sort(unlist(strsplit(name, "--")))
  paste(genes, collapse = "--")
}

all_fusions <- list(
  `STAR-Fusion` = sapply(sf$fusion_id, normalize_fusion_name),
  Arriba = sapply(arriba$fusion_id, normalize_fusion_name),
  FusionCatcher = sapply(fc$fusion_id, normalize_fusion_name)
)

# 取交集
library(VennDiagram)
venn.diagram(
  x = all_fusions,
  filename = "fusion_venn.png",
  fill = c("red", "blue", "green"),
  alpha = 0.5,
  cat.fontsize = 12
)

# 至少被2个工具检测到的融合
all_names <- unlist(all_fusions)
fusion_counts <- table(all_names)
consensus_fusions <- names(fusion_counts[fusion_counts >= 2])
cat("Consensus fusions (≥2 tools):\n")
print(consensus_fusions)

# ===== 过滤策略 =====
# 1. 读数支持过滤：junction reads ≥ 2, spanning reads ≥ 2
# 2. 已知融合优先：检查COSMIC/ChimerDB/FusionGDB
# 3. 阅读框检查：in-frame融合更可能有生物学意义
# 4. 黑名单过滤：去除已知假阳性（如平行基因读穿readthrough）
# 5. 基因表达过滤：低表达基因的融合可能是噪声

# 已知融合数据库检查
known_fusions <- read.table("known_fusions_database.txt", header = TRUE)
is_known <- consensus_fusions %in% known_fusions$fusion_name
cat("Known fusions:", sum(is_known), "/", length(consensus_fusions), "\n")

# ===== R中融合可视化 =====
# 使用chimeraviz包
BiocManager::install("chimeraviz")
library(chimeraviz)

# 从STAR-Fusion结果创建融合对象
fusions <- import_starfusion(
  "star_fusion_output/star-fusion.fusion_predictions.tsv",
  genome_version = "hg38"
)

# 概览图（所有融合的circos图）
plot_circle(fusions)

# 单个融合的详细视图
# 基因结构图
plot_fusion(fusions[[1]], 
            bamfile = "Aligned.sortedByCoord.out.bam",
            edb = EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86)

# 融合转录本图
plot_fusion_transcript(fusions[[1]], edb = EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86)

# ===== IGV手动验证 =====
# 在IGV中加载BAM文件，跳转到断点位置
# 查看是否有split reads和discordant pairs支持
# 检查reads的mapping quality和base quality

# ===== 实验验证方法 =====
# 1. RT-PCR验证
#    设计引物跨越融合断点
#    正向引物在5'基因，反向引物在3'基因
#    PCR产物测序确认断点序列

# 2. FISH（荧光原位杂交）
#    用于验证基因组水平的重排

# 3. Western Blot
#    检测融合蛋白的表达（需要针对融合特异的抗体）

# ===== 设计RT-PCR引物验证融合 =====
# 从融合序列中提取断点两侧序列
# 使用Primer3设计跨断点引物

# 示例：EML4-ALK融合
# 5'端引物（在EML4区域）: FORWARD_PRIMER
# 3'端引物（在ALK区域）: REVERSE_PRIMER
# PCR产物应为融合特异性的（正常样本无产物）

# 使用bedtools提取融合区域序列
bedtools getfasta -fi GRCh38.fa -bed fusion_breakpoints.bed -fo fusion_sequences.fa

# ===== 功能注释 =====
# 检查融合蛋白的结构域保留情况
# 5'基因保留了什么功能域？3'基因保留了什么？

# 使用Pfam/InterPro注释蛋白域
library(biomaRt)
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# 查询融合伴侣基因的蛋白域
gene_domains <- getBM(
  attributes = c("hgnc_symbol", "pfam", "pfam_start", "pfam_end", "interpro_description"),
  filters = "hgnc_symbol",
  values = c("EML4", "ALK"),
  mart = ensembl
)

# ===== 临床意义注释 =====
# 检查已知药物靶点
# OncoKB, CIViC, COSMIC融合数据库
# oncokb.org - 提供Level of Evidence
# civicdb.org - 社区驱动的证据数据库

# 检查阅读框
# In-frame融合：可能产生功能性嵌合蛋白
# Out-of-frame融合：通常产生截短蛋白或无义介导的mRNA降解(NMD)

# ===== 融合基因表达定量 =====
# 利用spanning和junction reads估计融合转录本的丰度
# 与野生型转录本的丰度比较

#!/bin/bash
# ============================================
# Gene Fusion Detection Pipeline
# ============================================

# 配置
THREADS=16
GENOME="GRCh38.primary_assembly.genome.fa"
GTF="gencode.v44.annotation.gtf"
STAR_INDEX="star_index/"
CTAT_LIB="ctat_genome_lib_build_dir/"
R1="clean_R1.fastq.gz"
R2="clean_R2.fastq.gz"
SAMPLE="sample01"
OUTDIR="fusion_results/${SAMPLE}"
mkdir -p $OUTDIR/{star_fusion,arriba,summary}

# === 1. STAR比对（共用） ===
STAR --runThreadN $THREADS \
  --genomeDir $STAR_INDEX \
  --readFilesIn $R1 $R2 \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMunmapped Within \
  --chimSegmentMin 10 \
  --chimOutType Junctions WithinBAM SoftClip \
  --chimJunctionOverhangMin 10 \
  --chimScoreDropMax 30 \
  --chimScoreJunctionNonGTAG 0 \
  --chimScoreSeparation 1 \
  --alignSJstitchMismatchNmax 5 -1 5 5 \
  --chimSegmentReadGapMax 3 \
  --chimMultimapNmax 20 \
  --twopassMode Basic \
  --outFileNamePrefix $OUTDIR/star_

# === 2. STAR-Fusion ===
STAR-Fusion \
  --genome_lib_dir $CTAT_LIB \
  -J $OUTDIR/star_Chimeric.out.junction \
  --output_dir $OUTDIR/star_fusion/ \
  --examine_coding_effect

# === 3. Arriba ===
arriba \
  -x $OUTDIR/star_Aligned.sortedByCoord.out.bam \
  -o $OUTDIR/arriba/fusions.tsv \
  -O $OUTDIR/arriba/fusions.discarded.tsv \
  -a $GENOME \
  -g $GTF \
  -b arriba/database/blacklist_hg38_GRCh38_v2.5.0.tsv.gz \
  -k arriba/database/known_fusions_hg38_GRCh38_v2.5.0.tsv.gz \
  -p arriba/database/protein_domains_hg38_GRCh38_v2.5.0.gff3

# Arriba可视化
samtools index $OUTDIR/star_Aligned.sortedByCoord.out.bam
draw_fusions.R \
  --fusions=$OUTDIR/arriba/fusions.tsv \
  --output=$OUTDIR/arriba/fusions.pdf \
  --alignments=$OUTDIR/star_Aligned.sortedByCoord.out.bam \
  --annotation=$GTF

# === 4. 汇总 ===
echo "=== STAR-Fusion Results ===" > $OUTDIR/summary/fusion_summary.txt
wc -l $OUTDIR/star_fusion/star-fusion.fusion_predictions.tsv >> $OUTDIR/summary/fusion_summary.txt
echo "=== Arriba Results ===" >> $OUTDIR/summary/fusion_summary.txt
wc -l $OUTDIR/arriba/fusions.tsv >> $OUTDIR/summary/fusion_summary.txt

echo "Pipeline complete for $SAMPLE"