1、输入文件：Tophat2等软件产生的比对结果文件（BAM格式），和参考基因组序列文件genome.fa
2、首先要对genome.fa文件建立索引，如：samtools faidx genome.fasta生成的索引文件以.fai后缀结尾。
3、对BAM文件进行排序，如：samtools sort abc.bam abc.sort；具体见3_BAMsort.sh
4、使用samtools的mpileup模块生成一个bcf文件，然后再使用bcftools命令对bcf文件进行处理，此处主要
使用view命令来进行SNP和Indelcalling，具体见4_call_snp.sh
1）mpileup 参数如下：
-g ：计算基因型的似然度，输出于BCF文件中
-u ：指定输出文件为BCF格式
-S ：产生每个样本中链偏见的P-value
-D ：产出每个样本的read的深度
-r ：指定特定的区域进行SNP，INDEL
-f ：指定建好索引的参考序列文件
-q ：指定可用比对的最低映射质量，默认为0
-Q ：指定可用reads的最低碱基质量，默认为13
其它参数参考官网指导手册
2）bcftools 参数如下：
view:这里主要使用它来处理bcf文件
-c ：指定用贝叶斯推理来鉴定变异，并且会自动调用-e参数（指定使用最大似然法）
-v ：仅仅输出变异的位点信息
-N ：忽略参考基因上不是A/T/C/G 的位点信息
-g：计算突变位点的每个样本的基因型
-b ：指定产出的文件格式是BCF的，默认为VCF格式。
其它参数参考官网指导手册
3）示范程序如下：4_call_snp.sh
#!/bin/bash -x
cd/home/cuckoo/data/liyan/RNAseq_train/SNPs
samtoolsmpileup -guSDf genome.fa s0924fb_sorted.bam sCal27_sorted.bam |/home/cuckoo/software/RNAseq/samtools-0.1.19/bcftools/bcftoolsview -cvNg - >sample_raw.vcf
5、过滤第四步产生的VCF文件：减少SNPs的假阳性率
1）过滤掉突变质量低于30的，reads深度小于5的SNPS位点：
输入文件：上一步的sample_raw.vcf文件
输出文件：过滤后的sample_filter.vcf文件
2）过滤掉reads深度大于100的SNPs位点：
/home/cuckoo/software/RNAseq/samtools-0.1.19/bcftools/bcftools viewsample_raw.vcf |
/home/cuckoo/software/RNAseq/samtools-0.1.19/bcftools/vcfutils.plvarFilter -D100 > sample_filter.vcf
3）第8列的DP4进行过滤：
vcf文件中的 DP4一列信息非常重要，提供了4个数据：1：比对结果和正链一致的reads数、
2：比对结果和负链一致的reads数、3：比对结果在正链的variant上的reads数、4：比对结果在负链的variant上的reads数。
可以设定 （value3+ value4）大于某一阈值，才算是variant。比如：
perl -ne'print $_ if /DP4=(d+),(d+),(d+),(d+)/ &&($3+$4)>=10 && ($3+$4)/($1+$2+$3+$4)>=0.8'smaple_raw.vcf > sample_filter.vcf
注：以上的过滤方法根据不同的实验设计和需求而异，但最终都会生成一个过滤后的文件：snp_filter.vcf，用于下一步的统计分析。
6、（可忽略这一步）对过滤后的VCF文件进行解析，统计突变的信息：染色体，位点，参考碱基，变异碱基，突变质量
具体见6_parse_vcf.pl文件：
输入文件：sample_fiter.vcf
输出文件：snp.txt
6、注释过滤后的VCF文件：工具-ANNOVAR
1、数据库准备：下载相应的注释文件到指定目录如：/data2/disk1/humandbAnnovar/，一般为annovar安装目录下的humandb目录。下载数据库文件地址：http://www.openbioinformatics.org/annovar/download/hg19_ALL.sites.2010_11.txt.gz
或是通过命令：./annotate_variation.pl -downdb -buildverhg19 -webfrom annovar refGene humandb/
2、转换数据格式：第五步过滤后vcf文件（sample_filter.vcf）换为ANNOVAR指定的格式
命令：convert2annovar.pl -formatvcf4 sample_filter.vcf >smaple_annovar.input
3、注释：使用table_annovar.pl进行注释，因为它可以一次性完成位点、基因、区域的注释工作。
1）table_annovar.pl使用方法：/home/cuckoo/software/annovar/table_annovar.plsample_annover.input /data2/disk1/humandbAnnovar/ -buildver hg19-out sample -remove -protocolrefGene,cytoBand,1000g2014oct_all,snp138 -operation g,r,f,f-nastring NA
2）table_annovar.pl参数说明：
输入文件：sample_annover.input
数据库位置：/data2/disk1/humandbAnnovar/，一般为annovar安装目录下的humandb目录。
-buildver：指定物种基因组的类型，hg19
-out：指定输出文件的前缀，sample；默认为txt格式。
-remove：指定删除所有的临时文件。
-protocol：指定参考文件的数据库来源，用逗号隔开;refGene,cytoBand,1000g2014oct_all,snp138四个必须包括。
-operation：指定与-protocol所对应的参考文件的类型，用逗号隔开。
-nastring：指定结果文件中的缺失值表示为NA
4、其它注释方法：annotate_variation.pl有三种注释模式:gene-based（默认）,region-based,filter-based
使用方法：
1）gene-based：annotate_variation.pl -buildver hg19 ex1.avinputhumandb/
2）region-based：
annotate_variation.pl -regionanno -buildver hg19 -dbtype cytoBandex1.avinput humandb/
annotate_variation.pl -regionanno -buildver hg19 -dbtype gff3-gff3dbfile tfbs.gff3 ex1.avinput humandb/
3）filter-based：
annotate_variation.pl -filter -dbtype 1000g2014oct_all -maf 0.01ex1.avinput humandb/
annotate_variation.pl -filter -buildver hg19 -dbtype snp138ex1.avinput humandb/
annotate_variation.pl -filter -dbtype ljb26_all -otherinfoex1.avinput humandb/
5、注意事项：
1）指导手册：table_annovar.pl和annotate_variation.pl的指导手册和参数说明可以通过输入相应的文件名获得。
2）每个文件的具体参数及使用方法请参考官网教程及指导手册。
7、统计第六步中结果文件的突变类型:转换和颠换，并且用R作条形图展示。
具体见7_snp_count.pl文件:
输入文件：snp.txt
输出文件：snp_type.txt
8、R语言画条形图展示突变类型，脚本8_snp_count.r
