「群体遗传学实战」第一课: 对SNP位点进行注释

文章

我们用于实战的数据集来自于2019年发表于NG的西瓜文章，它提供了GATK过滤后的SNP数据集用于分析，并且附录提供了完整的样本信息。该SNP数据集包括414个西瓜，2000万个SNP信息，数据大小为22G.

数据准备

根据文章提供的下载地址，我们分别下载西瓜的基因组，GFF注释文章和存放的VCF的数据集。

ftp://cucurbitgenomics.org/pub/cucurbit/genome/watermelon/97103/v2/97103_genome_v2.fa.gz
ftp://cucurbitgenomics.org/pub/cucurbit/genome/watermelon/97103/v2/97103_gene_gff_v2.gz
ftp://cucurbitgenomics.org/pub/cucurbit/reseq/watermelon/v2/watermelon_414acc_SNP.vcf.gz

需要注意的是，他们提供的SNP文件存在一些问题，不能用BCFtools进行解析，不过很容易解决，只需要运行如下命令

zgrep -v '##INFO' watermelon_414acc_SNP.vcf.gz | bgzip -c > watermelon_414acc_SNP2.vcf.gz &

使用SnpEff注释VCF文件

在文章的第二部分Genome variation map and phylogeny of Citrullus species.，作者对414个重测序结果分析得到SNP信息做了如下描述。

These accessions were sequenced to an average depth of 14.5× and coverage of 92.2% of the ‘97103’ genome. In total, 19,725,853 SNPs were identified, of which 1,100,803 were located in coding regions, causing 502,028 nonsynonymous mutations, 589,735 synonymous mutations, 1,031 start codon changes and 6,808 stop codon changes. Furthermore, 6,675,290 small indels were identified, of which 56,115 were located in coding regions.

很可惜，作者并没有在文章中指出他是使用什么软件对SNP进行注释，可能是自己写了一个脚本进行分析，而我为了偷懒，直接用一个现成的软件，snpEFF，对VCF文件中的SNP进行注释。

参考使用snpEff对VCF进行注释，我们需要先建立西瓜的注释数据库。

第一步，修改snpEff/snpEff.config，增加西瓜基因组信息

# watermelon
watermelon.genome : watermelon

第二步，在snpEff/data/新建一个watermelon文件夹，并存放基因组序列和GFF文件

mkdir -p snpEff/data/watermelon
cp  97103_genome_v2.fa.gz snpEff/data/watermelon/sequences.fa.gz
cp 97103_gene_gff_v2.gz  snpEff/data/watermelon/genes.gff.gz

第三步，运行构建命令

java -jar snpEff/snpEff.jar build -gff3 -v watermelon
``

最后，我们就可以对VCF进行注释了

```bash
java -jar /opt/biosoft/snpEff/snpEff.jar ann -no-utr -no-downstream -no-upstream -no-intergenic watermelon watermelon_414acc_SNP2.vcf.gz > watermelon_414acc_eff.vcf &

这一步会很久，最终会得到一个204G文件，这就是我们的注释结果。而另外的snpEff_genes.txt和snpEff_summary.html则是用来了解SNP的分布情况。

SNP分布图

根据SnpEff的注释结果（如下），落在编码区有1,168,139个。

Type (alphabetical order)		Count	Percent
DOWNSTREAM		7,325,679	20.41%
EXON		1,168,139	3.254%
INTERGENIC		16,068,984	44.769%
INTRON		3,719,155	10.362%
SPLICE_SITE_ACCEPTOR		1,317	0.004%
SPLICE_SITE_DONOR		1,415	0.004%
SPLICE_SITE_REGION		78,665	0.219%
UPSTREAM		7,530,074	20.979%
UTR_5_PRIME		3	0%

对于落在编码区域的SNP而言，大约有533,656个位点是非同义突变，638,879个位点是同义突变。

Type (alphabetical order)		Count	Percent
MISSENSE		533,656	45.277%
NONSENSE		6,120	0.519%
SILENT		638,879	54.204%

相比较于文章，我们的各个位置上的SNP信息都比较多，但是整体在一个数量级上，接下来就是利用这些SNP位点做更加深入的分析

• 系统发育分析
• 主成分分析
• 群体结构分析
• 基因流
• ...

这些都会在将来慢慢更新。