186_转录组融合基因验证¶

# 下载STAR-Fusion参考数据
wget https://data.broadinstitute.org/Trinity/CTAT_RESOURCE_LIB/GRCh38_gencode_v37_CTAT_lib.tar.gz
tar xzf GRCh38_gencode_v37_CTAT_lib.tar.gz

# 运行STAR-Fusion
STAR-Fusion \
    --genome_lib_dir GRCh38_gencode_v37_CTAT_lib/ \
    --left_fq sample_R1.fastq.gz \
    --right_fq sample_R2.fastq.gz \
    --output_dir star_fusion_output \
    --CPU 16 \
    --FusionInspector validate

# 主要输出
# star-fusion.fusion_predictions.tsv           # 融合预测结果
# star-fusion.fusion_predictions.abridged.tsv  # 精简版
# FusionInspector-validate/                     # 验证结果

# 查看结果
head -5 star_fusion_output/star-fusion.fusion_predictions.abridged.tsv
# FusionName  JunctionReadCount  SpanningFragCount  ...
# TMPRSS2--ERG  45  23  ...

# 安装Arriba
conda install -c bioconda arriba

# 需要先用STAR比对（特殊参数）
STAR \
    --runThreadN 16 \
    --genomeDir star_index/ \
    --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \
    --readFilesCommand zcat \
    --outSAMtype BAM SortedByCoordinate \
    --outSAMunmapped Within \
    --outBAMcompression 0 \
    --chimSegmentMin 10 \
    --chimOutType WithinBAM SoftClip \
    --chimJunctionOverhangMin 10 \
    --chimScoreDropMax 30 \
    --chimScoreJunctionNonGTAG 0 \
    --chimScoreSeparation 1 \
    --alignSJstitchMismatchNmax 5 -1 5 5 \
    --chimSegmentReadGapMax 3 \
    --outFileNamePrefix arriba_star/

# 运行Arriba
arriba \
    -x arriba_star/Aligned.sortedByCoord.out.bam \
    -o fusions.tsv \
    -O fusions.discarded.tsv \
    -a reference.fa \
    -g annotation.gtf \
    -b arriba/database/blacklist_hg38_GRCh38_v2.4.0.tsv.gz \
    -k arriba/database/known_fusions_hg38_GRCh38_v2.4.0.tsv.gz \
    -t arriba/database/known_fusions_hg38_GRCh38_v2.4.0.tsv.gz \
    -p arriba/database/protein_domains_hg38_GRCh38_v2.4.0.gff3

# 可视化融合
draw_fusions.R \
    --fusions=fusions.tsv \
    --alignments=arriba_star/Aligned.sortedByCoord.out.bam \
    --output=fusions.pdf \
    --annotation=annotation.gtf \
    --cytobands=arriba/database/cytobands_hg38_GRCh38_v2.4.0.tsv

import pandas as pd

# 读取各工具结果
star_fusion = pd.read_csv("star_fusion_output/star-fusion.fusion_predictions.abridged.tsv",
                           sep="\t")
arriba = pd.read_csv("fusions.tsv", sep="\t")
# fusioncatcher = pd.read_csv("fusioncatcher_output/final-list_candidate-fusion-genes.txt", sep="\t")

# 标准化融合基因名格式
def normalize_fusion(name, sep='--'):
    genes = sorted(name.replace('--', '__').replace('-', '__').split('__'))
    return '__'.join(genes)

star_fusions = set(star_fusion['#FusionName'].apply(lambda x: normalize_fusion(x)))
arriba_fusions = set()
for _, row in arriba.iterrows():
    fusion = normalize_fusion(f"{row['#gene1']}__{row['gene2']}")
    arriba_fusions.add(fusion)

# 交集
consensus = star_fusions & arriba_fusions
print(f"STAR-Fusion: {len(star_fusions)} fusions")
print(f"Arriba: {len(arriba_fusions)} fusions")
print(f"Consensus: {len(consensus)} fusions")
print("\nConsensus fusions:")
for f in sorted(consensus):
    print(f"  {f}")

# 使用FusionReport或Annovar进一步注释
# 检查是否在已知融合数据库中（COSMIC, ChimerDB, FusionGDB）

# R中的过滤策略
fusions <- read.delim("fusions.tsv", stringsAsFactors = FALSE)

# 过滤条件
filtered <- fusions[
  fusions$split_reads1 + fusions$split_reads2 >= 3 &   # 至少3个split reads
  fusions$discordant_mates >= 2 &                        # 至少2个discordant pairs
  fusions$confidence %in% c("high", "medium") &          # 高/中置信度
  !grepl("read_through", fusions$type) &                 # 排除read-through
  !grepl("healthy", fusions$tags),                       # 排除正常组织常见的
]

cat("过滤前:", nrow(fusions), "fusions\n")
cat("过滤后:", nrow(filtered), "fusions\n")

# 检查已知临床融合
known_fusions <- c("BCR--ABL1", "EML4--ALK", "TMPRSS2--ERG",
                    "EWSR1--FLI1", "SYT--SSX1", "PAX3--FOXO1",
                    "RET--CCDC6", "NTRK1--TPM3", "ROS1--CD74")

filtered$known <- paste0(filtered$gene1, "--", filtered$gene2) %in% known_fusions |
                  paste0(filtered$gene2, "--", filtered$gene1) %in% known_fusions

cat("已知临床融合:", sum(filtered$known), "\n")

# 提取融合断裂点序列用于引物设计
def get_fusion_sequence(gene1_seq, gene2_seq, breakpoint1, breakpoint2):
    """提取融合junction附近的序列"""
    left = gene1_seq[max(0, breakpoint1-150):breakpoint1]
    right = gene2_seq[breakpoint2:breakpoint2+150]
    return left + right

# RT-PCR引物设计指南
primer_design_notes = """
RT-PCR融合基因验证引物设计：
1. 正向引物设计在5'基因的最后1-2个外显子
2. 反向引物设计在3'基因的前1-2个外显子
3. 预期PCR产物大小：200-500bp
4. 引物跨越融合junction
5. 阳性对照：已知含该融合的细胞系
6. 阴性对照：正常组织cDNA
7. Sanger测序确认PCR产物序列

替代验证方法：
- FISH：荧光原位杂交，检测染色体水平的重排
- Western blot：检测融合蛋白（如果有抗体）
- 长读长测序：直接捕获完整融合转录本
"""
print(primer_design_notes)

# STAR-Fusion快速运行
STAR-Fusion --genome_lib_dir ctat_lib/ --left_fq R1.fq.gz --right_fq R2.fq.gz -O output/ --CPU 16

# Arriba快速运行
arriba -x Aligned.bam -o fusions.tsv -a ref.fa -g genes.gtf -b blacklist.tsv

# FusionCatcher运行
fusioncatcher -d data/ -i R1.fq.gz,R2.fq.gz -o fusioncatcher_out/ -p 16

# Arriba可视化
draw_fusions.R --fusions=fusions.tsv --alignments=Aligned.bam --output=plot.pdf