测序数据比对和变异检测

全基因组重测序是对已知基因组序列的物种进行不同个体的基因组测序，并在此基础上对个体或群体进行差异性分析。所以首先我们需要某个物种的全基因组序列和该物种的某个个体的基因组序列。

重测序数据分析流程

重测序数据分析流程

一 、参考基因组的下载

进入NCBI下载ASM1252v1的参考基因组

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/012/525/GCA_000012525.1_ASM1252v1/GCA_000012525.1_ASM1252v1_genomic.fna.gz

gunzip GCA_000012525.1_ASM1252v1_genomic.fna.gz

下载地址的获取

二、测试序列（即从测序数据）的准备：

下载链接：
链接：https://pan.baidu.com/s/1cpLI4RqEs7_Ulkkl76w75g
提取码：eg7g
复制这段内容后打开百度网盘手机App，操作更方便哦

三、创建项目目录

mkdir ~/Seqs/bwa_test

cp GCA_000012525.1_ASM1252v1_genomic.fna bwa_test/

cp test_* bwa_test/

bwa的介绍

1. 即Burrows-Wheeler-Alignment Tool。
2. 是一种能够将差异度较小的序列比对到一个较大的参考基因
   组上的软件包。它由三个不同的算法:
   – BWA-backtrack: 是用来比对 Illumina 的序列的,reads 长度最长
   能到 100bp。
   – BWA-SW: 用于比对 long-read ,支持的长度为 70bp-1Mbp;同时
   支持剪接性比对。
   – BWA-MEM: 和SW都支持较长的read长度,同时都支持剪接性比
   对(split alignments),但是BWA-MEM是更新的算法,也更快,
   更准确,且 BWA-MEM 对于 70bp-100bp 的 Illumina 数据来说,
   效果也更好些。

bwa的安装

sudo apt install bwa

bwa

bwa

bwa的使用

bwa的使用需要两种输入文件:
– 索引的Reference genome data(fasta格式 .fa, .fasta, .fna)
– Short reads data (fastaq格式 .fastaq, .fq)即重测序数据

四、参考基因组索引的建立

bwa index GCA_000012525.1_ASM1252v1_genomic.fna

五、reads比对到参考序列得到sam文件

bwa mem GCA_000012525.1_ASM1252v1_genomic.fna test_7942raw_1.fq.gz test_7942raw_2.fq.gz > test_bwa_7942.sam

2018-11-28 13-32-22屏幕截图.png

SAM文件

1.SAM格式主要应用于测序序列mapping到基因组上的结果表示,当然也可以表示任意的多重比对结果。
2.SAM是一种序列比对格式标准,由Heng Li (Sanger)制定,是以TAB为分割符的文本格式。
3.SAM的全称是sequence alignment/map format。
sam格式详解见https://en.wikipedia.org/wiki/SAM_%28file_format%29

less test_bwa_7942.sam

2018-11-28 13-37-46屏幕截图.png

BAM格式

1.sam是带有比对信息的序列文件(即告诉你这个reads在染色体上的位置等),用于储存序列数据。
2.BAM 也是存储序列比对信息的文件格式,但是是以二进制存储,可以节约空间,计算机可读。
3.二者都是fastq文件经过序列比对或者mapping后输出的格式;其储存的信息都是一致的。

SAMtools

samtools是一个用于操作sam和bam文件的工具合集。

sudo apt install samtools

samtools

2018-11-30 11-19-02屏幕截图.png
Official Document:http://www.htslib.org/doc/samtools.html

bam 与sam的格式转换，bam格式以二进制存储文件，更加节约空间。

为参考基因组建立索引,生成了prefix.fai文件

samtools faidx GCA_000012525.1_ASM1252v1_genomic.fna

sam文件转为bam文件

samtools view -bhS -t GCA_000012525.1_ASM1252v1_genomic.fna.fai -o test_bwa_7942.bam test_bwa_7942.sam

为bam文件排序,sort只能为bam文件排序,而不能为sam;不同版本samtools sort命令的-o参数不同

samtools sort test_bwa_7942.bam -o test_bwa_7942.bam.sorted

为排序的bam文件建立索引. *.bai文件

samtools index test_bwa_7942.bam.sorted

samtools tview:显示reads比对基因组的情况

samtools tview test_bwa_7942.bam.sorted GCA_000012525.1_ASM1252v1_genomic.fna

2018-11-30 11-31-34屏幕截图.png

测试参考基因组每个位点或一段区域的测序深度(不是常说的平均测序深度)

samtools depth test_bwa_7942.bam.sorted >depth.txt

2018-11-30 11-34-07屏幕截图.png

samtools flagstat:统计比对结果

samtools flagstat test_bwa_7942.bam.sorted

2018-11-30 11-35-37屏幕截图.png

-总共的reads数(QC-passed or failed)
-总体上reads的匹配率
-有多少reads是属于paired reads
-reads1中的reads数
-reads2中的reads数
-properly paired:正确配对的reads数量
-with itself and mate mapped:一对reads均比对上的reads数
-singletons:只有单条reads比对上的reads数

samtools rmdup:去除PCR重复

samtools rmdup test_bwa_7942.bam.sorted output.bam

samtools mpileup:生成bcf文件

该命令用于生成bcf文件,再使用bcftools进行SNP和Indel的分析

samtools mpileup -gf GCA_000012525.1_ASM1252v1_genomic.fna test_bwa_7942.bam.sorted > test_bwa_7942.bcf

mpileup生成的结果是pileup格式包含6行:

1. 参考序列名;
2. 位置;
3. 参考碱基;
4. 比对上的reads数;
5. 比对情况;
   – ‘.’代表与参考序列正链匹配。
   – ‘,’代表与参考序列负链匹配。
   – ‘ATCGN’代表在正链上的不匹配。
   – ‘atcgn’代表在负链上的不匹配。
   – ‘*’代表模糊碱基
   – ‘^{’代表匹配的碱基是一个read的开始;’}’后面紧跟的ascii码减去33代表比对质量;这两个符号修饰的是后面的碱基,其后紧跟的碱基(.,ATCGatcgNn)代表该read的第一个碱基。
   – ‘$’代表一个read的结束,该符号修饰的是其前面的碱基。
   – 正则式’+[0-9]+[ACGTNacgtn]+’代表在该位点后插入的碱基;比如上例中在scaffold_1的2847后插入了9个长度的碱基acggtgaag。表明此处极可能是indel。
   –正则式’-[0-9]+[ACGTNacgtn]+’代表在该位点后缺失的碱基;
   6.比对上的碱基的质量

六、基因变异检测软件

从bcf文件(samtools mpileup)中检测基因变异位点: bcftools，从bam文件中检测基因变异位点:GATK，基因变异的主要类型:SNP, indel等。bcftools是samtools中附带的软件,在samtools的安装文件夹中可以找到。

bcftools call -vm test_bwa_7942.bcf -o test_bwa_7942.variants.bcf

bcftools view -v snps,indels test_bwa_7942.variants.bcf>test_bwa_7942.snps.vcf

less test_bwa_7942.snps.vcf

变异位点的过滤——bcftools filter

bcftools filter -o test_bwa_7942.snps.filtered.vcf -i 'QUAL>20 &&DP>5' test_bwa_7942.snps.vcf

这里简单过滤:QUAL小于等于20,DP值小于等于5的位点,-i – include 只保留满足该条件的位点。

七、变异基因注释软件

Annovar安装和运行

wget http://www.openbioinformatics.org/annovar/download/0wgxR2rIVP/annovar.latest.tar.gz

tar -zvxf annovar.latest.tar.gz -C ~/BioSofts/

echo 'export PATH=~/BioSofts/annovar:$PATH'>>~/.bashrc

source ~/.bashrc

Annovar三种注释模式

1. gene-based annotation:判断SNV或CNV是否造成蛋白编码或氨基酸的改变,可用基因命名系统包括RefSeq,UCSC,ENSEMBL,GENCODE, AceView等。
2. region-based annotation:变异位于染色体哪个区域,预测转录因子结合位点、SD区域等
3. filter-based annotation:鉴定在特定数据库中记录的变异,如是否在dbSNP中被报道。
   http://annovar.openbioinformatics.org/en/latest/

Annovar主要程序和目录

-annotate_variation.pl #主程序,功能包括下载数据库,三种不同的注释
-coding_change.pl #可用来推断蛋白质序列
-convert2annovar.pl #将多种格式转为.avinput的程序
-retrieve_seq_from_fasta.pl #用于自行建立其他物种的转录本
-table_annovar.pl #注释程序,可一次性完成三种类型的注释
-variants_reduction.pl #可用来更灵活地定制过滤注释流程
-example/ #存放示例文件
-humandb/ #人类注释数据库

生成annovar输入文件

convert2annovar.pl -format vcf4 test_bwa_7942.snps.filtered.vcf >test_bwa_7942.snps.avinput

Avinput文件,每行代表一个位点,前5列依次为:
– chromosome
– start position
– end position
– the reference nucleotides
– the observed nucleotides
– 其余列:可有可无;如果有,在输出文件中会原样输出。

注释变异基因位点

annotate_variation.pl --geneanno --dbtype refGene --buildver 7942-genome-new test_bwa_7942.snps.avinput ~/BioSofts/annovar/humandb/

--geneanno: 采用gene-based annotation注释模式;
--dbtype :数据库为refGene;还有knowGene,ensGene等
--buildver:为参考基因组版本
最后为数据库所在目录
生成avinput.variant_function和avinput.exonic_variant_function后缀的两个结果文件。

Annovar结果解析

1.avinput.variant_function文件
2.包括所有突变信息:
3.第一列:variant effects, 变异分类,如intergenic, intronic, non-synonymous SNP, frameshift deletion, large-scale duplication等
4.第二列:基因名或Symbol
5.其余列:与avinput输入文件相同:染色体、start、end、ref_nt、obs_nt等

需要下载或自定义注释数据库

1.官方提供了一些基因组注释数据库
http://annovar.openbioinformatics.org/en/latest/user-guide/download/
2.如果没有,这需要自定义注释数据库