问题: 修改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
结果
结果示例,来源网络
网友评论