生信课程笔记10-变异的识别

宅在家两个多月，不知不觉已经是春天了，也许距离返校的日子更近了吧...

---

变异和突变

变异，指的是实际测序数据与国际规定的参考基因组之间的区别。很多变异其实只是造成人类多样性的原因。突变，指的是那些与疾病相关的变异。
举个例子：ENSEMBL等规定的人类参考基因组文件某位置是AAAAA，然后一个人实际测序得到的序列为AGCAA，那么相比于参考基因组，这个人就有2个变异位点。对于第2个位置，如果查看所有已知的测序，绝大部分人都是G，说明是参考基因组出现了问题，这个变异就不能称作突变。对于第3个位置，如果查看所有已知的测序，绝大部分人都是A，而恰好有一个人不是A，但他是个患者，那么这个变异就是突变了。

变异的类型

SNP（single nucleotide polymorphism）：单核苷酸多态性。个体间基因组DNA序列同一位置单个核苷酸变异（替换、插入或缺失）所引起的多态性。在人类基因组中SNP分布普遍并且密度较大，总数超过107， 平均每300bp（也有说1kbp）就有一个SNP。或称单核苷酸位点变异SNV。
INDEL（insertion-deletion）：插入和缺失。基因组上小片段（>50bp）的插入或缺失。
CNV（copy number variation）：基因组拷贝数变异。基因组中大片段的DNA形成非正常的拷贝数量。比如一个基因在染色体的一条染色单体上的数目为1，但是在染色体复制过程中，复制结束后该基因在染色单体数目由1变成了2或者n。它发生的频率远远高于染色体结构变异，并且整个基因组中覆盖的核苷酸总数大大超过SNP的总数。
SV（structure variation）：结构变异。染色体大片段的插入与缺失，染色体内部的某区域发生翻转颠换，两条染色体之间发生重组。

SNP

一般情况下只分析SNP，其它类型的变异分析有难度或不准确。
来自两个不同个体的DNA片段AAGCCTA和AAGCTTA为等位基因。几乎所有常见的SNP位点只有两个等位基因。
在人体中，SNP的发生机率大约是0.1%，也就是每1000个碱基对就可能有一个SNP（密度高）。对疾病发生和药物治疗有重大影响的SNP，估计只占数以百万计SNP的很小一部分。
SNP位点的分布是不均匀的，在非转录序列比在转录序列更常见。编码区的单核苷酸多态性——编码 SNP（coding SNP，cSNP）也有同义和非同义两种类型，非同义SNP会改变蛋白质的氨基酸序列。基因非编码区、基因间隔区的SNP仍然可能影响转录因子结合、剪接等过程。
从演化的观点来看，SNP具有相当程度的稳定性，即使经过代代相传，SNP所引起的改变却不大，因此可用以研究族群演化。

基于测序数据检测SNP的方法

基于测序数据检测结构变异SV的方法

---

以HISAT2+Samtools+bcftools为例，识别变异的流程

HISAT2

HISAT2 是一款利用改进的BWT算法进行序列比对的软件。由约翰霍普金斯大学计算生物学中心（CCB at JHU）开发，是TopHat的升级版本，速度提高了50倍。利用 HISAT2 + StringTie 流程，可以快速地分析转录组测序数据，获得每个基因和转录本的表达量。

首先需要构建参考基因组索引用于下一步的比对。HISAT2提供了两个脚本用于从基因组注释GTF文件中提取剪接位点和外显子位置，基于这些特征，可以使 RNA-Seq reads 比对更加准确。然后再进行reads mapping。

#（1）下载参考基因组和注释文件
# 例如：dmel.genome.fa.gz和dmel.annotation.gtf.gz
wget  dmel.genome.fa.gz
gzip  -d  dmel.genome.fa.gz  #或gunzip  dmel.genome.fa.gz
#（2）提取外显子和剪接位点
extract_splice_sites.py  dmel.annotation.gtf  >  dmel.ss
extract_exons.py  dmel.annotation.gtf  >  dmel.exon
extract_snps.py  snp142Common.txt  >  genome.snp
#（3）构建参考基因组索引
hisat2-build  --ss  dmel.ss  --exon  dmel.exon  (--snp  genome.snp)  dmel.genome.fa  dmel
# 得到8个索引文件，以dmel为共同的前缀，dmel.1~8.ht2
#（4）也可以下载现成的index
wget  bdgp6.tar.gz
tar  -zxvf  bdgp6.tar.gz
#（5）reads mapping
hisat2  -p  8  (--dta)  -x  dmel  -1  sample_1.fq.gz  -2  sample_2.fq.gz  -S  sample.sam  (&)
# hisat2  [options]*  -x   {-1   -2  | -U  | --sra-acc }  [-S ]
# -p  线程数
# -t  记录时间
# --dta   输出专门为转录本组装的比对结果（如果比对结果后续需要使用StringTie进行转录本组装，则需要加入--dta选项）
# -x  指定基因组索引的文件前缀
# -S  指定输出的SAM文件。默认输出到标准输出，比对结束后统计结果输出到标准错误输出。
# -f  表示输入文件格式为fasta（默认），-q表示输入文件格式为fastq。可以是gzip压缩的文件。
# -phred33  输入的FASTQ文件碱基质量值编码标准为phred33（默认）。

比对结果：

Samtools

SAM（sequence Alignment/mapping）数据格式是目前高通量测序中存放比对数据的标准格式。BAM是SAM的二进制格式。使用samtools将sam文件转化为bam文件，并进行排序。

#（1）将sam文件转化为bam文件
samtools view  -bS  sample.sam  >  sample_unsorted.bam
# -b：输出为bam格式
# -S：出入为sam格式
#（2）将bam文件进行排序
samtools sort  -@  8  sample_unsorted.bam  -o  sample.sorted.bam
# -o：输出文件的名字
# -@：线程数
#默认按照染色体位置进行排序，-n是根据read名进行排序，-t 根据TAG进行排序。
#（3）版本1.4后可直接简化
samtools sort  -@  8  -o  sample.sorted.bam  sample.sam
#（4）索引
samtools index  sample.sorted.bam
#得到sample.sorted.bam.bai索引文件

SAM文件：

bcftools

vcf格式（Variant Call Format）是存储变异位点的标准格式，用于记录variants（SNP / InDel）。BCF是VCF的二进制文件。

#（1）SNP calling
# mpileup命令：得到染色体上每个碱基的比对情况的汇总
bcftools mpileup  -Ob  -o  sample.bcf  -f  dmel.genome.fa  sample.sorted.bam
# 输入BAM文件sorted.bam
# -f  --fasta-ref  指定参考序列的fasta文件
# -O  指定输出文件的类型，压缩的BCF(b)，未压缩的BCF(u)，压缩的VCF(z)，未压缩的VCF(v) 
# -o  指定输出文件的名字sample.bcf
# call命令：执行SNP calling 
bcftools call  -vmO  z  -o  sample.vcf.gz  sample.bcf
# -v  只输出变异位点的信息，如果一个位点不是snp/indel则不会输出
# 两种calling算法：-c参数对应consensus-caller算法，-m参数对应multiallelic-caller算法，后者更适合多种allel和罕见变异的calling
#（2）tabix
tabix  -p  vcf  sample.vcf.gz
# 输入为压缩文件vcf.gz，生成的索引文件为sample.vcf.gz.tbi
#（3）index对vcf文件建立索引
bgzip  view.vcf  #输入的VCF文件必须是bgzip压缩后的文件
gunzip  view.vcf.gz  #解压缩
bcftools index  view.vcf.gz  #生成索引文件view.vcf.gz.csi
bcftools index  -t  view.vcf.gz  #生成索引文件view.vcf.gz.tbi
#（4）query通过表达式来指定输出格式
bcftools query  -f  '%CHROM/t%POS/t%REF/t%ALT[/t%SAMPLE=%GT]/n'  view.vcf.gz
# -f  通过一个表达式来指定输出格式
# %CHROM  代表VCF文件中染色体那一列，其他的列POS, ID, REF, ALT, QUAL, FILTER也是类似的写法
# []  对于FORMAT字段的信息用中括号括起来
# %SAMPLE  代表样本名称
# %GT  代表FORMAT字段中genotype的信息
# /t  代表制表符分隔
# /n  代表新的一行
#（5）sort按照染色体位置进行排序
bcftools sort  view.vcf.gz  -o  sort.view.vcf
#（6）filter过滤不可靠位点
bcftools filter  -O  z  -o  sample_filtered.vcf.gz  -s  LOWQUAL -i'%QUAL>10'  sample.vcf.gz
# -O  --output-type 输出的格式，z和v都行，压缩的VCF(z)，未压缩的VCF(v)
# -o  --output 输出文件的名称
# -s  --soft-filter 将过滤掉的位点用字符串注释
#（7）view命令用于VCF和BCF格式的转换
bcftools view  view.vcf.gz  -O  u  -o  view.bcf  
bcftools view  view.vcf.gz  -s  NA00001,NA00002  -o  subset.vcf 
bcftools view  view.vcf.gz  -k  -o  known.vcf 
# -O  指定输出文件的类型，压缩的BCF(b)，未压缩的BCF(u)，压缩的VCF(z)，未压缩的VCF(v) 
# -o  指定输出文件的名字
# -s  想要保留的样本信息，多个样本用逗号分隔；如果样本名称添加了^前缀，代表去除这些样本，比如-s ^NA00001,NA00002
# -k  表示筛选已知的突变位点，即ID那一列值不是.的突变位点
#（8）stats命令用于统计VCF文件的基本信息
bcftools stats  view.vcf  >  view.stats
bcftools stats  -F  dmel.genome.fa  -s  -  sample.vcf.gz  >  sample.vcf.gz.stats 
# -F  --fasta-ref  faidx indexed reference sequence（参考基因组） file to determine INDEL context
# -s  list of samples for sample stats, "-" to include all samples 统计的样本列表，-代表所有样本
#（9）plot-vcfstats命令进行可视化
plot-vcfstats  sample.vcf.gz.stats  -p  plots/sample.vcf.gz.stats
# -p  指定输出结果的目录
#（这个脚本位于bcftools安装目录的misc目录下，依赖latex 生成pdf 文件。这里需要matplotlib）

stats统计文件：