在实验流程中,由于气溶胶、或者实验操作不当等,可能会造成样本间存在污染,在设计的Panel中,加入一批人群频率为0.5左右的SNP位点,这种位点的检出丰度在理想状态下是0,50%,100%。如果检出的值偏差较大,则提示可能存在污染。
因此,在设计时,应考虑的是不能选择高GC区域,同时要避开容易发生CNV(或LOH)的区域。在大Panel中,应考虑挑选位点的均匀分布。
首先是考虑中国人群频率的获得,1000G的CHB比较旧,近期的研究有女娲基因组,共纳入2999个中国人的全基因组数据。目前只提供hg38数据下载(并且还下不了)。
另外还有WBBC(施一公当校长的那个西湖大学),共纳入14726个样本。目前提供了hg19和hg38两个版本的数据下载。因此这里采用WBBC的hg19数据。
只需要常染色体的0.49~0.51人群频率的位点,使用python过滤获得
import gzip
output = open("WBBC.AF.vcf", "w")
for i in range(1, 23):
with gzip.open("WBBC.chr{}.GRCh37.vcf.gz".format(str(i)), "rt") as f:
for line in f:
if line.startswith("#"):
if i == 1:
output.write(line)
else:
lines = line.split("\t")
ref = lines[3]
alt = lines[4]
info = lines[7]
if len(ref) != 1 or len(alt) != 1:
continue
AF = 0
for k in info.split(";"):
if k.startswith("AF="):
AF = float(k.replace("AF=", ""))
if AF >= 0.49 and AF <= 0.51:
output.write(line)
output.close()
然后过滤掉GC比例差的区域。使用bedtools对参考基因组获得GC比例在0.40~0.46的bed,获取此范围下的SNP位点。
bedtools makewindows -g grch37.fasta.fai -w 10000 > grch37.bin.bed
bedtools nuc -fi grch37.fasta -bed grch37.bin.bed |\
awk '{if(($5>=0.40) && ($5<=0.46)) print $0}' > grch37.normalGC.bed
bedtools intersect -a WBBC.AF.vcf -b grch37.normalGC.bed -header > WBBC.AF.GC.vcf
再然后过滤掉CNV易发区域,由于不知道有没有现成的数据库,这里直接去掉Decipher数据库中的相关基因区域。需要注意的是Decipher的坐标是hg38的,需要自行转换。
bedtools intersect =a WBBC.AF.GC.vcf -b decipher.bed -v -header > WBBC.AF.GC.decipher.vcf
获得位点列表,最后使用随机挑选若干位点即可。