利用二代测序数据完成全基因组变异检测小实例
2020-01-02 本文已影响0人
小明的数据分析笔记本
本篇文章为原文的翻译稿，稍微加上一点自己的理解。
本文使用到的数据：
- 参考基因组：大肠杆菌 NC_008253
- 二代测序数据：单端 wgsim软件模拟生成
操作系统及使用到的软件
- 以下操作基本在linux操作系统完成。分析过程使用到的数据非常小(If you are new to next-generation sequence analysis, you will soon find that one of the biggest obstacles is just finding and downloading sample data sets, and then downloading the very large reference genomes)，使用自己的笔记本电脑完全足够。windows系统需要安装虚拟机，或者win10可以直接安装linux系统。关于win10配置linux系统，可以参考洲更大神的文章「生信基础课」如何利用好手头的电脑，节省上千的服务器租用费
- 使用到的软件
wgsim —— 模拟生成二代测序数据
fastqc —— 检查原始测序数据的质量
fastp —— 原始数据质控及过滤
bowtie2 —— 将二代测序数据比对到参考基因组
qualimap —— 利用bam文件统计测序深度、覆盖度等
samtools —— sam格式BAM格式转化、索引、排序等
bcftools —— vcf文件操作
关于软件的安装方法同样可以参考这篇文章的内容「生信基础课」如何利用好手头的电脑，节省上千的服务器租用费
利用二代测序数据分析全基因组变异的基本流程
- 1、样本DNA提取
- 2、DNA测序
- 3、原始测序数据质量检查、过滤 （fastp）
- 4、高质量测序数据比对到参考基因组 (bowtie2)
- 5、比对数据评估及可视化 (samtools)
-
6、变异检测（SNP——single nucleotide polymorphisms; small insertions and deletions；但核苷酸多态性，插入缺失等）(samtools)
图片来源：https://www.ebi.ac.uk/sites/ebi.ac.uk/files/content.ebi.ac.uk/materials/2014/140217_AgriOmics/dan_bolser_snp_calling.pdf
本篇文章用到的数据下载
You can access all sample data files for this primer from the https://github.com/ecerami/samtools_primer
分析过程
本篇文章的分析过程只用到Ecoli参考基因组的前1000个碱基（参考基因组>后面的内容太长，我将其改为NC_008253）
- 使用samtools获取参考基因组任意位置的碱基
此处参考GATK4.0和全基因组数据分析实践（上）
samtools faidx NC_008253.fna
samtools faidx NC_008253.fna NC_008253:1-1000 > NC_008253_1K.fna
第一步：模拟生成单端测序数据
wgsim -N 1000 -S 1 NC_008253_1K.fna sim_reads.fastq /dev/null
- 参数-N指定生成参考序列的数量，本文用到1000条
- 参数-S指定随机数种子，这里设置为1，则生成的数据和原文一模一样
- 接下来是参考基因组
- 接下来是双端测序数据的文件名，如果生成单端，则第二个文件名为
/dev/null
运行以上命令的结果为
[wgsim] seed = 1
[wgsim_core] calculating the total length of the reference sequence...
[wgsim_core] 1 sequences, total length: 1000
NC_008253:1-1000 736 T G -
结果表示生成的测序数据中有一个SNP，在736这个位置，参考（ref）位置碱基为T，样本中为G
模拟生成双端数据
wgsim -N 4000 -1 150 -2 150 NC_008253.fna reads_1.fastq reads_2.fastq
第二步：检查测序数据的质量
fastqc sim_reads.fastq
关于fastqc的结果解读可以参考从零开始完整学习全基因组测序（WGS）数据分析：第3节 数据质控
第三步：根据一定的标准对原始数据进行过滤
本篇文章省略这一步
第四步：原始数据比对到参考基因组
# 构建索引
bowtie2-build NC_008253_1K.fna NC_008253_1K
# 比对
bowtie2 -x NC_008253_1K -U sim_reads.fastq -S sim_reads_aligned.sam
## 输出结果
1000 reads; of these:
1000 (100.00%) were unpaired; of these:
34 (3.40%) aligned 0 times
966 (96.60%) aligned exactly 1 time
0 (0.00%) aligned >1 times
96.60% overall alignment rate
SMA、BAM文件的格式原文介绍非常详细，加之自己也有不太理解的地方，这里不详细介绍
第五步：SAM格式转换BAM格式、排序、索引等
# SAM 转 BAM
samtools view -b -S -o sim_reads_aligned.bam sim_reads_aligned.sam
# 查看BAM文件
samtools view sim_reads_aligned.bam | less -S
samtools view sim_reads_aligned.bam | less -S | wc -l
# 根据一定的标准展示结果
samtools view -f 4 alignments/sim_reads_aligned.bam | less -S
samtools view -c -f 4 sim_reads_aligned.bam
samtools view -q 42 sim_reads_aligned.bam
#排序
samtools sort sim_reads_aligned.bam -o sim_reads_aligned.sorted.bam
#索引
samtools index sim_reads_aligned.sorted.bam
- 参数-f匹配指定标准（SAM flag）的测序数据，4指未能比对到参考基因组的reads的数量，
- 参数-F删除指定标准的测序数据
- 参数-c展示匹配指定标准的reads的数量
samtools view -c -f 4 sim_reads_aligned.bam输出结果为34，与bowtie2比对结果展示的数据相符合。 - 参数-q指定最小比对质量值，
samtools view -c -q 42 sim_reads_aligned.bam比对质量值大于42的有819条 - 排序的步骤比原文多了一个-o参数，因为这里使用的samtools版本和原文不一样
第六步：利用bam文件统计比对后的一些数据指标
这一步因为bam文件非常小，可以直接使用qualimap软件的windows版本
使用方法File——New Analysis——BAM QC，然后选择bam文件即可。关于qualimap软件linux版本的使用方法还需要继续摸索。
输出结果
image.png
1.png
2.png
3.png
第七步：变异检测
samtools mpileup -g -f NC_008253_1K.fna sim_reads_aligned.sorted.bam > sim_variants.bcf
### 这里提示
[warning] samtools mpileup option `g` is functional, but deprecated. Please switch to using bcftools mpileup in future.
意思是samtools mpileup option `g` 可能不在使用了，建议转为使用bcftools mpileup
bcftools call -c -v sim_variants.bcf > sim_variants.vcf
###这里提示
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
关于vcf文件的格式和内容自己还有很多不理解的地方，还需要多学习，可以参考GATK4.0和全基因组数据分析实践（下），「博客翻译」SNP过滤教程（一）
图片来源于原文的截图
第八步：可视化bam文件
samtools tview sim_reads_aligned.sorted.bam NC_008253_1K.fna
image.png
通过h，j，k，l，空格可以上下左右移动，按下字母g可以移动到指定位置
image.png
下一步计划
利用wgsim模拟生成多组数据，然后学习公众号宇宙实验媛分享的群体分化指数—Fst等与群体遗传学有关的文章。
