生信课程笔记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 <ht2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <sam>]
# -p 线程数
# -t 记录时间
# --dta 输出专门为转录本组装的比对结果（如果比对结果后续需要使用StringTie进行转录本组装，则需要加入--dta选项）
# -x 指定基因组索引的文件前缀
# -S 指定输出的SAM文件。默认输出到标准输出，比对结束后统计结果输出到标准错误输出。
# -f 表示输入文件格式为fasta（默认），-q表示输入文件格式为fastq。可以是gzip压缩的文件。
# -phred33 输入的FASTQ文件碱基质量值编码标准为phred33（默认）。
比对结果：
比对输出内容.PNG
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文件：
SAM文件.PNG
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统计文件：
stats结果.PNG
