#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 123A10132-D chr1 115256501 . T C . orientation;strand_bias AS_FilterStatus=strand_bias;AS_SB_TABLE=6559,6567|1,2078;DP=15208;ECNT=1;GERMQ=93;MBQ=20,20;MFRL=121,119;MMQ=60,60;MPOS=15;POPAF=7.30;ROQ=1;TLOD=3052.41 GT:AD:AF:DP:F1R2:F2R1:FAD:SB 0/1:13126,2079:0.237:15205:0,0:6564,2076:6579,2078:6559,6567,1,2078

<think>我们已经有32个VCF文件，每个文件对应一个病人的肿瘤-正常配对样本的变异调用结果。现在需要对这些VCF文件进行后续分析，包括合并、过滤、注释、可视化和生物学解释。由于用户是Windows系统的小白，我们将尽量使用图形界面工具或简单命令行工具，并提供详细步骤。 步骤概述： 1. 合并VCF文件（可选，但便于后续分析） 2. 筛选体细胞突变（肿瘤样本中突变而正常样本中无突变或低频率） 3. 变异注释（使用VEP） 4. 过滤（根据质量分数、深度等） 5. 可视化（使用IGV和maftools） 6. 生物学意义挖掘（突变频谱、驱动基因等） 详细步骤： 一、合并VCF文件（使用bcftools） 虽然我们有32个VCF文件，但为了后续分析方便，我们可以将它们合并成一个VCF文件。注意：合并前需要确保每个VCF文件都已经压缩并建立索引（.vcf.gz和.tbi文件）。如果还没有，请先压缩和索引。 1. 安装bcftools（Windows版）并添加到系统路径。 下载地址：https://github.com/samtools/bcftools/releases 选择Windows版本（如bcftools-1.16-win64.zip），解压后将bcftools.exe所在的目录添加到系统环境变量PATH中。 2. 压缩和索引每个VCF文件（如果还没有压缩）： 在存放VCF文件的文件夹中，打开命令提示符（cmd），执行以下命令（逐文件处理）： ```bash for %v in (*.vcf) do bgzip %v for %v in (*.vcf.gz) do bcftools index %v ``` 注意：如果使用PowerShell，循环语法不同。也可以在cmd中逐一手工执行。 3. 合并VCF文件： ```bash bcftools merge -m all -O z -o merged.vcf.gz *.vcf.gz ``` 参数说明： -m all: 合并所有记录 -O z: 输出为压缩的VCF（.vcf.gz） -o: 输出文件名 4. 为合并后的文件建立索引： ```bash bcftools index merged.vcf.gz ``` 二、筛选体细胞突变 在合并的VCF文件中，每个样本有两个列（肿瘤和正常）。我们需要筛选出在肿瘤样本中突变而在正常样本中为野生型（或极低频率）的位点。 使用bcftools过滤： ```bash bcftools filter -i 'FORMAT/GT[0]="0/1" || FORMAT/GT[0]="1/1"' merged.vcf.gz | bcftools filter -i 'FORMAT/GT[1]="0/0"' -O z -o somatic.vcf.gz ``` 解释： - 第一个bcftools filter：选择肿瘤样本（假设每个VCF中肿瘤是第一列，正常是第二列）基因型为杂合（0/1）或纯合突变（1/1）的位点。 - 第二个bcftools filter：从上述结果中进一步选择正常样本为野生型（0/0）的位点。 - 注意：如果样本顺序不一致，需要先查看合并后的VCF中的样本顺序（用`bcftools query -l merged.vcf.gz`），然后调整索引[0]和[1]的位置。 另一种方法是使用Mutect2的FILTER列，通常体细胞突变会标记为“PASS”。我们可以直接提取PASS的位点： ```bash bcftools view -f PASS merged.vcf.gz -O z -o somatic_pass.vcf.gz ``` 三、变异注释（使用VEP） 1. 安装VEP（Windows下可以通过WSL或直接使用perl脚本，但这里我们使用Windows版VEP，需要安装Strawberry Perl和VEP） 参考：https://www.ensembl.org/info/docs/tools/vep/script/vep_download.html#windows 步骤： a) 安装Strawberry Perl: http://strawberryperl.com/ b) 下载VEP安装包：https://github.com/Ensembl/ensembl-vep/archive/refs/heads/main.zip c) 解压到某个目录，例如D:\vep d) 打开命令提示符，进入vep目录，运行： ```bash perl INSTALL.pl ``` 在安装过程中，会提示下载缓存文件（cache），选择人类基因组版本（如GRCh38）和相应的缓存文件。 2. 注释somatic_pass.vcf.gz文件： ```bash perl vep.pl -i somatic_pass.vcf.gz --format vcf --output_file somatic_annotated.vcf --cache --offline --assembly GRCh38 --everything --force_overwrite ``` 参数说明： -i: 输入文件 --format: 输入格式 --output_file: 输出文件 --cache: 使用缓存 --offline: 离线模式 --assembly: 基因组版本 --everything: 开启所有注释选项 --force_overwrite: 覆盖已有输出文件 3. 输出文件为somatic_annotated.vcf，其中包含注释信息（在INFO字段中）。 四、进一步过滤 在注释后，我们可能希望根据注释信息进行过滤，例如保留错义突变、无义突变、剪接位点突变等。我们可以使用bcftools进行过滤。 例如，只保留具有“HIGH”或“MODERATE”影响的变异（根据VEP注释的CSQ字段中的IMPACT）： ```bash bcftools filter -i 'INFO/CSQ ~ "HIGH" || INFO/CSQ ~ "MODERATE"' somatic_annotated.vcf -O z -o somatic_annotated_filtered.vcf.gz ``` 注意：VEP注释信息在CSQ字段中，但上述过滤条件可能不够精确，因为CSQ字段包含多个转录本信息。更精确的方法是使用bcftools的+split-vep插件将CSQ字段拆分成多列，然后过滤。 步骤： 1. 拆分CSQ字段： ```bash bcftools +split-vep somatic_annotated.vcf.gz -O z -o somatic_annotated_split.vcf.gz ``` 2. 然后根据IMPACT字段过滤： ```bash bcftools filter -i 'INFO/CSQ_IMPACT="HIGH" || INFO/CSQ_IMPACT="MODERATE"' somatic_annotated_split.vcf.gz -O z -o somatic_annotated_filtered.vcf.gz ``` 五、可视化 1. 使用IGV查看变异： 下载并安装IGV：https://software.broadinstitute.org/software/igv/ 加载参考基因组（GRCh38），然后载入somatic_annotated_filtered.vcf.gz文件。可以查看特定基因的突变情况。 2. 使用maftools（在R中）进行汇总和可视化： 安装R和RStudio（参考之前的步骤）。 在RStudio中执行以下代码： ```r # 安装maftools if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("maftools") # 加载maftools library(maftools) # 将VCF转换为MAF格式（maftools需要MAF格式） # 注意：maftools可以直接读取VCF，但需要是未压缩的VCF。所以我们先解压，或者使用vcf2maf工具转换。 # 这里我们使用maftools的vcf2maf函数（注意：这个函数需要Perl和VEP，在Windows下可能复杂，所以我们选择直接读取VCF并转换） # 另一种方法：使用read.vcfR读取VCF，然后转换为MAF。 # 安装vcfR包 install.packages("vcfR") library(vcfR) # 读取VCF vcf <- read.vcfR("somatic_annotated_filtered.vcf.gz") # 转换为maftools需要的MAF格式 maf <- vcfR2maf(vcf, assembly="hg38") # 将MAF数据框写入文件 write.table(maf, file="somatic_annotated_filtered.maf", sep="\t", quote=FALSE, row.names=FALSE) # 使用maftools读取MAF somatic_maf <- read.maf(maf = "somatic_annotated_filtered.maf") # 绘制摘要图 plotmafSummary(maf = somatic_maf, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE) # 绘制oncoplot（瀑布图） oncoplot(maf = somatic_maf, top = 20) # 其他分析，如突变频谱 somatic_maf.titv = titv(maf = somatic_maf, plot = FALSE, useSyn = TRUE) plotTiTv(res = somatic_maf.titv) ``` 六、生物学意义挖掘 1. 驱动基因识别：通过Oncoplot图观察高频突变基因，或使用maftools中的函数进行驱动基因分析（如OncodriveCLUST）。 2. 通路富集：使用clusterProfiler包对突变基因进行KEGG或GO富集分析。 3. 突变特征分析：可以使用maftools的trinucleotideMatrix函数提取三核苷酸上下文，然后使用deconstructSigs包进行突变特征分析（需要安装R包deconstructSigs）。 示例代码（在R中）： ```r # 突变特征分析 if (!require("deconstructSigs", quietly = TRUE)) install.packages("deconstructSigs") library(deconstructSigs) # 从maf对象中提取三核苷酸上下文 somatic_maf.sig = trinucleotideMatrix(maf = somatic_maf, prefix = 'chr', add = TRUE, ref_genome = "BSgenome.Hsapiens.UCSC.hg38") # 注意：需要安装BSgenome.Hsapiens.UCSC.hg38 if (!require("BSgenome.Hsapiens.UCSC.hg38", quietly = TRUE)) BiocManager::install("BSgenome.Hsapiens.UCSC.hg38") # 重新运行trinucleotideMatrix somatic_maf.sig = trinucleotideMatrix(maf = somatic_maf, prefix = 'chr', add = TRUE, ref_genome = "BSgenome.Hsapiens.UCSC.hg38") # 转换为deconstructSigs输入格式 sigs.input = somatic_maf.sig$nmf_matrix # 将每个样本的突变频谱转换为96维向量 # 然后使用deconstructSigs进行突变特征分解 signatures = whichSignatures(tumor.ref = sigs.input, signatures.ref = signatures.cosmic, sample.id = "Sample1", contexts.needed = TRUE) # 绘制样本的突变特征组成图 plotSignatures(signatures) ``` 注意：由于样本量较小（32个病人），驱动基因分析和突变特征分析可能不够显著，但可以尝试。 七、其他注意事项 1. 在Windows下使用命令行工具，注意路径中不要有空格，最好使用英文目录。 2. 每一步操作都要检查输出文件是否正确生成。 3. 如果遇到内存不足，可以考虑分步骤处理或减少同时处理的样本数。 相关问题： 1. 如何确保在合并VCF文件时样本顺序正确？ 2. 在VEP注释中，如何只保留特定基因的变异？ 3. 如何用maftools比较不同种族间的突变差异？ 4. 如何从VCF文件中提取特定基因（如TP53）的突变？ 希望这些步骤能帮助你完成分析。</think>### 针对32名病人肿瘤/癌旁配对样本的VCF数据分析流程（Windows系统版） #### 一、准备工作 1. **软件安装**（全部支持Windows）： | 工具 | 作用 | 下载地址 | |------|------|----------| | **BCFtools** | VCF文件处理 | [https://github.com/samtools/bcftools/releases](https://github.com/samtools/bcftools/releases) | | **VEP** | 变异注释 | [https://www.ensembl.org/info/docs/tools/vep/script/vep_download.html](https://www.ensembl.org/info/docs/tools/vep/script/vep_download.html) | | **R + RStudio** | 统计分析与可视化 | [https://cran.r-project.org](https://cran.r-project.org) + [https://posit.co/download/rstudio-desktop](https://posit.co/download/rstudio-desktop) | | **IGV** | 突变可视化 | [https://software.broadinstitute.org/software/igv](https://software.broadinstitute.org/software/igv) | 2. **文件组织**： ```bash D:\cancer_analysis ├── input_vcfs # 存放32个VCF文件 ├── processed_data # 处理后的文件 ├── scripts # 分析脚本 └── results # 最终结果 ``` 3. **参考数据下载**： - 人类基因组参考：`GRCh38.fa`（[ENSEMBL下载](https://asia.ensembl.org/info/data/ftp/index.html)） - VEP注释数据库：`homo_sapiens_vep_104_GRCh38.tar.gz`（解压到`D:\vep_data`） --- #### 二、核心分析流程 ##### 步骤1：VCF文件预处理（BCFtools） ```bash # 1. 压缩和索引所有VCF文件 for %f in (input_vcfs\*.vcf) do ( bgzip -c "%f" > "processed_data\%~nf.vcf.gz" bcftools index "processed_data\%~nf.vcf.gz" ) # 2. 合并所有VCF文件 bcftools merge processed_data\*.vcf.gz -O z -o processed_data\merged.vcf.gz bcftools index processed_data\merged.vcf.gz ``` ##### 步骤2：筛选体细胞突变（肿瘤特异性） ```bash # 1. 提取肿瘤样本突变（样本名后缀_T） bcftools view -s $(bcftools query -l merged.vcf.gz | findstr "_T") processed_data\merged.vcf.gz > processed_data\tumor_only.vcf # 2. 筛选肿瘤突变/癌旁野生型 bcftools filter -i 'FORMAT/GT[0]="alt" && FORMAT/GT[1]="ref"' processed_data\tumor_only.vcf > processed_data\somatic_candidates.vcf # 3. 质量过滤（保留PASS标记） bcftools filter -i 'FILTER="PASS"' processed_data\somatic_candidates.vcf > processed_data\somatic_filtered.vcf ``` ##### 步骤3：变异注释（VEP） 1. 创建配置文件`D:\vep\vep.ini`： ```ini [DEFAULT] dir = D:\vep_data cache = yes assembly = GRCh38 [data] FASTA = D:\vep_data\GRCh38.fa ``` 2. 运行注释： ```bash perl vep.pl -i processed_data\somatic_filtered.vcf --format vcf -o results\somatic_annotated.vcf --cache --force_overwrite --everything ``` 输出文件包含关键注释信息： - `Consequence`: 变异类型（如错义突变） - `IMPACT`: 功能影响（HIGH/MODERATE） - `SYMBOL`: 基因符号 - `HGVSp`: 蛋白质变化 ##### 步骤4：R语言分析（RStudio） 1. **安装必要包**： ```r install.packages("tidyverse") BiocManager::install(c("maftools", "VariantAnnotation")) ``` 2. **突变频谱分析**： ```r library(maftools) library(ggplot2) # 导入VEP注释结果 maf_data = read.maf("results/somatic_annotated.vcf") # 突变频谱图 plotmafSummary(maf = maf_data, rmOutlier = TRUE) # 前20突变基因瀑布图 oncoplot(maf = maf_data, top = 20, fontSize = 10) # 突变类型分布 titv_plot = titv(maf = maf_data, plot = FALSE) plotTiTv(res = titv_plot) ``` 3. **驱动基因识别**： ```r # 识别显著突变基因 sig_genes = oncodrive(maf = maf_data, minMut = 5, pvalMethod = "zscore") head(sig_genes$results, 10) # 绘制驱动基因热图 somaticInteractions(maf = maf_data, top = 25, pvalue = 0.05) ``` ##### 步骤5：突变可视化（IGV） 1. 打开IGV → 加载基因组`GRCh38` 2. 导入文件`results/somatic_annotated.vcf` 3. 操作示例： - 输入基因名（如`TP53`）定位 - 右键突变点 → 查看注释详情 - 比较肿瘤/癌旁样本的覆盖深度 --- #### 三、生物学意义挖掘 1. **通路富集分析**（R）： ```r library(clusterProfiler) library(org.Hs.eg.db) # 提取突变基因列表 gene_list = [email protected]$Hugo_Symbol # KEGG通路富集 kegg_result = enrichKEGG(gene = gene_list, organism = 'hsa', pvalueCutoff = 0.05) dotplot(kegg_result, showCategory=15) # GO生物过程富集 go_result = enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, ont = "BP", pvalueCutoff = 0.01) barplot(go_result, showCategory=10) ``` 2. **突变特征分析**： ```r # 提取三核苷酸上下文 sig_matrix = trinucleotideMatrix(maf = maf_data, ref_genome = "BSgenome.Hsapiens.UCSC.hg38") # 突变特征分解（需安装deconstructSigs） library(deconstructSigs) signatures = whichSignatures(tumor.ref = sig_matrix$nmf_matrix, signatures.ref = signatures.cosmic, sample.id = "all") plotSignatures(signatures) ``` --- #### 四、结果解读要点 1. **关键文件**： - `somatic_annotated.vcf`：包含基因注释的体细胞突变 - `oncoplot.pdf`：展示高频突变基因的瀑布图 - `kegg_result.csv`：显著富集的癌症通路 2. **重点指标**： - **高频突变基因**（出现在>10%样本中） - **有害突变类型**： - `stop_gained`（无义突变） - `frameshift_variant`（移码突变） - `splice_acceptor_variant`（剪接受体突变） - **驱动基因**（oncodrive分析中$p<0.01$） 3. **临床关联建议**： - 关注已知癌症基因（如TP53、KRAS、PIK3CA） - 检查治疗靶点相关基因（如EGFR、BRCA1/2） - 分析肿瘤突变负荷（TMB）与免疫治疗相关性 --- #### 五、常见问题解决 1. **BCFtools报错**： - 确保VCF文件用`bgzip`压缩（非普通zip） - 检查样本名一致性：`bcftools query -l merged.vcf.gz` 2. **VEP注释失败**： ```bash # 重建FASTA索引 samtools faidx D:\vep_data\GRCh38.fa ``` 3. **R包安装失败**： - 在RStudio中运行：`options(download.file.method = "wininet")` - 尝试国内镜像：`options(repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")` 4. **内存不足处理**： ```bash # 分染色体处理大文件 bcftools view merged.vcf.gz -r chr1 -O z -o chr1.vcf.gz ```