全外显子测序(WES)由入门到精通记录

背景

作者是一个没有生信背景知识的程序员,近期对全外显子测序有些感兴趣,因此记录一下学习过程,和大家共同进步。

预备知识

WES基础知识
什么是基因
什么是DNA
二代测序基础知识
等位基因

软件安装

bwa

git clone https://github.com/lh3/bwa.git
cd bwa
make

由于bwa的makefile没有写intall的脚本,所以我们直接将编译的目录加入到环境变量中,

vim ~/.bashrc
export PATH=${PATH}:/usr/local/cuda/bin/:/home/fangl5/Downloads/bwa-0.7.17/

samtools

下载samtools-1.14.tar.bz2
解压后

    ./configure
    make
    sudo make install

java 11

sudo apt-get install openjdk-11-jdk

picard

wget https://github.com/broadinstitute/picard/releases/download/2.26.10/picard.jar

gatk

sudo apt-get install git-lfs
git clone https://github.com/broadinstitute/gatk.git
./gradlew bundle

将gatk加入到环境变量PATH中,和bwa中操作一致。

SOAPnuke

git clone https://github.com/BGI-flexlab/SOAPnuke.git
cd SOAPnuke 
make

加入环境变量中

下载数据

hg38.ga

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz

gatk参考文件

我直接将google cloud中文件夹下载下来了

gsutil -m cp -r \
  "gs://genomics-public-data/resources/broad/hg38/v0/" \
  .
gatk-tutotial 相关数据 google cloud地址
genomics-public-data
annovar 数据下载

http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/

实践fastq to vcf

# 大于2G的参考基因组要用参数-a bwtsw  
# 约40min
bwa index -a bwtsw hg38.fa
java -jar ~/Downloads/picard.jar CreateSequenceDictionary R=./hg38.fa O=./hg38.dict
samtools faidx hg38.fa
gatk IndexFeatureFile -I ./Mills_and_1000G_gold_standard.indels.hg38.vcf

https://app.terra.bio/#workspaces/help-gatk/GATK4-Germline-Preprocessing-VariantCalling-JointCalling/notebooks/launch/1-gatk-germline-variant-discovery-tutorial.ipynb

Q&A

核苷酸序列或氨基酸序列是什么?
WGS\WES\TS区别

http://www.genomesop.com/wp-content/uploads/2018/11/genomic_comparison3.png
WGS\WES\TS区别

参考资料

GATK4 流程分析- 从fastq到vcf
GATK_Discovery_Tutorial-Worksheet-AUS2016.pdf
基因组测序、外显子测序和靶向测序有什么样的区别,如何选择?
华大外显子测序
Exome sequencing data analysis for diagnosing a genetic disease
全外显子测序(WES)分析流程
使用fastp进行数据质控
GATK4 多个样本GenotypeGVCFs前用 CombineGVCFs还是GenomicsDBImport
intel 2015 面向 GATK 最佳实践流水线部署的基础设施
家系分析软件汇总
WGS数据分析入门篇 碱基矿工
全外显子测序(WES)分析流程 这一篇,用到了bed文件。
gatk spark加速

Whole-Genome-Sequencing
全外显子测序(wes)数据分析详细流程(小样本)
如何正确设置GATK VQSR的模型训练参数
https://blog.csdn.net/weixin_43217642/article/details/115291586

<think>好的,我现在要回答用户关于如何从外显子测序数据中得到VAF是否大于2%的问题。首先,我需要确认自己了解外显子测序WES)的基本流程和VAF的定义。VAF是变异等位基因频率,即某个变异位点上突变型等位基因占所有等位基因的比例。用户想知道如何从测序数据中计算这个值,并判断是否超过2%。 接下来,我应该分步骤思考如何处理数据。首先,原始数据通常是FASTQ格式,需要经过质控和比对。这部分用户可能已经处理过,但需要确认。然后,比对后的BAM文件需要进行变异检测,生成VCF文件。VCF文件里包含每个变异位点的信息,比如染色体位置、参考碱基、变异碱基等。这时候就需要关注VAF的计算了。 接下来,我需要解释VCF文件中的字段。在VCF的INFO或FORMAT部分,可能有AD(等位基因深度)和DP(总深度)字段。AD表示每个等位基因的测序深度,比如参考等位和变异等位。例如,AD=30,5表示参考等位有30个reads支持,变异等位有5个。这时候VAF的计算就是变异AD除以总DP,即5/(30+5)=14.29%。这时候用户需要检查VCF中是否已经有VAF直接给出,或者需要手动计算。 如果VCF中没有直接给出VAF,用户可能需要自己计算。这时候需要解析AD和DP的值,或者AF字段。有的流程可能会用AF来直接表示等位基因频率,但需要确认是否已经乘以100变成百分比。比如AF=0.02表示2%,这时候直接比较即可。如果AF是0.02,判断是否大于0.02,即2%。 然后,用户可能需要使用工具来过滤VCF,比如bcftools或者SnpSift。例如,用bcftools filter命令设置条件来筛选VAF大于2%的位点。这需要正确的命令行参数,比如-i 'FORMAT/AF>0.02'。同时,需要注意字段名称可能不同,比如有的用VAF,有的用AF,或者需要计算。 此外,数据质量的问题也需要考虑。比如低深度区域可能导致VAF估计不准确,假阳性变异等。这时候需要提醒用户注意质控步骤,确保DP足够高,比如至少20x,但具体阈值可能因研究而异。 最后,总结步骤:数据预处理、变异检测、提取或计算VAF、过滤并判断是否超过2%。同时,给出示例命令和注意事项,帮助用户实际操作。 在解释过程中,要确保使用正确的数学格式,比如VAF的计算公式用$$...$$,行内变量用$...$。同时保持结构清晰,分步骤说明,并强调可能遇到的问题,如字段名称差异、数据质量等。</think>外显子测序WES)数据中计算变异等位基因频率(Variant Allele Frequency, VAF)并判断是否大于2%,需经过以下步骤: --- ### **1. 原始数据处理与变异检测** - **测序数据质控**:原始数据(FASTQ文件)需通过工具(如FastQC)进行质量评估,过滤低质量reads。 - **序列比对**:使用BWA、Bowtie2等工具将reads比对到参考基因组,生成BAM文件。 - **变异检测**:通过GATK、VarScan2等工具识别单核苷酸变异(SNV)和小片段插入缺失(Indel),输出VCF文件。 --- ### **2. 从VCF文件中提取VAF信息** VCF文件(Variant Call Format)包含变异位点的详细信息,关键字段包括: - **`AD`(Allelic Depth)**:等位基因的测序深度(例如:`AD=30,5`表示参考等位支持30个reads,变异等位支持5个reads)。 - **`DP`(Total Depth)**:该位点的总测序深度(例如:`DP=35`)。 - **`AF`(Allele Frequency)**:部分工具直接提供变异等位频率(例如:`AF=0.1429`表示14.29%)。 #### **VAF计算公式**: $$ \text{VAF} = \frac{\text{变异等位基因的深度(AD\_alt)}}{\text{总深度(DP)}} \times 100\% $$ **示例**:若`AD=30,5`且`DP=35`,则: $$ \text{VAF} = \frac{5}{30+5} \times 100\% = 14.29\% $$ --- ### **3. 判断VAF是否大于2%** #### **方法1:直接解析VCF字段** - 若VCF中已包含`AF`字段,直接提取该值并判断是否大于0.02(即2%)。 - **示例命令(使用bcftools)**: ```bash bcftools view -i 'FORMAT/AF > 0.02' input.vcf > filtered.vcf ``` #### **方法2:手动计算VAF** - 若VCF无`AF`字段,需从`AD`和`DP`手动计算: ```bash # 使用bcftools提取AD和DP,并用awk计算VAF bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%AD\t%DP\n' input.vcf | \ awk -F'\t' '{ split($5, ad, ","); vaf = ad[2] / $6; if (vaf > 0.02) print $1,$2,$3,$4,vaf }' ``` --- ### **4. 注意事项** 1. **数据质量**:低深度(DP < 20)可能导致VAF估算不准确,需过滤低质量位点。 2. **肿瘤样本**:若为肿瘤-正常配对样本,需扣除正常样本的胚系变异(如使用MuTect2)。 3. **工具差异**:不同工具生成的VCF字段可能不同(如VarScan2使用`FREQ`字段)。 --- ### **总结步骤** 1. 从VCF中提取`AD`和`DP`(或直接使用`AF`字段)。 2. 计算VAF或直接读取现有值。 3. 通过脚本或工具筛选VAF > 2%的变异位点。 通过以上流程,可高效判断外显子测序数据中变异的VAF是否满足阈值要求。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值