首页归档Schema进度实践关于

数据的准备

此篇是遗传全外显子数据分析流程，适用于临床遗传全外显子数据分析。抛开历史包袱，使用在我个人认知中的更新的技术栈。

原始数据

这里使用SRR14724513，这是一个NA12878的数据，使用的是安捷伦的SureSelect v7探针，测序深度大概是50X。

使用这个数据是因为方便我们下游进行一个NA12878 Benchmark，以及他表明了使用的探针，可以比较容易找到对应的bed文件。

下载数据，并使用sratools将数据拆分为fastq。

wget https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR14724513/SRR14724513

使用sratoolkit中的fastq-dump拆分数据

fastq-dump --split-3 --gzip ./SRR14724513

下载bed文件

二代测序中WES的湿实验一般是探针捕获，我们需要去计算探针的捕获效率，也需要知道我们感兴趣的靶向区域的覆盖度、深度等质控信息，因此必须要有探针所在的区域信息，因此需要去获取bed文件。尽管都是"全外显子"，但每个厂家之间的探针还是存在差异的，比如，会对一些难以检测的区域进行特殊设计，如用侧翼探针进行覆盖，又或者是探针密度增加等。

安捷伦现在主推的V8+NCV和CREV4还会把一些临床关注的内含子加入到区域中。

对于上面的数据，我们下载安捷伦的V7 bed。安捷伦的探针可以在他的线上系统找到，但由于这个系统需要注册，这里在github里找到一份。

这份教程基于hg38，因此我下载了hg38版本。

现在，我们有原始数据和bed文件了。

这个bed文件中包含了track抬头，建议把抬头注释掉，或者删除。

参考数据库包

流程所需的参考数据统一打包，可以从HuggingFace下载（pzweuj/SchemaBio_Bundle），包含以下内容：

  • 参考基因组FASTA：GRCh37或GRCh38
  • BWA-MEM2索引：.bwt.2bit.64等索引文件
  • VEP Cache：合并版Cache（包含Ensembl和RefSeq转录本）
  • Schema Bundle：包含ClinVar VCF、cytoBand BED、GnomAD v4.1 filtered VCF、Pangolin VCF、EVOScore2 VCF、AlphaMissense v3 VCF
# 示例目录结构
ref_dir/
├── GRCh38.d1.vd1.fa
├── GRCh38.d1.vd1.fa.bwt.2bit.64
├── GRCh38.d1.vd1.fa.fai
├── GRCh38.d1.vd1.fa.pac
└── ...

cache_dir/
└── homo_sapiens_merged_vep_113_GRCh38/

schema_bundle/
├── hg38_clinvar_20260415.vcf.gz
├── hg38_cytoBand.bed.gz
├── hg38_gnomad.v4.1.filtered.vcf.gz
├── hg38_pangolin.vcf.gz
├── hg38_EVOScore2.vcf.gz
└── hg38_AlphaMissense.v3.vcf.gz

BED文件处理

不同的工具和分析步骤需要不同的bed文件。原始的探针bed文件可能包含"chr"前缀，或者线粒体标注为"chrM"而非"MT"，需要进行标准化处理。

去除chr前缀

在GRCh38中，有些bed文件使用"chr1"格式，有些使用"1"格式，需要统一。同时将"chrM"转换为"MT"：

# 去除chr前缀，chrM转为MT
sed 's/^chr//' input.bed | sed 's/^MT$/MT/' > fixed.bed

也可以使用脚本批量处理，自动修复零长度区域等问题。

外显子bed制作

除了探针bed，我们还需要一个外显子水平的bed，用于CNV分析的窗口划分。这里基于NCBI MANE Select注释来制作：

# 下载MANE GFF
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

从GFF中提取外显子区域并转为bed格式，然后使用bedtools merge去除重叠区域：

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

CNVKit靶向区域bed

CNV分析还需要一个与探针区域取交集后的靶向bed，用于区分target和antitarget：

bedtools intersect -b probe.bed -a Gencode.GRCh38.cnvkit.target.bed -wa -u > target.bed

样本指纹

样本指纹用于追踪样本身份，防止样本混淆。使用Pengelly等人发表的24-SNP身份鉴定面板，通过bcftools对这24个位点进行基因分型：

bcftools mpileup \
    -f GRCh38.d1.vd1.fa \
    -R pengelly_snp.txt \
    --annotate FORMAT/AD,FORMAT/DP \
    bam/SRR14724513.markdup.bam | \
bcftools call -m -Ov -o fingerprint.vcf

这个指纹信息会在后续的QC汇总报告中使用，也可以跨样本来追踪身份一致性。

原始数据质控

测序过程中，由于种种影响，原始数据里可能会存在低质量的接头序列、碱基、低质量reads等，一般需要先对数据进行一次过滤，获得高质量的数据。

这里使用fastp进行过滤。线程数上限为16。

fastp \
    -i rawdata/SRR14724513_1.fastq.gz \
    -I rawdata/SRR14724513_2.fastq.gz \
    -o cleandata/SRR14724513_R1.clean.fq.gz \
    -O cleandata/SRR14724513_R2.clean.fq.gz \
    -w 16 \
    -j qc/SRR14724513.fastp_stats.json \
    -h qc/SRR14724513.fastp_stats.html \
    --detect_adapter_for_pe

fastp会自动检测接头序列，生成JSON和HTML格式的质控报告。