前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >大肠杆菌全基因组重测序变异检测小实例(侧重变异过滤)

大肠杆菌全基因组重测序变异检测小实例(侧重变异过滤)

作者头像
用户7010445
发布2020-03-03 15:07:44
1.7K0
发布2020-03-03 15:07:44
举报

本文偏重对vcf文件的探索以及设置过滤标准

原文地址

Filtering and handling VCFs

fastq测序获取数据

未找到原文所用数据,本文使用GATK4.0和全基因组数据分析实践(上)文章中的大肠杆菌基因组作为参考序列,使用wgsim软件模拟生成双端150bp测序数据

代码语言:javascript
复制
wgsim -N 80000 -1 150 -2 150 ../Reference_genome/ecoli.fa sim_1_reads_R1.fastq sim_1_reads_R2.fastq

wgsim -N 80000 -1 150 -2 150 ../Reference_genome/ecoli.fa sim_2_reads_R1.fastq sim_2_reads_R2.fastq

wgsim -N 80000 -1 150 -2 150 ../Reference_genome/ecoli.fa sim_3_reads_R1.fastq sim_3_reads_R2.fastq

-N指定生成reads的条数 -1 -2生成reads的长度 接下来是参考序列 接下来是fastq文件的名字

使用samtools变异检测获取vcf文件

这一部分参考文章

  • GATK4.0和全基因组数据分析实践(上)
  • Variant calling tutorial

基本流程: bwa比对 samtools变异检测 完整代码

代码语言:javascript
复制
###构建索引
bwa index Reference_genome/ecoli.fa
bwa mem -t 4 -R '@RG\tID:foo\tSM:sample1' Reference_genome/ecoli.fa sim_reads/sim_1_reads_R1.fastq sim_reads/sim_1_reads_R2.fastq -o output_results/sim_1_aligned.sam
bwa mem -t 4 -R '@RG\tID:foo\tSM:sample2' Reference_genome/ecoli.fa sim_reads/sim_2_reads_R1.fastq sim_reads/sim_2_reads_R2.fastq -o output_results/sim_2_aligned.sam
bwa mem -t 4 -R '@RG\tID:foo\tSM:sample3' Reference_genome/ecoli.fa sim_reads/sim_3_reads_R1.fastq sim_reads/sim_3_reads_R2.fastq -o output_results/sim_3_aligned.sam

这里遇到的问题:

  • -R参数后面接的内容都是什么意思?
  • 学着在shell下写循环
代码语言:javascript
复制
cd output_results
#SAM装换为BAM
samtools view -S -b -o sim_1_aligned.bam sim_1_aligned.sam 
samtools view -S -b -o sim_2_aligned.bam sim_2_aligned.sam 
samtools view -S -b -o sim_3_aligned.bam sim_3_aligned.sam 
#排序
samtools sort sim_1_aligned.bam -o sim_1_aligned.sorted.bam
samtools sort sim_2_aligned.bam -o sim_2_aligned.sorted.bam
samtools sort sim_3_aligned.bam -o sim_3_aligned.sorted.bam
#索引
samtools index sim_1_aligned.sorted.bam 
samtools index sim_2_aligned.sorted.bam 
samtools index sim_3_aligned.sorted.bam
#变异检测
time samtools mpileup -g -t DP,AD -f ../Reference_genome/ecoli.fa sim_1_aligned.sorted.bam sim_2_aligned.sorted.bam sim_3_aligned.sorted.bam > sim_variants_3sample.bcf

###其一
time bcftools call -v -c sim_variants_3sample.bcf > sim_variants_3sample.vcf
###其二
time bcftools call -f GQ,GP -vmO z sim_variants_3sample.bcf -o sim_variants_3sample_1.vcf.gz

这样就得到了最终的vcf格式的文件。

这里遇到的问题:samtools加上bcftools检测变异的各个参数的含义还不太明白!

接下来重复原文内容

查看vcf文件中检测到多少没有经过过滤的变异
代码语言:javascript
复制
bcftools view -H sim_variants_3sample.vcf | wc -l
6918

通常获得的vcf文件都比较大,可以通过随机取样的方法获得小的vcf文件用于后续的分析

过滤vcf文件通常考虑四点:

  • Depth 深度(最小深度和最大深度)
  • Quality 质量值(>30)
  • Minor allele frequency 最小等位基因频率(MAF)
  • Missing data 缺失数据(如何过滤缺失数据需要具体情况具体分析,但是位点缺失率大于25%应该被舍弃)
计算等位基因频率
代码语言:javascript
复制
cd ../
mkdir vcf_handling
cd vcf_handling
vcftools --vcf ../output_results/sim_variants_3sample.vcf --freq2 --out sim_variant_AF
计算每个个体的平均深度
代码语言:javascript
复制
vcftools --vcf ../output_results/sim_variants_3sample.vcf  --depth --out sim_variant_depth
计算每个变异位点的平均深度
代码语言:javascript
复制
vcftools --vcf ../output_results/sim_variants_3sample.vcf --site-mean-depth --out sim_variant_sitedepth
计算位点平均质量值
代码语言:javascript
复制
vcftools --vcf ../output_results/sim_variants_3sample.vcf --site-quality --out sim_variant_sitequality
计算每个个体数据缺失率
代码语言:javascript
复制
vcftools --vcf ../output_results/sim_variants_3sample.vcf --missing-indv --out sim_variant_missingindv

将以上获得的五个文件导出到本地电脑,利用R语言进行可视化展示

五个文件分别是

代码语言:javascript
复制
sim_variant_missingindv.imiss
sim_variant_sitequality.lqual
sim_variant_sitedepth.ldepth.mean
sim_variant_depth.idepth
sim_variant_AF.frq
查看位点质量值分布
代码语言:javascript
复制
setwd("../../vcf_handling/")
library(tidyverse)
var_qual<-read_delim("sim_variant_sitequality.lqual",
                     delim="\t",
                     col_names=c("chr","pos","qual"),
                     skip=1)
library(ggplot2)
a<-ggplot(var_qual,aes(qual))+
  geom_density(fill="dodgerblue1")+
  theme_light()
a1<-a+xlim(0,100)
library(ggpubr)
ggarrange(a,a1,ncol=1,nrow=2)

image.png 从上图可以看出我们的位点质量值是偏低的,因为数据量比较小,位点质量值30代表检测出来的变异有千分之一的可能是错误的,推荐过滤变异的时候设置位点质量值大于30。

变异位点的平均深度
代码语言:javascript
复制
summary(var_depth$mean_depth)
var_depth<-read_delim("sim_variant_sitedepth.ldepth.mean",
                      delim="\t",col_names=c("chr","pos","mean_depth","var_depth"),
                      skip=1)
a<-ggplot(var_depth,aes(mean_depth))+
  geom_density(fill="dodgerblue1",color="black",alpha=0.3)+
  theme_light()
a

image.png 覆盖度的最小值设置(We could set our minimum coverage at the 5 and 95% quantiles这句话暂时还没看懂是什么意思!) 覆盖度最大值的设置:推荐设置为平均深度的两倍(Usually a good rule of thumb is something the mean depth * 2, so in this case we could set our maximum depth at 40x)

平均缺失率
代码语言:javascript
复制
INDV    N_DATA  N_GENOTYPES_FILTERED    N_MISS  F_MISS
sample1 6918    0       0       0
sample2 6918    0       0       0
sample3 6918    0       0       0

本次分析中的缺失率都为0,暂时还没想到是什么原因

最小等位基因频率
代码语言:javascript
复制
var_freq<-read_delim("sim_variant_AF.frq",delim="\t",
                     col_names=c("chr","pos","nalleles","nchr","a1","a2"),skip=1)
var_freq
var_freq$maf<-var_freq%>%
  select(a1,a2)%>%
  apply(1,function(z) min(z))
a<-ggplot(var_freq,aes(maf))+
  geom_density(fill="dodgerblue1",alpha=0.3)+
  theme_light()
a

image.png 这部分的解释自己还没有太看懂,留待后续分解

根据位点质量值和测序深度过滤我们的vcf文件

代码语言:javascript
复制
vcftools --vcf ../output_results/sim_variants_3sample.vcf --minQ 30 --min-meanDP 4 --max-meanDP 10 --minDP 4 --maxDP 10 --recode --stdout | gzip -c > sim_variants_filtered.vcf

参数含义: --vcf or --gzvcf输入未过滤的vcf文件 --minQ位点质量值 --min-meanDP位点最小平均深度 --max-meanDP位点最大平均深度 minDP the minimum depth allowed for a genotype maxDPthe maximum depth allowed for a genotype

剩下多少个变异
代码语言:javascript
复制
bcftools view -H sim_variants_filtered.vcf.gz | wc -l
2212

过滤掉了将近三分之二

这样一次完整的变异检测流程就完成了,当然这其中还有好多细节部分的知识点需要学习!

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2020-01-09,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 小明的数据分析笔记本 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 本文偏重对vcf文件的探索以及设置过滤标准
    • 原文地址
      • fastq测序获取数据
        • 使用samtools变异检测获取vcf文件
          • 接下来重复原文内容
            • 将以上获得的五个文件导出到本地电脑,利用R语言进行可视化展示
              • 根据位点质量值和测序深度过滤我们的vcf文件
              领券
              问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档