BSA分析（五）变异检测

❝

「因业务拓展，想组建一个数据分析团队（目前已有RNA-Seq、Chip-Seq、重测序与群体遗传、基因家族、比较基因组、宏基因组、微生物多样性16s/18s/ITS方向专业人员），欢迎有各种数据分析基础的朋友加入我们！——Bioinfor 生信云」

❞

基于二代测序数据进行基因组水平的 SNP、INDEL 变异检测。

二代 DNA 测序数据变异检测

构建基因组index

重测序分析很多软件不能直接读取基因组，需要对基因组提前构建好 index

• 软件：BWA、samtools
• 参考基因组：genome.fasta

# bwa 构建 index
bwa index genome.fasta

# 生成.fai 和 .dict文件, GATK软件依赖
 samtools faidx genome.fasta
samtools dict -o genome.dict genome.fasta

比对

• 软件：BWA、samtools
• 文件准备 参考基因组：genome.fasta 双端测序数据文件: S1_1.fq.gz 和 S1_2.fq.gz

bwa mem \
-t 2 \ # 线程数
-R '@RG\tID:S1\tSM:S1\tPL:illumina' \ # 添加read group信息
../index/genome.fasta \ # 参考基因组index路径
../S1_1.fq.gz \ # fq1数据路径
../S1_2.fq.gz \ # fq2数据路径
2>S1.bwa.log | \
samtools fixmate -m - - | \ # 添加 ms 标记
samtools sort -@ 2 \ # 线程数
-m 1G | \ # 每个线程内存分配
 samtools markdup - S1.sort.markdup.bam # 去除重复

# 构建bam index
samtools index S1.sort.markdup.bam

# 比对率统计
samtools flagstat S1.sort.markdup.bam > S1.sort.markdup.bam.flagstat

# 覆盖度统计
samtools coverage S1.sort.markdup.bam > S1.sort.markdup.bam.coverage

GATK 变异检测

GATK是 Broad 开发的用于测序数据的变异检测软件，后续推广到动植物研究中，是目前最广泛使用的变异检测软件。GATK有两种方式变异检测，一种是合并所有样品的gVCF文件，再通过GenotypeGVCF做变异检测;另一种是每个样品的GVCF文件先生成genomeDB文件再进行变异检测。两者的结果没有差异，根据自己的样品数量合理选择变异检测的方法。

• 软件：GATK
• 文件准备 参考基因组：genome.fasta 比对结果文件：S1.sort.markdup.bam

#变异检测
## 首先创建一个存放临时文件的目录
mkdir tmp

## 每个样品进行HaplotypeCaller变异检测
gatk --java-options "-Xmx10g -Djava.io.tmpdir=./tmp" \  #设置java参数
HaplotypeCaller \  
-R ./genome.fasta \  #参考基因组路径
-I ./S1.sort.markdup.bam  \ #bam文件路径
-ERC GVCF \ # 输出GVCF文件
-O S1.g.vcf  \ # 指定输出结果
1>S1.HC.log   2>&1  #写到日志文件

## 合并每个样品的GVCF文件
ls ./*.g.vcf > gvcf.list  \ #合并文件的列表
gatk  --java-options "-Xmx10g \ #限制使用内存
-Djava.io.tmpdir=./tmp"   \
CombineGVCFs \
-R ./genome.fasta \ 参考基因组
-V gvcf.list  \ #GVCF文件列表
-O all.merge.g.vcf  #输出的vcf文件

# 提取SNP。SNP(单核苷酸多态性)主要是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态性，包括单个碱基的转换、颠换等，为了定位目标性状，每组性状相关样本一起call 群体SNP。


## 群体变异检测
gatk  --java-options "-Xmx10g -Djava.io.tmpdir=./tmp"   \
GenotypeGVCFs -R ./genome.fasta \
-v all.merge.g.vcf \ #群体的GVCF文件
-O all.merge_raw.vcf  #输出文件

## 提取SNP
gatk  --java-options "-Xmx10g -Djava.io.tmpdir=./tmp"  \
SelectVariants  \
-R ./genome.fasta \
-V all.merge_raw.vcf \#合并后的VCF文件
--select-type SNP \#指定提取的类型SNP
-O all.raw.snp.vcf 

## 过滤SNP
gatk  --java-options "-Xmx10g -Djava.io.tmpdir=./tmp"  \
VariantFiltration \
-R ./genome.fasta \
-V all.raw.snp.vcf \
--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name 'SNP_filter' \ #符合其中任意一个条件的都过滤掉
-O all.filter.snp.vcf

## 提取过滤好的SNP
gatk  --java-options "-Xmx4g -Djava.io.tmpdir=./tmp"  \
SelectVariants  -R ./genome.fasta \
-V all.filter.snp.vcf --exclude-filtered \
 -O all.filtered.snp.vcf

# 提取INDEL。InDel是小型的Insertion和Deletion的总称，利用SelectVariants命令来提取 InDel 突变信息。


#提取InDel文件
gatk  --java-options "-Xmx10g -Djava.io.tmpdir=./tmp"  \
SelectVariants  -R ./genome.fasta \#参考基因组
-V all.merge_raw.vcf \#合并后的VCF文件
--select-type INDEL \#筛选提取类型
-O all.raw.indel.vcf  #输出文件

#过滤InDel
gatk  --java-options "-Xmx10g -Djava.io.tmpdir=./tmp"  \
VariantFiltration -R ./genome.fasta  \
-V all.raw.indel.vcf  \ #提取后的InDel文件
--filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 ||  ReadPosRankSum < -8.0"  \ #符合其中任意一个条件的都过滤掉
--filter-name 'INDEL_filter'  \ 
-O all.filter.indel.vcf

#提取过滤后的InDel
gatk  --java-options "-Xmx10g -Djava.io.tmpdir=./tmp"  \
SelectVariants  -R ./genome.fasta   \
-V all.filter.indel.vcf   \#过滤后的InDel文件
--exclude-filtered   \
-O all.filtered.indel.vcf

欢迎加入学习交流群