昨天游侠发布了一个求助帖《
英雄帖
|
寻找生信高手解决一个位点过滤的问题
》后，有很多朋友都积极响应，很快就有朋友提供了代码，在征得原创作者的同意后，发布前三位生信高手的脚本。
第一位
率先完成的是来自
迪安诊断
的
梁泽滨
具体用法为：
python trans_hom_call.py in.vcf out.vcf threshold
其中
threshold
为
allele
碱基占总
depth
比率，运行脚本前需要首先安装
pyvcf
包。
代码如下：
import sys
from vcf import VCFReader, VCFWriter
from vcf.model import make_calldata_tuple
def make_hom_calldata(format, data):
CallData =make_calldata_tuple(format.split(':'))
call_data ={k: getattr(data, k) for k in format.split(':')}
call_data.update(GT='0/0')
returnCallData(**call_data)
if __name__ == '__main__':
vcf_in,vcf_out, threshold = sys.argv[1:]
threshold =float(threshold)
vcf_reader =VCFReader(filename=vcf_in)
withopen(vcf_out, 'w') as fp:
vcf_writer = VCFWriter(fp, vcf_reader)
forrecord in vcf_reader:
forsample in record.samples:
if sample.data.AD and sample.data.DP and sample.data.AD[1] / sample.data.DP< threshold:
sample.data = make_hom_calldata(record.FORMAT, sample.data)
vcf_writer.write_record(record)
第二位
完成的是来自
智因东方
的
杨双浩
具体用法为：
perl yingxiongbang.pl yingxiongbang.vcf>log &
说明：
perl yingxiongbang.pl
usage:
perl yingxiongbang.pl *.vcf
exception:
chr1 15438825 . ATGTGTGTGTGTGTG ATGTGTGTG,ATGTG,ATGTGTG,A : sum(ATGTG,ATGTGTG,A)/AD<0.3=> 0/0
log
输出说明：
行信息：
chr1 12368705 . GT GTT,G 569.06 . AC=9,4;AF=0.300,0.133;AN=30;BaseQRankSum=0.087;DP=606;ExcessHet=24.4476;FS=0.000;InbreedingCoeff=-0.7074;MLEAC=9,4;MLEAF=0.30
########(2 2)4/32;
0/0:9,2,2:32:4:4,0,187,11,148,216
0/1:9,2,2:32:4:4,0,187,11,148,216 (old)
括号中为去掉
ref
深度之外的基因型深度数组，
4/32 :
表示括号中深度之和
/AD
如果
4/32 <0.3
则处理为
0/0
。
具体代码如下：
if (@ARGV==0)
{
print"usage:\n perl $0*.vcf\nexception:\n chr1 15438825 . ATGTGTGTGTGTGTG ATGTGTGTG,ATGTG,ATGTGTG,A :sum(ATGTG,ATGTGTG,A)/AD<0.3 => 0/0\n";
exit;
}
$vcf="$ARGV[0]";
open(I,$vcf);
open (O,">$vcf.out");
while(
)
{
chomp;
print"$_\n";
if(/^\#/)
{
print O"$_\n";
}
else
{
@a=split;
#print"@a\t";
foreach(@a)
{
if(/\//&& /:/)
{
#GT:AD:DP:GQ:PL 1/4:0,3,0,0,4:11:99
@b=split/\:/,$_;
$type=shift@b;
$refaltd=shift @b;
$pl=pop@b;
$gq=pop@b;
$dp=pop@b;
my$sumd=0;
#my$altd="";
my@c=split/\,/,$refaltd;
$refd=shift@c;
foreach$s (@c) {
$sumd+=$s;
#$altd.="$s,";
}
#$altd=~s/\,$//;
if($dp==0 || !defined $dp) {
printO "$_\t"; #update
}
else
{
$r=$sumd/$dp;
if($r<0.3)
{
print"########(@c)$sumd/$dp;\n0/0:$refaltd:$dp:$gq:$pl\n$_ (old)\n\n";
printO "0/0:$refaltd:$dp:$gq:$pl\t";
#die;
}
else
{
printO "$_\t";
}
}
}
else
{
printO "$_\t"
}
}
print O"\n";
}
}
第三位
完成的是来自
嘉宝仁和
的
仝微微
用法为：
perl yingxiongbang.pl in.VCF > out.VCF
脚本为：
#!usr/bin/perl -w
use strict;
die "usage:$0
\n" if @ARGV<1;
open AA,'
while(
){
chomp;
if(/^#/){
print"$_\n";
next;
}else{
#print"$_\n";
