挑选一批SNP位点用于分析污染

在实验流程中,由于气溶胶、或者实验操作不当等,可能会造成样本间存在污染,在设计的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

获得位点列表,最后使用随机挑选若干位点即可。