全外显子组测序（Whole Exome Sequencing, WES）技术简介

定义与原理

全外显子组测序（WES）是一种高通量测序技术，专注于对人类基因组中外显子区域（蛋白质编码区域，约占基因组的1-2%）进行深度测序。

外显子：基因中编码蛋白质的DNA片段，包含大部分已知的致病突变。

技术原理：通过探针杂交或PCR富集外显子区域，再结合高通量测序技术（如Illumina）进行测序。

Note

这个实战中大部分参考了https://bohrium.dp.tech/notebooks/2112719645中的教程，但是呢，我有以下几个改动：
1.原教程中进行种系突变分析的时候有些错误（对 tumor 组织调用了HaplotypeCaller），我在我的实战教程中做了更正，并把原教程中使用的 germline 进行了作图。
2.原教程中只对一个 normal 和一个 tumor 进行了体细胞突变的分析，在我的教程中将case1 中的其他的生物学重复以及技术重复都包括在内。

技术特点

优势局限性
高性价比：仅测序1-2%基因组，成本低覆盖不全：无法检测非编码区变异
高深度（通常>100×），适合低频突变检测捕获效率差异：部分外显子可能覆盖不足
临床适用性强：已知致病突变多位于外显子结构变异检测能力有限（如大片段缺失/重复）

Data description

case1和case2 都来源与乳腺癌病人，germline 是 normal 组织，生物学重复是 tumor 组织。
case 1：Germline（SRR3182423）、Biological replicate A （SRR3182433）
case 2：Germline（SRR3182419）、Biological replicate A （SRR3182436）

Download exon data and reference genome

下载参考基因组并构建索引：

wget -c ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.gz -P ./0.ref
gzip --decompress ./0.ref/human_g1k_v37.fasta.gz

bwa index -p ./human_g1k_v37 ../0.Ref_Genome/human_g1k_v37.fasta

通过 UCSC下载外显子数据

Table Browser

之后点击 Get output

Get output

点击 get BED

创建一个文件：exon_download_24chr.txt，来筛选出exon data中的染色体

chr1
chr10
chr11
chr12
chr13
chr14
chr15
chr16
chr17
chr18
chr19
chr2
chr20
chr21
chr22
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chrX
chrY

grep -Fwf exon_download_24chr.txt exon_download.bed > filtered_exon_download.bed

把 chr去掉，比如 chr1 改成 1。

awk '{gsub("chr", "", $0)}1' filtered_exon_download.bed > exon_data.bed

Download Rawdata

其中~/.aspera/connect/etc/asperaweb_id_dsa.openssh，需要写入自己的openssh文件的位置。使用ascp下载速度比 fast dump 快很多。

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn

# Case 1: Germline (SRR3182423) 和 Biological replicate A (SRR3182433)
ascp -k 1 -QT -l 500m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR318/003/SRR3182423/SRR3182423_1.fastq.gz ./SRR3182423_1.fastq.gz

ascp -k 1 -QT -l 500m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR318/003/SRR3182423/SRR3182423_2.fastq.gz ./SRR3182423_2.fastq.gz

ascp -k 1 -QT -l 500m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR318/003/SRR3182433/SRR3182433_1.fastq.gz ./SRR3182433_1.fastq.gz

ascp -k 1 -QT -l 500m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR318/003/SRR3182433/SRR3182433_2.fastq.gz ./SRR3182433_2.fastq.gz

# Case 2: Germline (SRR3182419) 和 Biological replicate A (SRR3182436)
ascp -k 1 -QT -l 500m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR318/009/SRR3182419/SRR3182419_1.fastq.gz ./SRR3182419_1.fastq.gz

ascp -k 1 -QT -l 500m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR318/009/SRR3182419/SRR3182419_2.fastq.gz ./SRR3182419_2.fastq.gz

ascp -k 1 -QT -l 500m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR318/006/SRR3182436/SRR3182436_1.fastq.gz ./SRR3182436_1.fastq.gz

ascp -k 1 -QT -l 500m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR318/006/SRR3182436/SRR3182436_2.fastq.gz ./SRR3182436_2.fastq.gz

Trim

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=25G
#SBATCH --cpus-per-task=5
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn

conda activate trim_galore

# Case 1: Germline (SRR3182423) 和 Biological replicate A (SRR3182433)
trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 1 \
-o ./ ./SRR3182423_1.fastq.gz ./SRR3182423_2.fastq.gz

trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 1 \
-o ./ ./SRR3182433_1.fastq.gz ./SRR3182433_2.fastq.gz

# Case 2: Germline (SRR3182419) 和 Biological replicate A (SRR3182436)
trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 1 \
-o ./ ./SRR3182419_1.fastq.gz ./SRR3182419_2.fastq.gz

trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 1 \
-o ./ ./SRR3182436_1.fastq.gz ./SRR3182436_2.fastq.gz

Mapping

Mapping之后的结果是 sam，再将 sam 转换成 bam 文件

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn


bwa mem -M -t 5 -R "@RG\tID:SRR3182419\tSM:SRR3182419\tLB:WXS\tPL:Illummina" ../4.bwa_index/human_g1k_v37 \
    ../3.Trimming/SRR3182419_1_val_1.fq.gz \
    ../3.Trimming/SRR3182419_2_val_2.fq.gz \
    > ./SRR3182419.sam

bwa mem -M -t 5 -R "@RG\tID:SRR3182423\tSM:SRR3182423\tLB:WXS\tPL:Illummina" ../4.bwa_index/human_g1k_v37 \
    ../3.Trimming/SRR3182423_1_val_1.fq.gz \
    ../3.Trimming/SRR3182423_2_val_2.fq.gz \
    -o ./SRR3182423.sam

bwa mem -M -t 5 -R "@RG\tID:SRR3182433\tSM:SRR3182433\tLB:WXS\tPL:Illummina" ../4.bwa_index/human_g1k_v37 \
    ../3.Trimming/SRR3182433_1_val_1.fq.gz \
    ../3.Trimming/SRR3182433_2_val_2.fq.gz \
    -o ./SRR3182433.sam

bwa mem -M -t 5 -R "@RG\tID:SRR3182436\tSM:SRR3182436\tLB:WXS\tPL:Illummina" ../4.bwa_index/human_g1k_v37 \
    ../3.Trimming/SRR3182436_1_val_1.fq.gz \
    ../3.Trimming/SRR3182436_2_val_2.fq.gz \
    -o ./SRR3182436.sam

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o convert.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn

samtools sort -@ 5 ./SRR3182419.sam -o ./SRR3182419.bam
samtools sort -@ 5 ./SRR3182423.sam -o ./SRR3182423.bam
samtools sort -@ 5 ./SRR3182433.sam -o ./SRR3182433.bam
samtools sort -@ 5 ./SRR3182436.sam -o ./SRR3182436.bam

对 Mapping 进行质控

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn


samtools stats --reference ../0.Ref_Genome/human_g1k_v37.fasta ../5.Mapping/SRR3182419.bam > ./SRR3182419.stat
plot-bamstats -p ./SRR3182419 ./SRR3182419.stat

samtools stats --reference ../0.Ref_Genome/human_g1k_v37.fasta ../5.Mapping/SRR3182423.bam > ./SRR3182423.stat
plot-bamstats -p ./SRR3182423 ./SRR3182423.stat

samtools stats --reference ../0.Ref_Genome/human_g1k_v37.fasta ../5.Mapping/SRR3182433.bam > ./SRR3182433.stat
plot-bamstats -p ./SRR3182433 ./SRR3182433.stat

samtools stats --reference ../0.Ref_Genome/human_g1k_v37.fasta ../5.Mapping/SRR3182436.bam > ./SRR3182436.stat
plot-bamstats -p ./SRR3182436 ./SRR3182436.stat

也可以画图看一下质量分布

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn

plot-bamstats -p ./SRR3182419 ./SRR3182419.stat

plot-bamstats -p ./SRR3182423 ./SRR3182423.stat

plot-bamstats -p ./SRR3182433 ./SRR3182433.stat

plot-bamstats -p ./SRR3182436 ./SRR3182436.stat

去除 PCR 扩增引入的重复序列

gatk MarkDuplicatesSpark -I ./4.align/SRR3182419.bam -O ./6.gatk/SRR3182419_dedup.bam --remove-all-duplicates true
gatk MarkDuplicatesSpark -I ./4.align/SRR3182423.bam -O ./6.gatk/SRR3182423_dedup.bam --remove-all-duplicates true
gatk MarkDuplicatesSpark -I ./4.align/SRR3182433.bam -O ./6.gatk/SRR3182433_dedup.bam --remove-all-duplicates true
gatk MarkDuplicatesSpark -I ./4.align/SRR3182436.bam -O ./6.gatk/SRR3182436_dedup.bam --remove-all-duplicates true

准备gatk BaseRecalibrator（基因碱基质量重矫正）需要用到的输入文件：

-O ./6.gatk/table/SRR3182419_recall.table
--known-sites ./0.ref/dbsnp132_20101103.vcf
--reference ./0.ref/human_g1k_v37.fasta

在测序过程中，碱基质量分数（Base Quality Score）可能会受到测序仪错误、化学反应等因素的影响，导致不准确。因此，BaseRecalibrator 通过分析已知变异数据库和测序误差模式，对 BAM/CRAM 文件中的碱基质量分数进行重新校准，从而提高变异检测的准确性。

为人类全基因组构建 samtools fasta 索引和 dictionary

下载变异位点并给变异位点创建索引

使用samtools faidx
gatk IndexFeatureFile

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=25G
#SBATCH --cpus-per-task=5
#SBATCH -t 3-00:00:00
#SBATCH -o establish.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn


samtools faidx ./human_g1k_v37.fasta

picard CreateSequenceDictionary R=./human_g1k_v37.fasta \
O=./human_g1k_v37.dict

wget -c https://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/technical/reference/dbsnp132_20101103.vcf.gz -P ./
gunzip ./dbsnp132_20101103.vcf.gz

gatk IndexFeatureFile --input ./dbsnp132_20101103.vcf

最终生成的文件

校准之前需要先生成BQSR（Base Quality Score Recalibration）校准表

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn

gatk BaseRecalibrator -I ../SRR3182419_dedup.bam \
-O ./SRR3182419_recall.table \
--known-sites ~/workspace/WES/0.Ref_Genome/dbsnp132_20101103.vcf \
--reference ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta

gatk BaseRecalibrator -I ../SRR3182423_dedup.bam \
-O ./SRR3182423_recall.table \
--known-sites ~/workspace/WES/0.Ref_Genome/dbsnp132_20101103.vcf \
--reference ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta

gatk BaseRecalibrator -I ../SRR3182433_dedup.bam \
-O ./SRR3182433_recall.table \
--known-sites ~/workspace/WES/0.Ref_Genome/dbsnp132_20101103.vcf \
--reference ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta

gatk BaseRecalibrator -I ../SRR3182436_dedup.bam \
-O ./SRR3182436_recall.table \
--known-sites ~/workspace/WES/0.Ref_Genome/dbsnp132_20101103.vcf \
--reference ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta

得到的文件

生成校准表之后应用 BQSR

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o apply.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn


gatk ApplyBQSR -I ./SRR3182419_dedup.bam \
-bqsr ./table/SRR3182419_recall.table \
-O ./SRR3182419_bqsr.bam \
-R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta

gatk ApplyBQSR -I ./SRR3182423_dedup.bam \
-bqsr ./table/SRR3182423_recall.table \
-O ./SRR3182423_bqsr.bam \
-R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta

gatk ApplyBQSR -I ./SRR3182433_dedup.bam \
-bqsr ./table/SRR3182433_recall.table \
-O ./SRR3182433_bqsr.bam \
-R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta

gatk ApplyBQSR -I ./SRR3182436_dedup.bam \
-bqsr ./table/SRR3182436_recall.table \
-O ./SRR3182436_bqsr.bam \
-R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta

结果文件：

结果文件

寻找 germline mutations

使用gatk HaplotypeCaller找germline mutations
原教程中对 tumor 和 normal 组织都调用了 HaplotypeCaller，我认为是错误的，因为在肿瘤组织中你是无法区分种系突变（germline mutations）和体细胞突变的（somatic mutations）。

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn

# Case 1: 正常样本SRR3182423
gatk HaplotypeCaller \
   -I ../7.gatk/SRR3182423_bqsr.bam \
   -O ./SRR3182423_raw.vcf \
   -R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta \
   -ERC GVCF \
   --dbsnp ~/workspace/WES/0.Ref_Genome/dbsnp132_20101103.vcf \
   -L ~/workspace/WES/0.Ref_Genome/exon_data.bed

# Case 2: 正常样本SRR3182419
gatk HaplotypeCaller \
   -I ../7.gatk/SRR3182419_bqsr.bam \
   -O ./SRR3182419_raw.vcf \
   -R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta \
   -ERC GVCF \
   --dbsnp ~/workspace/WES/0.Ref_Genome/dbsnp132_20101103.vcf \
   -L ~/workspace/WES/0.Ref_Genome/exon_data.bed

结果文件：

首先使用GenomicDBImport 再使用GenotypeGVCFs进行变异检测

GenomicsDBImport 会把多个 GVCF 文件中的变异数据导入到一个数据库中，这样在处理多个样本的数据时，可以更高效地进行联合分析。数据库会按照染色体（例如 1 到 22，X 和 Y）划分，每个染色体的变异信息存储在对应的数据库文件中。

GenotypeGVCFs 通过从多个样本的 GVCF 文件中获取变异信息，推断每个样本的基因型（例如 0/0、0/1、1/1）。它利用参考基因组、变异数据以及已知的变异信息来推测每个样本在给定位置的等位基因组合，并将结果输出为 VCF 格式。
gatk GatherVcfs 是 GATK 工具包中的一个命令，用于合并多个 VCF 文件。在处理大规模基因组数据时，经常会生成多个染色体（或样本）对应的 VCF 文件，使用这个命令可以将它们合并为一个统一的大文件。

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn


for chr in {1..22} X Y; do
   gatk GenomicsDBImport \
      -R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta \
      -V ./SRR3182423_raw.vcf \
      -V ./SRR3182419_raw.vcf \
      -L ${chr} \
      --genomicsdb-workspace-path ./gvcfs_chr${chr}.db

   gatk GenotypeGVCFs \
      -R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta \
      -V gendb://./gvcfs_chr${chr}.db \
      -O ./gvcfs_chr${chr}.vcf
done

gatk GatherVcfs \
   $(for i in {1..22} X Y; do echo "-I ./gvcfs_chr${i}.vcf"; done) \
   -O ./merged_genotypes.vcf

再进行硬过滤。

VQSR 主要用于大样本，这里只有两个样本，只需要进行硬过滤即可。

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn



gatk VariantFiltration \
  -R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta \
  -V ./merged_genotypes.vcf \
  -O ./hard_filtered.vcf \            # 输出文件名您可自定义
  --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || SOR > 3.0 || ReadPosRankSum < -8.0 || MQRankSum < -12.5" \
  --filter-name "LowQuality" \
  --mask-extension 5 \                # 对过滤位点前后5bp也标记
  --mask-name "NearBadVariant"

gatk SelectVariants \
  -R ~/workspace/WES/0.Ref_Genome/human_g1k_v37.fasta \
  -V ./hard_filtered.vcf \
  -O ./final_pass.vcf \               # 最终结果
  --exclude-filtered true             # 只保留FILTER=PASS的变异

硬过滤完成之后，进行注释并转换成 MAF 文件

# 设置 ANNOVAR 的路径（如果你的 ANNOVAR 在别的路径下，需要修改此变量）
ANNOVAR_DIR=~/workspace/WES/tools/annovar

# 输入 VCF 文件名（可以改成你的输入文件路径）
VCF_FILE=./final_pass.vcf

# 运行 convert2annovar.pl 将 VCF 文件转换为 ANNOVAR 格式
echo "Converting VCF to ANNOVAR input format..."
$ANNOVAR_DIR/convert2annovar.pl \
    -format vcf4 \
    "$VCF_FILE" \
    -allsample \
    --outfile final_pass

# 遍历所有生成的 .avinput 文件，并运行 table_annovar.pl 进行注释
echo "Annotating variants using ANNOVAR..."
for file in final_pass.*.avinput; do
    echo "Processing: $file"
    $ANNOVAR_DIR/table_annovar.pl \
        "$file" \
        $ANNOVAR_DIR/humandb/ \
        -buildver hg19 \
        -out "${file%.avinput}" \
        -remove \
        -protocol refGene \
        -operation g \
        -nastring . \
        -otherinfo
done

echo "Annotation process completed!"

# 遍历每个样本的注释文件进行后处理
for germline in SRR3182419 SRR3182423; do
    echo "Processing germline sample: $germline"

    # 移除Chr开头的行，并添加样本信息
    grep -v '^Chr' ./final_pass.${germline}.hg19_multianno.txt | \
    cut -f 1-24 | \
    awk -v T="$germline" '{print $0"\t"T}' > ./${germline}.annovar.vcf

    echo "Processed ${germline}.annovar.vcf"
done

# 创建表头（只需执行一次）
echo "Creating header..."
head -1 ./final_pass.SRR3182419.hg19_multianno.txt | \
awk '{print $0 "\tTumor_Sample_Barcode"}' > ./header

# 合并所有样本结果
echo "Merging results..."
cat ./header ./*annovar.vcf > ./annovar_merge.vcf

echo "Annotation and merging process completed!"

代码：

rm(list=ls())
setwd("~/workspace/WES/8.Germline_Mutations/gvcf/")
library("maftools")

# germline-------------------------------------------------------------------------


maf_data <- annovarToMaf(
  annovar = "./germline_annovar_merge.vcf",
  MAFobj = T,
  refBuild = "hg19",
  tsbCol = "Tumor_Sample_Barcode",
  table = "refGene"
)


oncoplot(maf = maf_data,
         #top = 30, # 绘制前 30 个变异
         fontSize=1,
         sampleOrder = maf_data@clinical.data$Tumor_Sample_Barcode, # 样本的排序顺序
         showTumorSampleBarcodes = T) # 显示肿瘤样本的名称

plotmafSummary(maf_data)

出图：

oncoplot

plotmafSummary

使用 gatk Mutect2 寻找Somatic mutations

这里仅仅展示一对样品的处理代码

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o Filter.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn


# Step 1: 使用 Mutect2 检测体细胞突变
gatk Mutect2 -R ../0.Ref_Genome/human_g1k_v37.fasta \
-I ../7.gatk/SRR3182436_bqsr.bam -tumor $(basename "../7.gatk/SRR3182436_bqsr.bam" _bqsr.bam) \
-I ../7.gatk/SRR3182419_bqsr.bam -normal $(basename "../7.gatk/SRR3182419_bqsr.bam" _bqsr.bam) \
-L ../0.Ref_Genome/exon_data.bed \
-O ./SRR3182436_mutect2.vcf

# Step 2: 使用 FilterMutectCalls 过滤突变
gatk FilterMutectCalls -R ../0.Ref_Genome/human_g1k_v37.fasta \
-V ./SRR3182436_mutect2.vcf \
-O ./SRR3182436_somatic.vcf

# Step 3: 使用 Perl 过滤 VCF 文件
cat ./SRR3182436_somatic.vcf | perl -alne '{
    if(/^#/) { print }  # 如果是注释行，直接输出
    else {
        next unless $F[6] eq "PASS";  # 过滤掉 FILTER 列不是 "PASS" 的行
        next if $F[0] =~ /_/;         # 过滤掉染色体名称包含下划线的行
        print;                        # 输出符合条件的行
    }
}' > ./SRR3182436_filter.vcf

最终得到两个肿瘤的 VCF 文件：

SRR3182436_filter.vcf
SRR3182433_filter.vcf

使用annovar对 VCF 结果进行注释，看突变发生在哪一个基因上

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn

~/workspace/WES/tools/annovar/table_annovar.pl \
../SRR3182436_filter.vcf \
~/workspace/WES/tools/annovar/humandb/ \
-buildver hg19 \
-out ./SRR3182436 \
-remove \
-protocol refGene \
-operation g \
-nastring . \
-vcfinput

~/workspace/WES/tools/annovar/table_annovar.pl \
../SRR3182433_filter.vcf \
~/workspace/WES/tools/annovar/humandb/ \
-buildver hg19 \
-out ./SRR3182433 \
-remove \
-protocol refGene \
-operation g \
-nastring . \
-vcfinput

image.png

我们可以看到注释结果中出现了哪个基因：CHD5，突变的位置，氨基酸突变的信息等

注释文件的前三行

将 VCF 转变成MAF 文件

#!/bin/bash
#SBATCH -J rawdata
#SBATCH -p dna
#SBATCH -N 1
#SBATCH --mem=40G
#SBATCH --cpus-per-task=8
#SBATCH -t 3-00:00:00
#SBATCH -o raw_data.out
#SBATCH --mail-type=ALL
#SBATCH --mail-user=24211510018@m.fudan.edu.cn


# 以 "Chr" 开头的行移除 添加两列数据 分别是 tumor 和 normal 的编号
grep -v '^Chr' ../annotation/SRR3182433.hg19_multianno.txt | cut -f 1-24 | awk -v T=SRR3182433 -v N=SRR3182423_germline '{print $0"\t"T"\t"N}' > ./SRR3182433.annovar.vcf
grep -v '^Chr' ../annotation/SRR3182436.hg19_multianno.txt | cut -f 1-24 | awk -v T=SRR3182436 -v N=SRR3182419_germline '{print $0"\t"T"\t"N}' > ./SRR3182436.annovar.vcf

# 生成一个 header 文件
# header追加两列 Tumor_Sample_Barcode： 和 Matched_Norm_Sample_Barcode
head -1 ../annotation/SRR3182433.hg19_multianno.txt | awk '{print $0 "\tTumor_Sample_Barcode\tMatched_Norm_Sample_Barcode"}' > ./header

# 合并 SRR3182433.annovar.vcf 和 SRR3182436.annovar.vcf 成 annovar_merge.vcf
cat ./header ./*annovar.vcf > ./annovar_merge.vcf

下游在 RStudio 上使用 maftools 包进行分析

rm(list=ls())
setwd("~/workspace/WES/9.Mutect2/maf/")
library("maftools")

# -------------------------------------------------------------------------


## 读入 annovar_merge.vcf 转成 maf 格式变量 annovar.laml
annovar.laml <- annovarToMaf(annovar = "./annovar_merge.vcf",
                             refBuild='hg19',
                             tsbCol = 'Tumor_Sample_Barcode',
                             table = 'refGene',
                             MAFobj = T)

save(annovar.laml,file = 'annovar_laml.Rdata')


# -------------------------------------------------------------------------



## 绘制突变频率汇总图
postscript('plotmafSummary_annovar.eps', width = 5, height = 5, horizontal = FALSE)
plotmafSummary(maf = annovar.laml,
               rmOutlier = TRUE, # 移除异常值
               showBarcodes = T, # 显示条形码来表示每个样本的突变频率
               textSize = 0.4, # 设置文本的大小
               addStat = 'median', # 添加统计信息，这里选择显示中位数
               dashboard = TRUE, # 显示仪表盘式图表
               titvRaw = FALSE) # 不显示硬件替代碱基转换比率（Transition/Transversion）

# Variants per sample (Median: 404) 样品突变数量的中位数。
# TMB tumor mutations burden
# variant classification summary 是每一种类型的 mutation，除以样本数量
# top 10 mutated genes 这个图表示，每个基因上有几个突变，如图所示，TTN 有 4 个突变
dev.off()


# -------------------------------------------------------------------------

oncoplot(maf = annovar.laml, top = 10)


# -------------------------------------------------------------------------

## 生成癌症样本的分析图：可视化肿瘤样本中的基因突变信息
postscript('oncoplot_top30_annovar.eps', width = 5, height = 5, horizontal = FALSE)
oncoplot(maf = annovar.laml,
         top = 30, # 绘制前 30 个变异
         fontSize=0.5,
         sampleOrder = annovar.laml@clinical.data$Tumor_Sample_Barcode, # 样本的排序顺序
         showTumorSampleBarcodes = T) # 显示肿瘤样本的名称
dev.off()


# -------------------------------------------------------------------------

## 可视化 KDM4A 基因突变
postscript('KDM4A_annovar.eps', width = 10, height = 3, horizontal = FALSE)
maftools::lollipopPlot(maf = annovar.laml,
                       gene = 'KDM4A', # 指定要绘制的基因，这里是 "SP140"
                       AACol = 'AAChange.refGene', # 显示的氨基酸改变信息的列，这里是 "AAChange.refGene"。
                       labelPos = 'all', # 指定标签显示的位置，这里设置为 "all"，表示显示所有标签。
                       labPosSize = 0.8,
                       repel = TRUE,
                       showDomainLabel = FALSE)

plotsummary

oncoplot

lollipopPlot

将所有的生物学重复和技术重复都考虑上

# 体细胞突变检测与过滤流程
for tumor_id in SRR3182434 SRR3182435 SRR3182444; do
    # 1. Mutect2检测
    gatk Mutect2 \
        -R ../0.Ref_Genome/human_g1k_v37.fasta \
        -I ../7.gatk/${tumor_id}_bqsr.bam -tumor ${tumor_id} \
        -I ../7.gatk/SRR3182423_bqsr.bam -normal SRR3182423 \
        -L ../0.Ref_Genome/exon_data.bed \
        -O ./${tumor_id}_mutect2.vcf

    # 2. GATK过滤
    gatk FilterMutectCalls \
        -R ../0.Ref_Genome/human_g1k_v37.fasta \
        -V ./${tumor_id}_mutect2.vcf \
        -O ./${tumor_id}_somatic.vcf

    # 3. Perl过滤（严格按您的要求）
    cat ./${tumor_id}_somatic.vcf | perl -alne '{
        if(/^#/){print}
        else{
            next unless $F[6] eq "PASS";
            next if $F[0] =~ /_/;
            print
        }
    }' > ./${tumor_id}_filter.vcf  # 已改为_filter.vcf

    # 简洁日志
    echo "完成 ${tumor_id}：输入变异 $(grep -vc '^#' ${tumor_id}_somatic.vcf) → 输出 $(grep -vc '^#' ${tumor_id}_filter.vcf)"
done

注释

for sample in SRR3182434 SRR3182435 SRR3182444; do
    ~/workspace/WES/tools/annovar/table_annovar.pl \
        ../${sample}_filter.vcf \
        ~/workspace/WES/tools/annovar/humandb/ \
        -buildver hg19 \
        -out ./${sample} \
        -remove \
        -protocol refGene \
        -operation g \
        -nastring . \
        -vcfinput
done

处理成maftools 能够读取的形式

for tumor in SRR3182433 SRR3182434 SRR3182435 SRR3182444; do
    # 移除Chr开头的行并添加样本信息
    grep -v '^Chr' ../annotation/${tumor}.hg19_multianno.txt | \
    cut -f 1-24 | \
    awk -v T="$tumor" -v N="SRR3182423_germline" '{print $0"\t"T"\t"N}' > ./${tumor}.annovar.vcf
done

# 处理Case 2的样本（匹配正常样本SRR3182419）
for tumor in SRR3182436; do
    # 移除Chr开头的行并添加样本信息
    grep -v '^Chr' ../annotation/${tumor}.hg19_multianno.txt | \
    cut -f 1-24 | \
    awk -v T="$tumor" -v N="SRR3182419_germline" '{print $0"\t"T"\t"N}' > ./${tumor}.annovar.vcf
done

# 创建表头（只需执行一次）
head -1 ../annotation/SRR3182433.hg19_multianno.txt | \
awk '{print $0 "\tTumor_Sample_Barcode\tMatched_Norm_Sample_Barcode"}' > ./header

# 合并所有样本结果
cat ./header ./*annovar.vcf > ./annovar_merge.vcf

出图：

plotSummary

oncoplot

在 oncoplot 我们可以看到一共有 5 个样本，除了SRR3182436是 case2 的样品，其他的都是 case1 的样品，包括 3 个生物学重复和 1 个技术重复。

lollipopPlot

参考：

https://bohrium.dp.tech/notebooks/2112719645

WES 实战

WES 实战

全外显子组测序（Whole Exome Sequencing, WES）技术简介

定义与原理

Note

技术特点

Data description

Download exon data and reference genome

Download Rawdata

Trim

Mapping

对 Mapping 进行质控

去除 PCR 扩增引入的重复序列

准备gatk BaseRecalibrator（基因碱基质量重矫正）需要用到的输入文件：

为人类全基因组构建 samtools fasta 索引和 dictionary

下载变异位点并给变异位点创建索引

校准之前需要先生成BQSR（Base Quality Score Recalibration）校准表

生成校准表之后应用 BQSR

寻找 germline mutations

首先使用GenomicDBImport 再使用GenotypeGVCFs进行变异检测

使用 gatk Mutect2 寻找Somatic mutations

使用annovar对 VCF 结果进行注释，看突变发生在哪一个基因上

将 VCF 转变成MAF 文件

下游在 RStudio 上使用 maftools 包进行分析

参考：

推荐阅读更多精彩内容

WES 实战

全外显子组测序（Whole Exome Sequencing, WES）技术简介

定义与原理

Note

技术特点

Data description

Download exon data and reference genome

Download Rawdata

Trim

Mapping

对 Mapping 进行质控

去除 PCR 扩增引入的重复序列

准备gatk BaseRecalibrator（基因碱基质量重矫正）需要用到的输入文件：

为人类全基因组构建 samtools fasta 索引 和 dictionary

下载变异位点并给变异位点创建索引

校准之前需要先生成BQSR（Base Quality Score Recalibration） 校准表

生成校准表之后应用 BQSR

寻找 germline mutations

首先使用GenomicDBImport 再使用GenotypeGVCFs进行变异检测

使用 gatk Mutect2 寻找Somatic mutations

使用annovar对 VCF 结果进行注释，看突变发生在哪一个基因上

将 VCF 转变成MAF 文件

下游在 RStudio 上使用 maftools 包进行分析

参考：

推荐阅读更多精彩内容

为人类全基因组构建 samtools fasta 索引和 dictionary

校准之前需要先生成BQSR（Base Quality Score Recalibration）校准表