在2019年的一篇NG文章"Resequencing of 414 cultivated and wild watermelon accessions identifies selection for fruit quality traits"的方法部分，作者提到他們使用位于四倍兼并位點(diǎn)( fourfold degenerated sites)上的SNP構(gòu)建系統(tǒng)發(fā)育樹。原本的19,725,853個(gè)SNP經(jīng)過這一過濾，最終只剩下89,914 個(gè)SNP。

Phylogenetic and population analyses. Bi-allelic SNPs with a missing data rate less than 15% and a minor allele count greater than three were kept for population genomic analyses. Additionally, only SNPs at fourfold degenerated sites (89,914 SNPs) were used to construct a neighbor-joining phylogenetic tree using MEGA7 with 500 bootstraps61.

所謂四倍兼并位點(diǎn)(4DTV), 指無論怎么突變，最終氨基酸都一樣的位點(diǎn)，例如"TCA/TCT/TCC/TCG"中的第三位就是一個(gè)4DTV，只要在這個(gè)位點(diǎn)的堿基，無論怎么突變最終都是絲氨酸(Ser)。

我稍微檢索了一下，發(fā)現(xiàn)并沒有相應(yīng)的腳本，于是就有了這道練習(xí)題。

給定一個(gè)用SnpEff注釋的VCF文件，請(qǐng)過濾出其中的位于四倍兼并位點(diǎn)上的SNP。

我花了一個(gè)晚上時(shí)間寫了代碼，使用方法為python calc_4dTv_in_eff_vcf.py input.vcf output.vcf ref.fa，三個(gè)參數(shù)分別是輸入的vcf，輸出的vcf文件名，以及參考基因組序列。

需要注意的是，通常SnpEff會(huì)對(duì)一個(gè)位點(diǎn)注釋多個(gè)信息，而我只選擇其中第一條注釋信息進(jìn)行過濾。

#!/usr/bin/env python3from sys import argv from pysam import VariantFile from pysam import FastaFilefile_in = argv[1] file_out = argv[2] fafile = argv[3]codon = set(["TC", "CT", "CC", "CG", "AC", "GT", "GC", "GG"]) rev_dict = dict(A='T',T='A', C='G', G='C')bcf_in = VariantFile(file_in) bcf_out = VariantFile(file_out, "w", header = bcf_in.header) fa_in = FastaFile(fafile)for rec in bcf_in.fetch():ann = rec.info['ANN']info = rec.info['ANN'][0].split('|')# only use synonymouse variantsif info[1] != "synonymous_variant":continue# only the 3rd position can be 4dTvif int(info[9][2:-3]) % 3 != 0:continue# determine the strand by the REF column and mutation# if the ref is not same as the mutation siteif rec.ref == info[9][-3]:pre = fa_in.fetch(rec.chrom, rec.pos-3, rec.pos-1)else:tmp = fa_in.fetch(rec.chrom, rec.pos, rec.pos+2)tmp.upper()pre = rev_dict[tmp[1]] + rev_dict[tmp[0]]if pre not in codon:continuebcf_out.write(rec)

總結(jié)

以上是生活随笔為你收集整理的「生信练习题」从SnpEff注释得到的VCF中过滤4DTV位点的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。