美文网首页
VCF文件添加ID

VCF文件添加ID

作者: Zukunft_Lab | 来源:发表于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