首页归档Schema进度实践关于

CNV 检测

临床遗传全外显子一般需要分析"外显子"和"大片段"两个水平。按照以往逻辑，会使用XHMM、ExomeDepth等工具，但我个人是All In CNVkit的。

基线准备

从既往经验来说，使用同批次产生的同性别数据作为基线样本，会有较好的表现（一般需要≥8个样本）。

但是显然，这对于样本量少的公司，无法实现。而且，这样做对于后台流程也会比较复杂。因此，我们还是采取固定基线的方法，采取若干的正常表型数据及格的样本进行CNV基线建立。建议男女性样本各20个。

当然，由于这次我们使用的是公共数据，因此只能虚空演示。

首先是窗口划分，与WGS、CNV-seq等划分固定大小的窗口（bins）的方式不同，WES建议以外显子进行窗口划分，即每个bins实际是不同大小的。这样做的好处是最终我们可以去评估单个外显子水平的CNV；后面也可以对窗口进行segment，从而获得大片段水平的CNV。

显然，我们需要先去获得一个外显子的bed。

外显子bed

这里基于Gencode来制作一个GRCh38的外显子bed。首先下载gff文件（GRCh37可以通过Gencode下载）

wget https://ftp.ncbi.nlm.nih.gov/refseq/MANE/MANE_human/release_1.3/MANE.GRCh38.v1.3.refseq_genomic.gff.gz
gunzip MANE.GRCh38.v1.3.refseq_genomic.gff.gz

然后处理为bed格式

output_bed = open("GRCh38.NCBI.ManeSelect.Exon.bed", "w", encoding="utf-8")
output_bed.write("#Chrom\tStart\tEnd\tRegion\tGene\tEnsembl\tRefseq\tStrand\n")
with open("MANE.GRCh38.v1.3.refseq_genomic.gff", "r") as f:
    for line in f:
        if line.startswith("#"):
            continue
        ensembl = refseq = "."
        exon_id = "exon1"
        gene = ""
        lines = line.rstrip().split("\t")
        if (lines[2] == "exon") and ("tag=MANE Select" in line):
            info_fields = lines[8].split(";")
            for i in info_fields:
                if i.startswith("Dbxref"):
                    dbxrefs = i.lstrip("Dbxref=").split(",")
                    for d in dbxrefs:
                        if d.startswith("Ensembl:"):
                            ensembl = d.split("Ensembl:")[1]
                elif i.startswith("transcript_id="):
                    refseq = i.split("transcript_id=")[1]
                elif i.startswith("ID=exon-"):
                    exon_id = "exon" + i.split("-")[-1]
                elif i.startswith("gene="):
                    gene = i.split("gene=")[1]

            output_list = [lines[0], lines[3], lines[4], exon_id, gene, ensembl, refseq, lines[6]]
            output_bed.write("\t".join(output_list) + "\n")
output_bed.close()

最后，建议使用bedtools进行一次merge，用以去掉一些同源重复区域。

bedtools merge -i GRCh38.NCBI.ManeSelect.Exon.bed > GRCh38.NCBI.ManeSelect.Exon.merge.bed

CNVKit分析流程

CNVKit的分析分为几个步骤：首先生成antitarget区域，然后分别计算target和antitarget的覆盖度，最后合并生成结果。

生成antitarget区域

antitarget是靶向区域之间的区域，用于估计背景噪声。需要排除不可映射区域：

# 获取可访问区域，排除不可映射区域
cnvkit.py access GRCh38.d1.vd1.fa -x hg38_excludeble.bed -o access.bed

# 生成antitarget区域
cnvkit.py antitarget target.bed -g access.bed -o SRR14724513.antitarget.bed

计算覆盖度

分别计算target和antitarget区域的覆盖度，使用 -p 参数指定并行线程数：

# target覆盖度
cnvkit.py coverage bam/SRR14724513.markdup.bam target.bed \
    -o SRR14724513.targetcoverage.cnn -p 16

# antitarget覆盖度
cnvkit.py coverage bam/SRR14724513.markdup.bam SRR14724513.antitarget.bed \
    -o SRR14724513.antitargetcoverage.cnn -p 16

合并结果

使用预先构建好的参考基线，将target和antitarget的覆盖度合并：

cnvkit.py fix \
    SRR14724513.targetcoverage.cnn \
    SRR14724513.antitargetcoverage.cnn \
    cnvkit_reference.cnn \
    -o SRR14724513.cnvkit.cnr

.cnr文件包含了每个bin的log2 ratio值，是后续CNV判断的基础。

外显子水平

CBS基因水平分割

在得到CNVKit的.cnr文件后，使用CBS（Circular Binary Segmentation）算法在基因水平上进行分割，评估每个基因的拷贝数状态：

python3 gene_level_segment.py \
    -i SRR14724513.cnvkit.cnr \
    -o SRR14724513.cnv.gene.txt \
    --cn-del-threshold 1.5 \
    --cn-dup-threshold 2.5

CBS算法会将连续的bins按照log2 ratio的显著变化点进行分段，然后对每段计算平均拷贝数。--cn-del-threshold和--cn-dup-threshold用于定义拷贝数变异的阈值：

  • CN < 1.5：缺失（Deletion）
  • CN > 2.5：重复（Duplication）
  • 1.5 ≤ CN ≤ 2.5：正常

报告会包含每个基因的拷贝数估计、置信度评分、变异类型等信息。

大片段水平

窗口合并

为了获得大片段水平的CNV，需要将小的外显子bins合并为固定大小的窗口（如20kb），然后在窗口水平上进行评估：

python3 merge_bins_for_cnv.py \
    -i SRR14724513.cnvkit.cnr \
    -o SRR14724513.cnv.region.txt \
    --del-threshold 1.5 \
    --dup-threshold 2.5 \
    --bin-size 20000 \
    --keep-antitarget

--bin-size 20000指定合并后的窗口大小为20kb。--keep-antitarget保留antitarget区域的信息用于辅助判断。合并时使用加权平均，权重为每个bin的probe数量。

常见的窗口大小选择：

  • 10kb：较细粒度，适合检测较小的CNV
  • 20kb：平衡粒度，推荐默认值
  • 50kb：较粗粒度，适合检测大片段CNV

CNV注释

使用CNVAnno对检出的CNV进行基因注释：

python3 cnvanno.py \
    SRR14724513.cnv.gene.txt \
    -d /app/CNVAnno/data \
    -g GRCh38 \
    -f tsv \
    -o SRR14724513.cnvanno.txt

CNVAnno会将CNV区域与基因数据库进行比对，注释受影响的基因及其关联的疾病信息。也可以使用在线工具AutoCNV进行辅助分析。

实际上CNV的ClinGen规则判定受主观因素影响，建议有能力的可以自行根据ClinGen规则搭建适用于自己的注释流程。