技术介绍

  • 第二代测序技术，1990s - 2010s

第二代测序技术从原理上分为三种方法：

  • Roche 454测序法
  • IIIumina Solexa/Hiseq测序法
  • ABI Solid测序法

三种二代测序技术的对比：

image.png

总结

二代测序技术总体而言都有着同一种特性：

  • 需荧光或化学发光物质
  • 需聚合酶或连接酶
  • 需购买昂贵的试剂耗材和光学系统
  • 需强大的图形分析计算能力

WES分析流程

WES分析流程分为5个步骤：

  1. WES分析流程
  2. 数据质量控制
  3. 序列比对分析
  4. 变异检测
  5. 变异注释

一. WES分析流程

原始测序数据——数据质量控制——序列比对——变异检测——变异注释——完成

二. 数据质量控制

数据质控标准：

  1. 去掉reads中的接头
  2. 去除低质量（BP<20）碱基的reads
  3. 去除序列头尾的N碱基
  4. 去除头尾N碱基后若剩余reads长度小于40bp（双端），则丢弃该对序列

原始数据下载——sratoolkit软件

1
2
3
4
5
6
# 安装
conda install sra-tools
# 找到sra数据，下载srr list
prefetch SRR1139956
# sra转为fastq
nohup fastq-dump -gzip --split-3 -O ./ SRR1139956.sra &

基本参数：

fastq-dump [options] [path]

  1. –gzip: 使用gzip压缩输出文件
  2. –split-3: 设置文件拆分，如果是双端的序列，则拆分为*_1.fastq和*_2.fastq
  3. -O: 设置输出文件路径

数据质量控制常用的软件是FASTQC，在这里不再赘述，详情请看https://www.jinhenghaoblog.top/posts/20250828165237.html

序列比对分析

进行序列比对分析所用到的软件为BWA：

1
2
3
4
5
6
7
8
# 下载参考数据
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz
gunzip Homo_sapiens_assembly38.fasta.gz
# 建立索引
bwa index Homo_sapiens_assembly38.fasta
# bwa进行比对
nohup bwa mem -t 5 db/Homo_sapiens_assembly38.fasta test.fastq 1>test.sam 2>test.log
# 这里的test.fastq是经过过滤得到的fastq文件

基本参数：

bwa mem [options] [idxbase] [in1.fq] [in2.fq]

  1. idxbase: 指定比对参考基因组序列（建好索引）
  2. in1.fq、in2.fq：输入fastq文件
  3. -t : 设置线程，这里设置的是5，具体看计算机配置

SAM文件格式转换及排序

1
2
3
4
5
# sam转bam
samtool view -bS -h test.sam > test.bam

# 排序bam
samtools sort -@ 5 test.bam -o test.sort.bam

基本参数：

samtools sort [options] [in.bam] [out.prefix]

  1. -n: 根据reads名字排序
  2. -o: 输出文件
  3. -@: 设置线程，这里设置的是5，具体看计算机配置

标记/删除PCR重复的reads——GATK

1
2
3
4
5
# 建立输出文件夹
mkdir gatk_markdup
cd gatk_markdup
# 标记PCR重复的reads
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates -I test.sort.bam -M test.metrics 1>test_log.mark 2 >&1

基本参数：

gatk MarkDuplicates -I [INPUT] -O [OUTPUT] -M [METRICS]

  1. MarkDuplicates: 鉴定重复的reads，是GATK中的一个工具
  2. -I ： 输入文件
  3. -O ： 输出文件
  4. -M： 复制度量写入的文件

变异检测

变异检测——标记flagstat

1
2
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" FixMateInformation -I test_marked.bam -O test_marked_fixed.bam -SO coordinate 1>test_log.fix 2>&1
samtools index test_marked_fixed.bam

基本参数：

gatk FixMateInformation -I [INPUT] -O [OUTPUT] -M [METRICS]

  1. FixMateInformation : 是GATK中的一个工具
  2. -I：输入文件
  3. -O：输出文件
  4. -SO： 输出文件的排序方式，这里选择的是coordinate

变异检测——GATK找变异

1
2
3
4
5
6
7
8
9
10
# 找变异
# 下载参考
nohup wget ftp://gsapubftp-anonymous@ftp,broadinstitute.org/bundle/hg38/dbsnp 146.hg38.vcf.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills and 1000G gold standard.indels.hg38.vcf.gz &

indel="Mills_and_1000G_gold_standard.indels.hg38.vcf.gz"
snp="dbsnp_146.hg38.vcf.gz"

# 运行命令
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" BaseRecalibrator -R /analysis/Database/hg38/hg38.fa -I test_marked_fixed.bam --known-sites $snp --known-sites $indel -O test_recal.table

基本参数：

gatk BaseRecalibrator -I [INPUT] -O [OUTPUT] -M [METRICS]

  1. –known-sites: 已知的多态位点数据库（可以指定多个）
  2. -R： 参考序列文件
  3. -I：输入文件
  4. -O：输出文件

变异检测——矫正bam并生成VCF

1
2
3
4
5
6
# ————————矫正bam
gatk --java-options "-XmX20G -Djava.io.tmpdir=./" ApplyBQSR -R /analysis/Database/hg38/hg38.fa -I test_marked_fixed.bam -bqsr test_recal_table -O test_bqsr.bam 1>test_log.ApplyBQSR 2>&1

# ——————calling variation得到vcf
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller -R /analysis/Database/hg38/hg38.fa -I test_bqsr.bam --dbsnp $snp -O test_raw.vcf 1>test_log.HC 2>&1
samtools index test.bqsr.bam

基本参数：

gatk HaplotypeCaller -i [INPUT] -O [OUTPUT] -M [METRICS]

  1. -R：参考序列文件
  2. -I：输入文件
  3. –dbsnp：dbSNP文件
  4. -O：输出文件

变异注释

变异注释使用的软件为Annovar软件

Annovar软件可实现三种不同的注释方式：

  • 确定SNP或CNV是否导致蛋白质编码变化和确定受影响的氨基酸
  • 识别特定基因组区域的变异
  • 鉴定特定数据库中记录的变异

Annovar——下载数据库

1
perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar refGene humandb/

基本参数：

Perl annotate_variation.pl [options]

  1. -buildver: 基因组对应版本，这里是hg19
  2. -webfrom annovar：指定从annovar库里下载，若annovar库中没有，则会从UCSC中下载
  3. refGene：数据库名称
  4. humandb/：数据库存储路径

Annovar——变异注释

1
2
3
4
# 输入文件格式转换
perl conver2annovar.pl -format vcf4 -allsample test.indel.vcf.gz -outfile anno
# 注释
perl annotate_variation.pl -generanno -buildver hg38 test.avinput humandb/

基本参数：

Perl table_annovar.pl [options]

  1. -remove: 删除所有的临时文件
  2. -geneanno：指定注释方法。这里表示gene-based，-regionanno表示region-based，–filer表示filter-based；
  3. -buildver: 指定参考基因组