美文网首页
VCF文件添加ID

VCF文件添加ID

作者: 糖异生的鱼 | 来源:发表于2022-10-04 11:14 被阅读0次

问题: 修改VCF文件中ID列以“.”标识为“染色体号-SNP的位置”

示例图片,来源网络

1.0 python

脚本
#author_zukunft
#data_202210

import vcf
import argparse

parser = argparse.ArgumentParser(description='add ID to VCF_file: ".">>"chr*-POS"')
parser.add_argument('--input', type=str, nargs='?',
                    help='input file')

parser.add_argument('--output', type=str, nargs='?',
                    help='output file')

args = parser.parse_args()

input_file = args.input
output_file = args.output


vcf_reader = vcf.Reader(filename=input_file)
vcf_writer = vcf.Writer(open("tsp.2-8.rb17.dp5_indv2_maf5.impute.clear_site.snp.nochr.ID.vcf", 'w'), vcf_reader)

if __name__ == "__main__":
    line=0
    for record in vcf_reader:
        # print(record.ID)
        record.ID = str(record.CHROM) + "_"+ str(record.POS) #自由修改
        # print(record.ID)
        line+=1
        # print("正在执行第" + line + "SNP")
        vcf_writer.write_record(record)

vcf_writer.close()
print("done,The numbers of SNPs:")
print(line)
运行
tools_py=/path/vcf_add_ID.py #vcf_add_ID.py见下面

input_vcf=/path/**.vcf.gz  #**.vcf也可
output_vcf=/path/**.ID.vcf.gz

time python3 $tools_py --input $input_vcf --output $output_vcf

2.0 awk

#严格来讲,这是一个投机取巧的操作,用来查看vcf头文件的行数以便后面操作;**.vcf是最原始文件 输出header 
input_vcf=/path/**.vcf.gz
output_vcf=/path/**.awk.body.ID.vcf

touch/path/snp.list.txt #空文件
vcftools --gzvcf $input_vcf --snps /path/snp.list.txt --recode --recode-INFO-all --out /path/header
wc -l path/header

#如果头文件上面输出的合适 可以不多这一步 n:头文件的行数
zless /path/**.vcf.gz | head -n *** > /path/header

#提取除header以外的内容 进行ID操作 输出到不包括header的临时文件(body)
less $input_vcf | grep "^#" -v | awk '{$3=$1"-"$2;print $0}' | sed 's/ /\t/g' > $output_vcf

#合并header + body
finally_output_vcf=/path/**.awk.ID.vcf
cat /path/header $output_vcf > $finally_output_vcf

结果

结果示例,来源网络

相关文章

网友评论

      本文标题:VCF文件添加ID

      本文链接:https://www.haomeiwen.com/subject/lutpartx.html