手把手教你构建人+EBV融合基因组：10X单细胞转录组定量参考文件制作全攻略

在单细胞转录组研究中，参考基因组的质量直接决定了基因定量的准确性。对于研究 Epstein-Barr 病毒（EBV）相关疾病（如 Burkitt 淋巴瘤、鼻咽癌等）的科研人员来说，仅使用人类参考基因组往往无法捕捉病毒基因的表达信息。今天，我们就来详细介绍如何构建包含人类基因组与 EBV 基因组的融合参考文件，为 10X 单细胞转录组分析提供精准“地图”。

image.png

单细胞转录组、空间转录组、转录组、宏基因组、16S扩增子、基因组、ChIP-seq、miRNA-seq等组学、TCGA/GEO数据挖掘、研究学习交流《群》：

为什么需要人+EBV融合参考基因组？

EBV 作为一种常见的疱疹病毒，在感染人体后可潜伏在细胞内，与多种恶性肿瘤的发生发展密切相关。在单细胞转录组测序中，若参考基因组仅包含人类序列，EBV 基因的 reads 会因无法匹配而被过滤，导致病毒-宿主相互作用的关键信息丢失。
构建人+EBV 融合参考基因组，能同时定量人类基因和 EBV 基因的表达，为解析病毒感染机制、筛选疾病标志物提供更全面的视角。

核心工具与版本选择

本次构建基于 10X Genomics 官方工具 Cell Ranger，关键版本选择如下：

• 人类基因组版本：GRCh38（primary assembly，基于 Ensembl Release 109）
• 基因注释文件：GENCODE v44（primary assembly 注释）
• EBV 基因组：参考序列 GCF_002402265.1（ASM240226v1）
• Cell Ranger 版本：9.0.1（支持自定义参考基因组构建）

详细构建步骤

一、准备工作：创建工作目录

首先创建用于存储源文件和结果的目录，避免文件混乱：

# 定义基因组版本和构建目录
genome="GRCh38"
version="hsa-EBV"
build="GRCh38-GENCODEv44_build"
source="reference_sources"

# 创建目录
mkdir -p "$build" "$source"

二、下载核心参考文件

参考文件的选择直接影响后续分析准确性，需重点关注版本兼容性：

1. 人类基因组 FASTA 文件

选择 Ensembl Release 109 的 GRCh38 主装配序列，而非更新的 Release 110。原因是 Release 110 将假常染色体区域（PAR）去遮蔽，可能导致基因映射歧义，而 Release 109 保持了 PAR 区域的遮蔽状态，更适合定量分析：

fasta_url="http://ftp.ensembl.org/pub/release-109/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz"
fasta_in="${source}/Homo_sapiens.GRCh38.dna.primary_assembly.fa"

# 若文件未下载，则通过 curl 下载并解压
if [ ! -f "$fasta_in" ]; then
    curl -sS "$fasta_url" | zcat > "$fasta_in"
fi

2. 人类基因注释 GTF 文件

选择 GENCODE Release 44 的主装配注释文件，包含更全面的基因结构信息（外显子、转录本等）：

gtf_url="http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.primary_assembly.annotation.gtf.gz"
gtf_in="${source}/gencode.v44.primary_assembly.annotation.gtf"

# 若文件未下载，则通过 curl 下载并解压
if [ ! -f "$gtf_in" ]; then
    curl -sS "$gtf_url" | zcat > "$gtf_in"
fi

3. EBV 基因组及注释文件

从 NCBI 下载 EBV 参考基因组（GCF_002402265.1）及其基因注释 GTF 文件，存储于 reference_sources 目录。

wget -O "reference_sources/GCF_002402265.1_ASM240226v1_genomic.fna.gz" 
wget -O "GCF_002402265.1_ASM240226v1_genomic.gtf.gz" 
curl -sS https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/402/265/GCF_002402265.1_ASM240226v1/GCF_002402265.1_ASM240226v1_genomic.fna.gz | zcat > "reference_sources/GCF_002402265.1_ASM240226v1/GCF_002402265.1_ASM240226v1_genomic.fna"
curl -sS https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/402/265/GCF_002402265.1_ASM240226v1/GCF_002402265.1_ASM240226v1_genomic.fna.gz | zcat > "reference_sources/GCF_002402265.1_ASM240226v1/GCF_002402265.1_ASM240226v1_genomic.fna"

if [ ! -f "$gtf_in" ]; then
    curl -sS "$gtf_url" | zcat > "$gtf_in"
fi

三、格式标准化：让文件“兼容”Cell Ranger

不同数据库的文件格式存在差异，需统一格式以确保 Cell Ranger 正确识别：

1. 修正人类基因组 FASTA 头部信息

Ensembl 与 GENCODE 的 FASTA 头部格式不同，需将 Ensembl 格式转换为 GENCODE 兼容格式（如染色体名称添加“chr”前缀）：

# 输入 FASTA 头部（Ensembl 格式）：
# >1 dna:chromosome chromosome:GRCh38:1:1:248956422:1 REF
# 输出 FASTA 头部（GENCODE 格式）：
# >chr1 1

fasta_modified="$build/$(basename "$fasta_in").modified"
cat "$fasta_in" /
    | sed -E 's/^>(/S+).*/>/1 /1/' /  # 简化头部信息
    | sed -E 's/^>([0-9]+|[XY]) />chr/1 /' /  # 常染色体/性染色体加“chr”
    | sed -E 's/^>MT />chrM /' /  # 线粒体染色体命名为 chrM
    > "$fasta_modified"

2. 优化基因注释 GTF 文件

GTF 文件中基因 ID 包含版本号（如 ENSG00000223972.5），需去除版本号并保留版本信息，同时过滤低质量注释：

# 去除 ID 中的版本号（如 ENSG00000223972.5 → ENSG00000223972），并添加版本标签
gtf_modified="$build/$(basename "$gtf_in").modified"
ID="(ENS(MUS)?[GTE][0-9]+)/.([0-9]+)"  # 匹配 Ensembl 基因/转录本/外显子 ID 格式
cat "$gtf_in" /
    | sed -E 's/gene_id "'"$ID"'";/gene_id "/1"; gene_version "/3";/' /
    | sed -E 's/transcript_id "'"$ID"'";/transcript_id "/1"; transcript_version "/3";/' /
    | sed -E 's/exon_id "'"$ID"'";/exon_id "/1"; exon_version "/3";/' /
    > "$gtf_modified"

3. 筛选高质量基因注释

仅保留蛋白编码基因、长链非编码 RNA（lncRNA）及免疫相关基因（如 IG、TR 家族），排除通读转录本（readthrough transcript）以减少干扰：

# 定义需保留的基因类型
BIOTYPE_PATTERN="(protein_coding|lncRNA|IG_C_gene|IG_D_gene|IG_J_gene|TR_C_gene|TR_D_gene|TR_J_gene|TR_V_gene)"
GENE_PATTERN="gene_type /"${BIOTYPE_PATTERN}/""
TX_PATTERN="transcript_type /"${BIOTYPE_PATTERN}/""
READTHROUGH_PATTERN="tag /"readthrough_transcript/""

# 生成基因允许列表
cat "$gtf_modified" /
    | awk '$3 == "transcript"' /  # 仅保留转录本条目
    | grep -E "$GENE_PATTERN" /  # 筛选符合类型的基因
    | grep -E "$TX_PATTERN" /  # 筛选符合类型的转录本
    | grep -Ev "$READTHROUGH_PATTERN" /  # 排除通读转录本
    | sed -E 's/.*(gene_id "[^"]+").*//1/' /
    | sort | uniq /
    > "${build}/gene_allowlist"

# 基于允许列表过滤 GTF，并排除 Y 染色体 PAR 区域冗余基因
gtf_filtered="${build}/$(basename "$gtf_in").filtered"
grep -E "^#" "$gtf_modified" > "$gtf_filtered"  # 保留注释头
grep -Ff "${build}/gene_allowlist" "$gtf_modified" /
    | awk -F "/t" '$1 != "chrY" || $1 == "chrY" && $4 >= 2752083 && $4 < 56887903 && !/ENSG00000290840/' /
    >> "$gtf_filtered"

四、整合人类与 EBV 基因组及注释

将人类和 EBV 的基因组序列、基因注释文件分别合并，形成融合文件：

# 修正 EBV 基因组头部格式
sed -E "s/>NC_007605.1 Human gammaherpesvirus 4, complete genome/>NC_007605/" reference_sources/GCF_002402265.1_ASM240226v1_genomic.fna > $build/GCF_002402265.1_ASM240226v1_genomic.modified.fna

# 优化 EBV 注释文件，因为一个基因有多个外显子，而10x要求gene_id唯一
awk -f gtf.awk reference_sources/GCF_002402265.1_ASM240226v1_genomic.gtf > ${build}/GCF_002402265.1_ASM240226v1_genomic.gtf.fix

# 合并人类与 EBV 的 GTF 和 FASTA
cat  ${build}/GCF_002402265.1_ASM240226v1_genomic.gtf.fix $gtf_filtered > combined.gtf
cat reference_sources/GCF_002402265.1_ASM240226v1_genomic.fna $fasta_modified > combined.fna

gtf.awk

{
    # 匹配 gene_id 的模式，捕获基因ID内容
    if (match($0, /gene_id "([^"]+)"/, arr)) {
        gene_id = arr[1]
        # 检查 gene_id 是否在数组中，若在则计数加1，不在则初始化为1
        if (gene_id in gene_count) {
            gene_count[gene_id]++
        } else {
            gene_count[gene_id] = 1
        }
        # 根据计数情况修改 gene_id
        if (gene_count[gene_id] > 1) {
            gsub(/gene_id "[^"]+"/, "gene_id /"" gene_id "_" gene_count[gene_id] "/"")
        }
    }
    print $0
}

五、使用 Cell Ranger 构建最终参考文件

通过 Cell Ranger 的 mkref 工具将融合文件打包为 10X 兼容的参考基因组：

mkdir -p $version

# 构建参考文件
cellranger mkref /
    --ref-version="$version" /
    --genome="$genome" /
    --fasta=combined.fna /
    --genes=combined.gtf /
    --nthreads=30 /  # 多线程加速
    --output-dir $version

结果与应用

成功构建的人+EBV 融合参考基因组（$version 文件夹）可直接用于 10X 单细胞转录组数据的定量分析（通过 cellranger count 命令）。相比传统人类参考基因组，该融合参考能：

1. 精准定量 EBV 基因（如 EBNA、LMP 家族）在单细胞中的表达水平；
2. 捕捉病毒感染细胞的特异性基因表达特征；
3. 解析病毒-宿主相互作用的分子机制，为疾病诊断和治疗提供新靶点。