新建
探索全外显子测序数据分析流程：从数据预处理到变异分析
孙楠
推荐镜像 :wes-java17-gatk:2023.09.11
推荐机型 :c32_m128_cpu
赞 3
5
4
目录
数据集
wes_0.rawdata(v2)
wes_0.ref(v1)
探索全外显子测序数据分析流程：从数据预处理到变异分析
代码
文本
全外显子测序（Whole Exome Sequencing，WES）是一种高通量测序技术，用于对生物体的所有外显子区域进行测序。通过全外显子测序，可以有效地检测基因组中与疾病相关的变异，如突变、缺失、插入等。全外显子测序数据分析涉及多个层面的内容，如变异检测和注释： 分析外显子测序数据的首要任务是检测样本中存在的基因组变异，包括单核苷酸变异（SNV）、插入缺失（Indel）等。这些变异可能与疾病相关，需要进行详细的注释，包括变异类型、位置、影响的基因、功能等信息。
更多细节推荐观看 B站视频（WES讲解）
全外显子测序数据预处理（Step 0 ～ Step 6）是确保从原始测序数据到高质量测序数据转化的关键步骤。本文将演示如何下载人类全基因组作为参考基因组，以及从 SRA 数据库检索测序数据，来成功运行 WES 数据分析流程。这一流程将涵盖从原始数据到最终分析结果的处理步骤，让您深入了解WES数据的处理过程。
首先，您将下载 GRCh37 版本的人类全基因组，以便将测序数据与参考基因组进行比对和分析。接着，您将从SRA数据库中检索登录号为 SRP070662 的测序数据，这些数据将作为我们的分析对象。在数据准备阶段，您将使用数据质控工具（FastQC）来评估测序数据的质量、碱基分布和序列长度等，以确保数据质量符合分析要求。接下来，您将使用去接头工具（TrimGalore）来处理原始数据，去除接头序列和低质量碱基，以减少后续分析中的误差。随后，您将利用比对工具（BWA）将经过质控和去接头的测序数据与 GRCh37 参考基因组进行比对，从而获取测序数据的位置信息。在得到比对结果后，您将使用 GATK 执行去重步骤、对测序数据进行标记碱基质量值（BQSR）的校正，提升数据质量，使得后续的变异检测等分析更加可靠和准确。
初步了解推荐观看 B站视频（二代基因测序之全外显子测序数据预处理流程）
在全外显子测序数据预处理流程的基础上，我们将进入变异检测阶段（Step 7 ～ Step 8）。在这一阶段，我们将利用 GATK 工具对经过上面质控、去接头、比对和质量校正的测序数据进行分析，以检测样本中的遗传变异。通过使用 gatk HaplotypeCaller 命令，我们将识别单核苷酸变异（SNVs）和小插入缺失（Indels），并生成原始的变异信息。随后，我们使用 gatk Mutect2 命令来识别样本中的突变，同时区分肿瘤样本和正常样本之间的变异。之后，通过 gatk FilterMutectCalls 命令进行过滤，排除掉低质量的突变，以及可能的假阳性结果。最终，我们筛选出具有生物学意义的高质量突变，得到一个经过过滤和注释的变异列表，为后续的生物学解释和临床研究提供有价值的信息。整个变异检测过程将帮助我们更深入地理解样本中的遗传变异，为疾病研究和治疗提供重要参考。
初步了解推荐观看 B站视频（二代基因测序之全外显子测序体细胞突变分析GATK流程）
现在，让我们开始这个流程之旅吧！无论您是新手还是有经验的分析师，这个流程都将为您提供实用的指导，帮助您成功地运行 WES 数据分析流程。让我们一起探索如何从原始数据中获取有价值的信息，为进一步的生物信息学研究和解释奠定坚实的基础。
代码
文本
目标
掌握 WES 常用分析流程，包括质控、去接头、序列比对、去重、BQSR校正、变异检测variant calling等。
目录
代码
文本
📖 上手指南
若无特殊说明，本文档所有命令需要在 Linux Shell 中执行。若读者想独立运行下面 Linux 命令，点击文档最上面的“Open in Bohrium”，连接镜像“wes-java17-gatk:2023.09.11”，然后点击“SSH连接到节点”，并按步骤粘贴命令到 Linux Shell 执行。
镜像中包含如下软件包，可以通过粘贴如下命令到 Linux Shell 检查各个软件包的安装：
fastqc --version
multiqc --version
trim_galore --version
bwa
samtools --version
gatk --help
gatk HaplotypeCaller --help
picard -h
代码
文本
mkdir 0.rawdata 0.ref 1.QC 2.clean_QC 3.bwa_index 4.align 5.stat 6.gatk 7.Germline_Mutations 8.Mutect2 tools
mkdir 1.QC/fast_QC 1.QC/multi_QC
mkdir 2.clean_QC/fast_QC 2.clean_QC/multi_QC 2.clean_QC/trim_galore_QC
mkdir 6.gatk/table
mkdir 7.Germline_Mutations/gvcf/
mkdir 8.Mutect2/annotation 8.Mutect2/maf
代码
文本
将随附数据“wes_0.rawdata”移动到“0.rawdata”文件夹中，随附数据“wes_0.ref”移动到“0.ref”文件夹中。
代码
文本
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
代码
文本
2. 下载人类外显子数据
在 UCSC table browser 中，选择要下载的数据。这里的选择包括：
- Assembly = Feb. 2009 (GRCh37/hg19)：选择基因组的特定版本，即2009年2月的GRCh37/hg19版本。
- Track = Ensembl Genes：选择要检索的基因组注释数据的特定跟踪（Ensembl基因）。
- Output Format = BED：选择要将数据输出为BED格式，BED是一种常用的基因组坐标文件格式。
- Output = Coding Exons：选择只输出编码外显子数据。
移除第一列中的"chr"：
下载的BED文件中的第一列通常包含染色体名称，以"chr"开头（例如：chr1、chr2等）。在某些情况下，需要将这些"chr"前缀移除。为此，使用了一个命令行工具 awk，它允许对文本文件进行处理。
代码
文本
# 筛选出 exon_download.bed 文件中第一列属于 exon_download_24chr.txt 列表的行
cd 0.ref
grep -Fwf exon_download_24chr.txt exon_download.bed > filtered_exon_download.bed
代码
文本
其中 exon_download_24chr.txt 文件只包含 24 条染色体信息（我们只分析这 24 个染色体分子）： chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr2 chr20 chr21 chr22 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chrX chrY
代码
文本
awk '{gsub("chr", "", $0)}1' filtered_exon_download.bed > exon_data.bed
cd ..
代码
文本
exon_data.bed 每一行表示一个外显子的信息：
数据解释：
- 第一列：染色体编号，例如："1" 表示第1号染色体。
- 第二列：外显子的起始位置，基于染色体上的坐标（起始坐标）。
- 第三列：外显子的结束位置，基于染色体上的坐标（终止坐标）。
- 第四列：外显子的标识符，可能包括基因ID、外显子的编号等。
- 第五列：分数， 此处为 "0"，表示没有特定的分数信息。
- 第六列：外显子的方向，"+" 表示正链，"-" 表示负链。
代码
文本
3. 下载 WES 数据
从 SRA 数据库检索到登录号是 SRP070662 的测序数据，选择 'Send results to Run selector' 按钮，进入该 id 的详细数据列表。再选择 Assay Type 是 WXS 的测序数据，一共显示有 30 个样本：6 个case，每个 case 有 5 个样本，分别是 Germline，3 个 Biological replicate：A、B、C，1 个 Technical replicate。
我们选择下面 4 套数据集（都双端测序 PAIRED，Average Spot Length: 151）：
- case 1：Germline（id是SRR3182423）、Biological replicate A （id是SRR3182433）
- case 2：Germline（id是SRR3182419）、Biological replicate A （id是SRR3182436）
有很多工具或者方法可以下载，如 sratoolkit或 sra-explorer 提到的工具。我们已经将数据下载到目录中，读者无需下载（从下面 4 个链接可以看到每个样本对应 3 个文件，只需要下载*_1.fastq和*_2.fastq）。
代码
文本
ls -lh 0.rawdata
代码
文本
代码
文本
fastqc ./0.rawdata/*.fastq --outdir ./1.QC/fast_QC
代码
文本
运行上面命令在 ./1.QC/fast_QC/ 文件夹为每个 .fastq 生成 1 个 html 文件和 1 个 zip 压缩文件。
代码
文本
ls -lh 1.QC/fast_QC
代码
文本
其中 html 文件是质控报告，zip 文件是 html 文件图表对应的类。样本数较多，可以用 multiqc 对所有样本的结果整合：
代码
文本
multiqc ./1.QC/fast_QC/*.zip -o ./1.QC/multi_QC
代码
文本
ls -lh 1.QC/multi_QC
代码
文本
质控报告部分截图如下：
从完整的质控报告看到一些 reads 的质量值较差或测到了接头，因此需要去除掉这些质量较差的 reads 或接头，使用软件 trim- galore。
代码
文本
Step 2：去接头
利用工具 Trim Galore
代码
文本
trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 24 \
-o ./2.clean_QC/trim_galore_QC ./0.rawdata/SRR3182419_1.fastq ./0.rawdata/SRR3182419_2.fastq \
>> ./2.clean_QC/trim_galore_QC/SRR3182419_trim.log 2>&1
代码
文本
trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 24 \
-o ./2.clean_QC/trim_galore_QC ./0.rawdata/SRR3182423_1.fastq ./0.rawdata/SRR3182423_2.fastq \
>> ./2.clean_QC/trim_galore_QC/SRR3182423_trim.log 2>&1
代码
文本
trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 24 \
-o ./2.clean_QC/trim_galore_QC ./0.rawdata/SRR3182433_1.fastq ./0.rawdata/SRR3182433_2.fastq \
>> ./2.clean_QC/trim_galore_QC/SRR3182433_trim.log 2>&1
代码
文本
trim_galore --paired -q 28 --phred33 --length 30 --stringency 3 --gzip --cores 24 \
-o ./2.clean_QC/trim_galore_QC ./0.rawdata/SRR3182436_1.fastq ./0.rawdata/SRR3182436_2.fastq \
>> ./2.clean_QC/trim_galore_QC/SRR3182436_trim.log 2>&1
代码
文本
ls -lh 2.clean_QC/trim_galore_QC
代码
文本
- SRR3182419_1_val_1.fq.gz 表示经过质量修剪后的有效数据，即去除了低质量的部分。
代码
文本
再质控：
代码
文本
fastqc ./2.clean_QC/trim_galore_QC/*.fq.gz --outdir ./2.clean_QC/fast_QC
代码
文本
multiqc ./2.clean_QC/fast_QC/*.zip -o ./2.clean_QC/multi_QC
代码
文本
ls -lh 2.clean_QC/multi_QC
代码
文本
质控报告部分截图如下：
如何查看碱基质量，推荐阅读微信文章（从零开始完整学习全基因组测序数据分析：第2节 FASTA和FASTQ） 和 微信文章（从零开始完整学习全基因组测序数据分析：第3节 数据质控）。
代码
文本
Step 3: 构建参考基因组的 BWA 索引
利用工具命令 bwa index
代码
文本
bwa index -p ./3.bwa_index/human_g1k_v37 ./0.ref/human_g1k_v37.fasta
代码
文本
ls -lh 3.bwa_index
代码
文本
生成的 BWA 索引文件说明：
- human_g1k_v37.amb: 它存储了参考基因组序列的一些元信息，如序列长度、名称等。
- human_g1k_v37.ann: 它存储了一些额外的注释信息，有助于在比对过程中进行一些元数据的操作。
- human_g1k_v37.bwt: 这是Burrows-Wheeler变换的数据，是BWA索引的核心部分之一。BWT是一种数据转换算法，用于将原始DNA序列进行变换，以便在比对时能够更高效地进行匹配。
- human_g1k_v37.pac: 这是FM索引的数据，也是BWA索引的核心部分之一。FM索引是一种数据结构，用于实现快速查找和比对DNA序列。
- human_g1k_v37.sa: 这是后缀数组的数据，也是BWA索引的核心部分之一。后缀数组是一种用于存储DNA序列的子串信息的数据结构，用于在比对时快速查找匹配。
代码
文本
bwa mem -M -t 24 -R "@RG\tID:SRR3182419\tSM:SRR3182419\tLB:WXS\tPL:Illummina" ./3.bwa_index/human_g1k_v37 \
./2.clean_QC/trim_galore_QC/SRR3182419_1_val_1.fq.gz \
./2.clean_QC/trim_galore_QC/SRR3182419_2_val_2.fq.gz \
> ./4.align/SRR3182419.sam
代码
文本
bwa mem -M -t 24 -R "@RG\tID:SRR3182423\tSM:SRR3182423\tLB:WXS\tPL:Illummina" ./3.bwa_index/human_g1k_v37 \
./2.clean_QC/trim_galore_QC/SRR3182423_1_val_1.fq.gz \
./2.clean_QC/trim_galore_QC/SRR3182423_2_val_2.fq.gz \
-o ./4.align/SRR3182423.sam
代码
文本
bwa mem -M -t 24 -R "@RG\tID:SRR3182433\tSM:SRR3182433\tLB:WXS\tPL:Illummina" ./3.bwa_index/human_g1k_v37 \
./2.clean_QC/trim_galore_QC/SRR3182433_1_val_1.fq.gz \
./2.clean_QC/trim_galore_QC/SRR3182433_2_val_2.fq.gz \
-o ./4.align/SRR3182433.sam
代码
文本
bwa mem -M -t 24 -R "@RG\tID:SRR3182436\tSM:SRR3182436\tLB:WXS\tPL:Illummina" ./3.bwa_index/human_g1k_v37 \
./2.clean_QC/trim_galore_QC/SRR3182436_1_val_1.fq.gz \
./2.clean_QC/trim_galore_QC/SRR3182436_2_val_2.fq.gz \
-o ./4.align/SRR3182436.sam
代码
文本
运行上面 4 个命令在 ./4.align 文件夹生成了 4 个文件：
- SRR3182419.sam
- SRR3182423.sam
- SRR3182433.sam
- SRR3182436.sam
sam 文件是纯文本格式，占用存储空间，可通过 samtools 转换成二进制的 bam 文件，记录了每一条 reads 比对到参考基因组的结果。两者记录的信息是一样的。
代码
文本
samtools sort -@ 10 ./4.align/SRR3182419.sam -o ./4.align/SRR3182419.bam
samtools sort -@ 10 ./4.align/SRR3182423.sam -o ./4.align/SRR3182423.bam
samtools sort -@ 10 ./4.align/SRR3182433.sam -o ./4.align/SRR3182433.bam
samtools sort -@ 10 ./4.align/SRR3182436.sam -o ./4.align/SRR3182436.bam
代码
文本
ls -lh 4.align
代码
文本
代码
文本
Step 5: 比对结果进行质控
利用工具 samtools stats、plot-bamstats
代码
文本
samtools stats --reference ./0.ref/human_g1k_v37.fasta ./4.align/SRR3182419.bam > ./5.stat/SRR3182419.stat
#plot-bamstats -p ./5.stat/SRR3182419 ./5.stat/SRR3182419.stat
代码
文本
samtools stats --reference ./0.ref/human_g1k_v37.fasta ./4.align/SRR3182423.bam > ./5.stat/SRR3182423.stat
#plot-bamstats -p ./5.stat/SRR3182423 ./5.stat/SRR3182423.stat
代码
文本
samtools stats --reference ./0.ref/human_g1k_v37.fasta ./4.align/SRR3182433.bam > ./5.stat/SRR3182433.stat
#plot-bamstats -p ./5.stat/SRR3182433 ./5.stat/SRR3182433.stat
代码
文本
samtools stats --reference ./0.ref/human_g1k_v37.fasta ./4.align/SRR3182436.bam > ./5.stat/SRR3182436.stat
#plot-bamstats -p ./5.stat/SRR3182436 ./5.stat/SRR3182436.stat
代码
文本
1. 删除重复序列
利用工具命令 gatk MarkDuplicates
代码
文本
gatk MarkDuplicatesSpark -help
代码
文本
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
代码
文本
ls -lh 6.gatk
代码
文本
代码
文本
查看删除重复序列前后 reads 数目变化：
代码
文本
samtools view ./4.align/SRR3182419.bam |wc -l
samtools view ./6.gatk/SRR3182419_dedup.bam |wc -l
代码
文本
代码
文本
代码
文本
samtools faidx ./0.ref/human_g1k_v37.fasta
代码
文本
运行上面命令在 ./0.ref 文件夹生成了 1 个文件：
- human_g1k_v37.fasta.fai
代码
文本
picard CreateSequenceDictionary R=./0.ref/human_g1k_v37.fasta \
O=./0.ref/human_g1k_v37.dict
代码
文本
运行上面命令在 ./0.ref 文件夹生成了 1 个文件：
- human_g1k_v37.dict
代码
文本
wget -c https://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/technical/reference/dbsnp132_20101103.vcf.gz -P ./0.ref
gunzip ./0.ref/dbsnp132_20101103.vcf.gz
代码
文本
dbsnp132_20101103.vcf.gz 包含了人类基因组中已知的单核苷酸多态性（Single Nucleotide Polymorphism，SNP）的信息。VCF 文件是一种常用的基因组变异数据格式。
给变异位点创建索引：
代码
文本
gatk IndexFeatureFile --input ./0.ref/dbsnp132_20101103.vcf
代码
文本
运行上面命令在 ./0.ref 文件夹生成了 1 个文件：
- dbsnp132_20101103.vcf.idx
代码
文本
ls -lh 0.ref
代码
文本
代码
文本
wget -c https://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/technical/reference/dbsnp132_20101103.vcf.gz -P ./0.ref
代码
文本
gatk BaseRecalibrator -I ./6.gatk/SRR3182419_dedup.bam \
-O ./6.gatk/table/SRR3182419_recall.table \
--known-sites ./0.ref/dbsnp132_20101103.vcf \
--reference ./0.ref/human_g1k_v37.fasta
代码
文本
gatk BaseRecalibrator -I ./6.gatk/SRR3182423_dedup.bam \
-O ./6.gatk/table/SRR3182423_recall.table \
--known-sites ./0.ref/dbsnp132_20101103.vcf \
--reference ./0.ref/human_g1k_v37.fasta
代码
文本
gatk BaseRecalibrator -I ./6.gatk/SRR3182433_dedup.bam \
-O ./6.gatk/table/SRR3182433_recall.table \
--known-sites ./0.ref/dbsnp132_20101103.vcf \
--reference ./0.ref/human_g1k_v37.fasta
代码
文本
gatk BaseRecalibrator -I ./6.gatk/SRR3182436_dedup.bam \
-O ./6.gatk/table/SRR3182436_recall.table \
--known-sites ./0.ref/dbsnp132_20101103.vcf \
--reference ./0.ref/human_g1k_v37.fasta
代码
文本
ls -lh 6.gatk/table
代码
文本
代码
文本
gatk ApplyBQSR -I ./6.gatk/SRR3182419_dedup.bam \
-bqsr ./6.gatk/table/SRR3182419_recall.table \
-O ./6.gatk/SRR3182419_bqsr.bam \
-R ./0.ref/human_g1k_v37.fasta
代码
文本
gatk ApplyBQSR -I ./6.gatk/SRR3182423_dedup.bam \
-bqsr ./6.gatk/table/SRR3182423_recall.table \
-O ./6.gatk/SRR3182423_bqsr.bam \
-R ./0.ref/human_g1k_v37.fasta
代码
文本
gatk ApplyBQSR -I ./6.gatk/SRR3182433_dedup.bam \
-bqsr ./6.gatk/table/SRR3182433_recall.table \
-O ./6.gatk/SRR3182433_bqsr.bam \
-R ./0.ref/human_g1k_v37.fasta
代码
文本
gatk ApplyBQSR -I ./6.gatk/SRR3182436_dedup.bam \
-bqsr ./6.gatk/table/SRR3182436_recall.table \
-O ./6.gatk/SRR3182436_bqsr.bam \
-R ./0.ref/human_g1k_v37.fasta
代码
文本
ls -lh ./6.gatk/*_bqsr.*
代码
文本
代码
文本
mkdir 7.Germline_Mutations
代码
文本
gatk HaplotypeCaller -I ./6.gatk/SRR3182419_bqsr.bam \
-O ./7.Germline_Mutations/SRR3182419_raw.vcf \
-R ./0.ref/human_g1k_v37.fasta \
-ERC GVCF \
--dbsnp ./0.ref/dbsnp132_20101103.vcf \
-L ./0.ref/exon_data.bed
代码
文本
gatk HaplotypeCaller -I ./6.gatk/SRR3182423_bqsr.bam \
-O ./7.Germline_Mutations/SRR3182423_raw.vcf \
-R ./0.ref/human_g1k_v37.fasta \
-ERC GVCF \
--dbsnp ./0.ref/dbsnp132_20101103.vcf \
-L ./0.ref/exon_data.bed
代码
文本
gatk HaplotypeCaller -I ./6.gatk/SRR3182433_bqsr.bam \
-O ./7.Germline_Mutations/SRR3182433_raw.vcf \
-R ./0.ref/human_g1k_v37.fasta \
-ERC GVCF \
--dbsnp ./0.ref/dbsnp132_20101103.vcf \
-L ./0.ref/exon_data.bed
代码
文本
gatk HaplotypeCaller -I ./6.gatk/SRR3182436_bqsr.bam \
-O ./7.Germline_Mutations/SRR3182436_raw.vcf \
-R ./0.ref/human_g1k_v37.fasta \
-ERC GVCF \
--dbsnp ./0.ref/dbsnp132_20101103.vcf \
-L ./0.ref/exon_data.bed
代码
文本
ls -lh 7.Germline_Mutations/*.vcf*
代码
文本
上面执行了四个 HaplotypeCaller 命令，每个命令针对一个样本（SRR3182419、SRR3182423、SRR3182433、SRR3182436）检测原始变异，生成一个 Genomic VCF（gVCF） 文件。
代码
文本
mkdir ./7.Germline_Mutations/gvcf
代码
文本
for chr in {1..22} X Y
do
echo <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8889em;vertical-align:-0.1944em;"></span><span class="mord mathnormal">c</span><span class="mord mathnormal">h</span><span class="mord mathnormal" style="margin-right:0.02778em;">r</span><span class="mord mathnormal" style="margin-right:0.03588em;">g</span><span class="mord mathnormal">a</span><span class="mord mathnormal">t</span><span class="mord mathnormal" style="margin-right:0.03148em;">k</span><span class="mord mathnormal">G</span><span class="mord mathnormal">e</span><span class="mord mathnormal">n</span><span class="mord mathnormal">o</span><span class="mord mathnormal">mi</span><span class="mord mathnormal" style="margin-right:0.02778em;">csD</span><span class="mord mathnormal" style="margin-right:0.05017em;">B</span><span class="mord mathnormal" style="margin-right:0.07847em;">I</span><span class="mord mathnormal">m</span><span class="mord mathnormal">p</span><span class="mord mathnormal" style="margin-right:0.02778em;">or</span><span class="mord mathnormal">t</span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mbin">−</span><span class="mspace" style="margin-right:0.2222em;"></span></span><span class="base"><span class="strut" style="height:1.0361em;vertical-align:-0.2861em;"></span><span class="mord mathnormal" style="margin-right:0.00773em;">R</span><span class="mord">./0.</span><span class="mord mathnormal">re</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span><span class="mord">/</span><span class="mord mathnormal">h</span><span class="mord mathnormal">u</span><span class="mord mathnormal">ma</span><span class="mord"><span class="mord mathnormal">n</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight" style="margin-right:0.03588em;">g</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.2861em;"><span></span></span></span></span></span></span><span class="mord">1</span><span class="mord"><span class="mord mathnormal" style="margin-right:0.03148em;">k</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-left:-0.0315em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight" style="margin-right:0.03588em;">v</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord">37.</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span><span class="mord mathnormal">a</span><span class="mord mathnormal">s</span><span class="mord mathnormal">t</span><span class="mord mathnormal">a</span></span></span></span>(ls ./7.Germline_Mutations/*raw.vcf | awk '{print "-V "$0" "}') -L ${chr} --genomicsdb-workspace-path ./7.Germline_Mutations/gvcf/gvcfs_chr<span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8889em;vertical-align:-0.1944em;"></span><span class="mord"><span class="mord mathnormal">c</span><span class="mord mathnormal">h</span><span class="mord mathnormal" style="margin-right:0.02778em;">r</span></span><span class="mord">.</span><span class="mord mathnormal">d</span><span class="mord mathnormal">b</span><span class="mord mathnormal" style="margin-right:0.03588em;">g</span><span class="mord mathnormal">a</span><span class="mord mathnormal">t</span><span class="mord mathnormal" style="margin-right:0.03148em;">k</span><span class="mord mathnormal">G</span><span class="mord mathnormal">e</span><span class="mord mathnormal">n</span><span class="mord mathnormal">o</span><span class="mord mathnormal">t</span><span class="mord mathnormal" style="margin-right:0.03588em;">y</span><span class="mord mathnormal">p</span><span class="mord mathnormal">e</span><span class="mord mathnormal">G</span><span class="mord mathnormal" style="margin-right:0.22222em;">V</span><span class="mord mathnormal" style="margin-right:0.13889em;">CF</span><span class="mord mathnormal">s</span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mbin">−</span><span class="mspace" style="margin-right:0.2222em;"></span></span><span class="base"><span class="strut" style="height:1.0361em;vertical-align:-0.2861em;"></span><span class="mord mathnormal" style="margin-right:0.00773em;">R</span><span class="mord">./0.</span><span class="mord mathnormal">re</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span><span class="mord">/</span><span class="mord mathnormal">h</span><span class="mord mathnormal">u</span><span class="mord mathnormal">ma</span><span class="mord"><span class="mord mathnormal">n</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight" style="margin-right:0.03588em;">g</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.2861em;"><span></span></span></span></span></span></span><span class="mord">1</span><span class="mord"><span class="mord mathnormal" style="margin-right:0.03148em;">k</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-left:-0.0315em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight" style="margin-right:0.03588em;">v</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord">37.</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span><span class="mord mathnormal">a</span><span class="mord mathnormal">s</span><span class="mord mathnormal">t</span><span class="mord mathnormal">a</span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mbin">−</span><span class="mspace" style="margin-right:0.2222em;"></span></span><span class="base"><span class="strut" style="height:0.8889em;vertical-align:-0.1944em;"></span><span class="mord mathnormal" style="margin-right:0.22222em;">V</span><span class="mord mathnormal" style="margin-right:0.03588em;">g</span><span class="mord mathnormal">e</span><span class="mord mathnormal">n</span><span class="mord mathnormal">d</span><span class="mord mathnormal">b</span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mrel">:</span><span class="mspace" style="margin-right:0.2778em;"></span></span><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mord">//7.</span><span class="mord mathnormal">G</span><span class="mord mathnormal" style="margin-right:0.02778em;">er</span><span class="mord mathnormal">m</span><span class="mord mathnormal" style="margin-right:0.01968em;">l</span><span class="mord mathnormal">in</span><span class="mord"><span class="mord mathnormal">e</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3283em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight" style="margin-right:0.10903em;">M</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal">u</span><span class="mord mathnormal">t</span><span class="mord mathnormal">a</span><span class="mord mathnormal">t</span><span class="mord mathnormal">i</span><span class="mord mathnormal">o</span><span class="mord mathnormal">n</span><span class="mord mathnormal">s</span><span class="mord">/</span><span class="mord mathnormal" style="margin-right:0.03588em;">gv</span><span class="mord mathnormal">c</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span><span class="mord">/</span><span class="mord mathnormal" style="margin-right:0.03588em;">gv</span><span class="mord mathnormal">c</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span><span class="mord"><span class="mord mathnormal">s</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight">c</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal">h</span><span class="mord mathnormal" style="margin-right:0.02778em;">r</span></span></span></span>{chr}.db -O ./7.Germline_Mutations/gvcf/gvcfs_chr${chr}.vcf
done
代码
文本
ls -lh ./7.Germline_Mutations/gvcf
代码
文本
代码
文本
gatk GatherVcfs \
<span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mopen">(</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span><span class="mord mathnormal" style="margin-right:0.02778em;">or</span><span class="mord mathnormal">iin</span><span class="mord"><span class="mord">1..22</span></span><span class="mord mathnormal" style="margin-right:0.07847em;">X</span><span class="mord mathnormal" style="margin-right:0.22222em;">Y</span><span class="mpunct">;</span><span class="mspace" style="margin-right:0.1667em;"></span><span class="mord mathnormal">d</span><span class="mord mathnormal">oec</span><span class="mord mathnormal">h</span><span class="mord mathnormal">o</span><span class="mord">"</span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mbin">−</span><span class="mspace" style="margin-right:0.2222em;"></span></span><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mord mathnormal" style="margin-right:0.07847em;">I</span><span class="mord">./7.</span><span class="mord mathnormal">G</span><span class="mord mathnormal" style="margin-right:0.02778em;">er</span><span class="mord mathnormal">m</span><span class="mord mathnormal" style="margin-right:0.01968em;">l</span><span class="mord mathnormal">in</span><span class="mord"><span class="mord mathnormal">e</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3283em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight" style="margin-right:0.10903em;">M</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal">u</span><span class="mord mathnormal">t</span><span class="mord mathnormal">a</span><span class="mord mathnormal">t</span><span class="mord mathnormal">i</span><span class="mord mathnormal">o</span><span class="mord mathnormal">n</span><span class="mord mathnormal">s</span><span class="mord">/</span><span class="mord mathnormal" style="margin-right:0.03588em;">gv</span><span class="mord mathnormal">c</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span><span class="mord">/</span><span class="mord mathnormal" style="margin-right:0.03588em;">gv</span><span class="mord mathnormal">c</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span><span class="mord"><span class="mord mathnormal">s</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.1514em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight">c</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal">h</span><span class="mord mathnormal" style="margin-right:0.02778em;">r</span></span></span></span>{i}.vcf";done) \
-O ./7.Germline_Mutations/gvcf/merge.vcf
代码
文本
上面 VCF 文件包含了 germline mutations 信息。
代码
文本
mkdir 8.Mutect2
代码
文本
对给定的肿瘤和正常样本进行突变检测，分析目标区域是指定的外显子区域，并将检测到的 somatic mutation 信息输出到指定的 VCF 文件中：
代码
文本
gatk Mutect2 -R ./0.ref/human_g1k_v37.fasta \
-I ./6.gatk/SRR3182433_bqsr.bam -tumor <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mopen">(</span><span class="mord mathnormal">ba</span><span class="mord mathnormal">se</span><span class="mord mathnormal">nam</span><span class="mord mathnormal">e</span><span class="mord">"./6.</span><span class="mord mathnormal" style="margin-right:0.03588em;">g</span><span class="mord mathnormal">a</span><span class="mord mathnormal">t</span><span class="mord mathnormal" style="margin-right:0.03148em;">k</span><span class="mord">/</span><span class="mord mathnormal" style="margin-right:0.00773em;">SRR</span><span class="mord">318243</span><span class="mord"><span class="mord">3</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3361em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight">b</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal" style="margin-right:0.03588em;">q</span><span class="mord mathnormal" style="margin-right:0.02778em;">sr</span><span class="mord">.</span><span class="mord mathnormal">bam</span><span class="mord"><span class="mord">"</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3361em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight">b</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal" style="margin-right:0.03588em;">q</span><span class="mord mathnormal" style="margin-right:0.02778em;">sr</span><span class="mord">.</span><span class="mord mathnormal">bam</span><span class="mclose">)</span><span class="mspace"> </span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mbin">−</span><span class="mspace" style="margin-right:0.2222em;"></span></span><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mord mathnormal" style="margin-right:0.07847em;">I</span><span class="mord">./6.</span><span class="mord mathnormal" style="margin-right:0.03588em;">g</span><span class="mord mathnormal">a</span><span class="mord mathnormal">t</span><span class="mord mathnormal" style="margin-right:0.03148em;">k</span><span class="mord">/</span><span class="mord mathnormal" style="margin-right:0.00773em;">SRR</span><span class="mord">318242</span><span class="mord"><span class="mord">3</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3361em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight">b</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal" style="margin-right:0.03588em;">q</span><span class="mord mathnormal" style="margin-right:0.02778em;">sr</span><span class="mord">.</span><span class="mord mathnormal">bam</span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mbin">−</span><span class="mspace" style="margin-right:0.2222em;"></span></span><span class="base"><span class="strut" style="height:0.6944em;"></span><span class="mord mathnormal">n</span><span class="mord mathnormal" style="margin-right:0.02778em;">or</span><span class="mord mathnormal">ma</span><span class="mord mathnormal" style="margin-right:0.01968em;">l</span></span></span></span>(basename "./6.gatk/SRR3182423_bqsr.bam" _bqsr.bam) \
-L ./0.ref/exon_data.bed \
-O ./8.Mutect2/SRR3182433_mutect2.vcf
代码
文本
对突变检测结果进行筛选和过滤，以去除那些可能是测序或分析误差引起的假阳性突变，并提高突变检测的可靠性。将经过处理后的突变信息输出到指定的 VCF 文件中：
代码
文本
gatk FilterMutectCalls -R ./0.ref/human_g1k_v37.fasta \
-V 8.Mutect2/SRR3182433_mutect2.vcf \
-O ./8.Mutect2/SRR3182433_somatic.vcf
代码
文本
对 Mutect2 生成的 VCF 文件进行了两个层次的进一步的过滤和处理：首先根据 PASS 标志和染色体名称进行过滤，然后将筛选后的结果保存为一个新的 VCF 文件，以获取更可靠和有用的突变信息。
代码
文本
cat ./8.Mutect2/SRR3182433_somatic.vcf | perl -alne '{if(/^#/){print}else{next unless <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mord mathnormal" style="margin-right:0.13889em;">F</span><span class="mopen">[</span><span class="mord">6</span><span class="mclose">]</span><span class="mord mathnormal">e</span><span class="mord mathnormal" style="margin-right:0.03588em;">q</span><span class="mord">"</span><span class="mord mathnormal" style="margin-right:0.13889em;">P</span><span class="mord mathnormal">A</span><span class="mord mathnormal" style="margin-right:0.05764em;">SS</span><span class="mord">"</span><span class="mpunct">;</span><span class="mspace" style="margin-right:0.1667em;"></span><span class="mord mathnormal">n</span><span class="mord mathnormal">e</span><span class="mord mathnormal">x</span><span class="mord mathnormal">t</span><span class="mord mathnormal">i</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span></span></span></span>F[0] =~/_/;print } }' > ./8.Mutect2/SRR3182433_filter.vcf
代码
文本
gatk Mutect2 -R ./0.ref/human_g1k_v37.fasta \
-I ./6.gatk/SRR3182436_bqsr.bam -tumor <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mopen">(</span><span class="mord mathnormal">ba</span><span class="mord mathnormal">se</span><span class="mord mathnormal">nam</span><span class="mord mathnormal">e</span><span class="mord">"./6.</span><span class="mord mathnormal" style="margin-right:0.03588em;">g</span><span class="mord mathnormal">a</span><span class="mord mathnormal">t</span><span class="mord mathnormal" style="margin-right:0.03148em;">k</span><span class="mord">/</span><span class="mord mathnormal" style="margin-right:0.00773em;">SRR</span><span class="mord">318243</span><span class="mord"><span class="mord">6</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3361em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight">b</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal" style="margin-right:0.03588em;">q</span><span class="mord mathnormal" style="margin-right:0.02778em;">sr</span><span class="mord">.</span><span class="mord mathnormal">bam</span><span class="mord"><span class="mord">"</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3361em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight">b</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal" style="margin-right:0.03588em;">q</span><span class="mord mathnormal" style="margin-right:0.02778em;">sr</span><span class="mord">.</span><span class="mord mathnormal">bam</span><span class="mclose">)</span><span class="mspace"> </span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mbin">−</span><span class="mspace" style="margin-right:0.2222em;"></span></span><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mord mathnormal" style="margin-right:0.07847em;">I</span><span class="mord">./6.</span><span class="mord mathnormal" style="margin-right:0.03588em;">g</span><span class="mord mathnormal">a</span><span class="mord mathnormal">t</span><span class="mord mathnormal" style="margin-right:0.03148em;">k</span><span class="mord">/</span><span class="mord mathnormal" style="margin-right:0.00773em;">SRR</span><span class="mord">318241</span><span class="mord"><span class="mord">9</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3361em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight">b</span></span></span></span><span class="vlist-s"></span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span><span class="mord mathnormal" style="margin-right:0.03588em;">q</span><span class="mord mathnormal" style="margin-right:0.02778em;">sr</span><span class="mord">.</span><span class="mord mathnormal">bam</span><span class="mspace" style="margin-right:0.2222em;"></span><span class="mbin">−</span><span class="mspace" style="margin-right:0.2222em;"></span></span><span class="base"><span class="strut" style="height:0.6944em;"></span><span class="mord mathnormal">n</span><span class="mord mathnormal" style="margin-right:0.02778em;">or</span><span class="mord mathnormal">ma</span><span class="mord mathnormal" style="margin-right:0.01968em;">l</span></span></span></span>(basename "./6.gatk/SRR3182419_bqsr.bam" _bqsr.bam) \
-L ./0.ref/exon_data.bed \
-O ./8.Mutect2/SRR3182436_mutect2.vcf
代码
文本
tools/gatk/gatk FilterMutectCalls -R ./0.ref/human_g1k_v37.fasta \
-V 8.Mutect2/SRR3182436_mutect2.vcf \
-O ./8.Mutect2/SRR3182436_somatic.vcf
代码
文本
cat ./8.Mutect2/SRR3182436_somatic.vcf | perl -alne '{if(/^#/){print}else{next unless <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:1em;vertical-align:-0.25em;"></span><span class="mord mathnormal" style="margin-right:0.13889em;">F</span><span class="mopen">[</span><span class="mord">6</span><span class="mclose">]</span><span class="mord mathnormal">e</span><span class="mord mathnormal" style="margin-right:0.03588em;">q</span><span class="mord">"</span><span class="mord mathnormal" style="margin-right:0.13889em;">P</span><span class="mord mathnormal">A</span><span class="mord mathnormal" style="margin-right:0.05764em;">SS</span><span class="mord">"</span><span class="mpunct">;</span><span class="mspace" style="margin-right:0.1667em;"></span><span class="mord mathnormal">n</span><span class="mord mathnormal">e</span><span class="mord mathnormal">x</span><span class="mord mathnormal">t</span><span class="mord mathnormal">i</span><span class="mord mathnormal" style="margin-right:0.10764em;">f</span></span></span></span>F[0] =~/_/;print } }' > ./8.Mutect2/SRR3182436_filter.vcf
代码
文本
ls -lh 8.Mutect2
代码
文本
上面 VCF 文件是 tumor 样本包含的 somatic mutations 信息。
代码
文本
通过 GATK 的 mutect 流程找到 somatic mutations 只有突变的染色体和坐标信息，需要借助 vcf 注释软件 annovar 快速了解突变发生在哪个基因或者哪号外显子。
代码
文本
cd tools
# 下载 annovar 软件
wget -c http://www.openbioinformatics.org/annovar/download/0wgxR2rIVP/annovar.latest.tar.gz
tar -zxvf annovar.latest.tar.gz
cd annovar
# 下载基因组版本是 hg19 对应的数据库文件，存到文件夹 humandb/ 中
nohup ./annotate_variation.pl -downdb -webfrom annovar gnomad_genome --buildver hg19 humandb/ >down.Log 2>&1 &
ls -lh tools/annovar/humandb/
代码
文本
代码
文本
开始注释：
代码
文本
mkdir ./8.Mutect2/annotation/
代码
文本
tools/annovar/table_annovar.pl ./8.Mutect2/SRR3182436_filter.vcf ./tools/annovar/humandb \
-buildver hg19 \
-out ./8.Mutect2/annotation/SRR3182436 \
-remove \
-protocol refGene \
-operation g \
-nastring . \
-vcfinput
代码
文本
tools/annovar/table_annovar.pl ./8.Mutect2/SRR3182433_filter.vcf ./tools/annovar/humandb/ \
-buildver hg19 \
-out ./8.Mutect2/annotation/SRR3182433 \
-remove \
-protocol refGene \
-operation g \
-nastring . \
-vcfinput
代码
文本
ls -lh 8.Mutect2/annotation/
代码
文本
每个样本输出3个文件：
SRR3182433.hg19_multianno.txt中记录：突变记录位点注释上的信息（下面展示该文件前两行信息）：
代码
文本
mkdir 8.Mutect2/maf
代码
文本
# 以 "Chr" 开头的行移除
grep -v '^Chr' ./8.Mutect2/annotation/SRR3182433.hg19_multianno.txt | cut -f 1-24 | awk -v T=SRR3182433 -v N=SRR3182423_germline '{print $0"\t"T"\t"N}' >./8.Mutect2/maf/SRR3182433.annovar.vcf
grep -v '^Chr' ./8.Mutect2/annotation/SRR3182436.hg19_multianno.txt | cut -f 1-24 | awk -v T=SRR3182436 -v N=SRR3182419_germline '{print $0"\t"T"\t"N}' >./8.Mutect2/maf/SRR3182436.annovar.vcf
# header追加两列 Tumor_Sample_Barcode： 和 Matched_Norm_Sample_Barcode
head -1 ./8.Mutect2/annotation/SRR3182433.hg19_multianno.txt | awk '{print $0 "\tTumor_Sample_Barcode\tMatched_Norm_Sample_Barcode"}' > ./8.Mutect2/maf/header
# 合并 SRR3182433.annovar.vcf 和 SRR3182436.annovar.vcf 成 annovar_merge.vcf
cat ./8.Mutect2/maf/header ./8.Mutect2/maf/*annovar.vcf >./8.Mutect2/maf/annovar_merge.vcf
代码
文本
ls -lh 8.Mutect2/maf
代码
文本
📖 上手指南
上面的注释结果文件 annovar_merge.vcf 读入 R 中（将Kernel切换成 R），转成 maf 格式的文件进行后续操作。
代码
文本
[ ]
install.packages("BiocManager")
library(BiocManager)
BiocManager::install("maftools")
rm(list = ls())
require(maftools)
options(stringsAsFactors = F) # 设置R环境的选项，将字符型变量视为非因子变量，以避免意外的类型转换。
setwd("/Users/sunnan/Downloads/")
## 读入 annovar_merge.vcf 转成 maf 格式变量 annovar.laml
annovar.laml <- annovarToMaf(annovar = "/Users/sunnan/Downloads/annovar_merge.vcf",
refBuild='hg19',
tsbCol = 'Tumor_Sample_Barcode',
table = 'refGene',
MAFobj = T)
save(annovar.laml,file = 'annovar_laml.Rdata')
代码
文本
变量 annovar.laml 如下：
代码
文本
[ ]
## 绘制突变频率汇总图
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）
dev.off()
代码
文本
该突变频率汇总图包括 6 个小图，每个小图解释如下：
Variant Classification：这个小图展示了不同变异分类的数量分布。纵轴显示有 7 种变异分类。
Variant Type：这个小图展示了不同变异类型的数量分布。变异类型可以是单核苷酸变异（SNV）、插入、缺失等。
SNV Class：这个小图展示了单核苷酸变异（SNV）的亚类别数量分布。
Variants per sample (Median: 404)：这个小图展示了 2 个 tumor 样本中的变异数目分布。标题中的 "Median: 404" 表示中位数变异数目为 404，这个值可以用来描述样本间的变异负荷。
Variant Classification summary：这个小图是变异分类的汇总图。它将不同变异分类（例如，高风险、低风险等）的变异数目进行了汇总并展示。
Top 10 mutated genes：这个小图展示了变异频率最高的前 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()
代码
文本
上面结果图片是一个肿瘤学图，标题 "Altered in 2 (100%) of 2 samples." 意味着在所分析的 2 个样本中，这 30 个特定基因在这 2 个 tumor 样本中都发生了变化，发生频率是 100%。横轴表示 tumor 样本，纵轴表示基因，每个方块表示相应的基因在相应样本中是否有突变。图表的颜色和填充模式表示不同的突变类型：
绿色：Missense Mutation（错义突变），这是一种单个核苷酸发生改变，从而导致蛋白质编码中的一个氨基酸发生变化的突变。这可能会影响蛋白质的结构和功能，可能会对蛋白质的正常功能产生影响，甚至导致疾病发生。
红色：Nonsense Mutation（无义突变），这种突变会导致一个编码氨基酸的密码子变成了一个终止密码子，导致蛋白质合成过程提前终止。这通常导致缺失或截短的蛋白质产物，对于正常细胞功能可能产生严重的影响。
蓝色：Frame_Shift_Del（移码缺失），这是指DNA序列中删除了一个或多个碱基，导致蛋白质编码框架发生错位，从而导致编码氨基酸序列发生改变。这通常会导致终止密码子的出现，产生缺失的或截短的蛋白质。
紫色：Frame_Shift_Ins（移码插入），这种变异是指DNA序列中插入了一个或多个碱基，导致蛋白质编码框架发生错位，从而导致编码氨基酸序列发生改变。这同样会导致终止密码子的出现，产生缺失的或截短的蛋白质。
黑色：Multi Hit（多次命中），此类突变可能表示在同一位点发生了多个不同的变异，这可能导致一个或多个碱基发生了改变。这种情况可能对蛋白质功能产生严重影响，因为多次命中会导致复杂的变异模式。
TMB 是 "Tumor Mutational Burden"（肿瘤突变负荷）的缩写。它是用于衡量肿瘤组织中突变数量的指标，通常以每兆碱基对（mutants per megabase，简称 mut/Mb）来表示。高 TMB 意味着肿瘤细胞中存在较多的突变，低 TMB 则表示突变数量较少。TMB 455 意味着在每兆碱基对（Mb）的 DNA 序列中，平均有 455 个突变。
代码
文本
[ ]
## 可视化 SP140 基因突变
postscript('SP140_annovar.eps', width = 10, height = 3, horizontal = FALSE)
maftools::lollipopPlot(maf = annovar.laml,
gene = 'SP140', # 指定要绘制的基因，这里是 "SP140"
AACol = 'AAChange.refGene', # 显示的氨基酸改变信息的列，这里是 "AAChange.refGene"。
labelPos = 'all') # 指定标签显示的位置，这里设置为 "all"，表示显示所有标签。
dev.off()
代码
文本
该 lollipop 图可视化了 SP140 基因的突变。基因 SP140 在所考虑的样本中具有 50% 的体细胞突变率，并提供了该基因的 RefSeq 参考转录本编号 NM_007237。横轴表示基因的各个氨基酸位置，纵轴表示基因在不同样本中的突变信息。每个突变都用一个小的“棒棒糖”（lollipop）符号表示，其中糖棒的一端位于相应的氨基酸位置，另一端指向具体的突变类型，如图中的 Q199H 表示在位置 199，氨基酸由 Q 突变成了 H，突变类型是 Missense Mutation（错义突变）。Sp100、SAND、PHD、Bromo_SP100C_like 是蛋白质结构或结构域的名称。
下面还可以再看 TP53 这个基因的 lollipop 图：
代码
文本
[ ]
## lollipopPlot for TP53
postscript('TP53_annovar.eps', width = 10, height = 3, horizontal = FALSE)
maftools::lollipopPlot(maf = annovar.laml,
gene = 'TP53',
AACol = 'AAChange.refGene',
labelPos = 'all')
dev.off()
代码
文本
代码
文本
通过本文中的详细指导，您已经了解了全外显子测序数据预处理和变异检测的关键步骤。从下载参考基因组、准备测序数据，到质控、比对、过滤和注释，每个环节都为您揭示了如何从原始数据中提取有关遗传变异的宝贵信息。无论您是初学者还是有经验的分析师，这个流程都将为您提供坚实的基础，让您能够更深入地理解全外显子测序数据的内容和意义。随着您逐步掌握每个步骤，您将能够更准确地识别和解释样本中的遗传变异，为生物学研究和临床应用提供有力的支持。无论您的目标是揭示疾病机制、寻找治疗靶点还是进行个性化医疗，这个流程都将引导您朝着正确的方向前进。让我们一同迈出这个分析之旅的第一步，开始探索并挖掘全外显子测序数据所蕴含的丰富信息吧！
代码
文本
已赞3
本文被以下合集收录
生物信息学 Notebooks Collection
liyongge
更新于 2024-09-13
33 篇75 人关注
序列分析
liyongge
更新于 2024-04-01
4 篇1 人关注
推荐阅读
公开
探索全外显子测序数据分析流程：从数据预处理到变异分析副本bohrfb06fb
发布于 2024-05-13
公开
一篇Notebook带你走进单细胞转录组分析liyongge
发布于 2023-06-20
15 赞16 转存文件1 评论