Thank you Dear For this guidance
Dear All I have a combined VCF file of few individuals. In this VCF instead of having normal CHR1, chr2 notions for chromosomes it is having the chromosome information as
gi|996703411|ref|NW_015379183.1|, gi|996703411|ref|NW_015379175.1
In this notion NW_015379183.1 corresponds to a specific Chromosome. The same is true for its positions, If I have the chromosome numbers for all gi|996703411|ref|NW_015379183.1| sort of notions how I can replace the chromosome names to the original names.
3 answers
use bcftools annotate https://samtools.github.io/bcftools/bcftools.html
with
--rename-chrs file rename chromosomes according to the map in file, with "old_name new_name\n" pairs separated by whitespaces, each on a separate line.
Please try not to use "dear", it sounds inappropriate in a professional scientific forum. Thank you!
As said above bcftools will be the best way to do it. If you want a code version here is mine
Assuming that you have a match_table.txt file like this (separate by tab) :
gi|996703411|ref|NW_015379183.1| chr1
gi|996703411|ref|NW_015379184.1| chr2
gi|996703411|ref|NW_015379185.1| chr3
gi|996703411|ref|NW_015379186.1| chr4
gi|996703411|ref|NW_015379187.1| chr5
Coding version in python :
###Create a dictionnary containing your match_table.txt
match_dict={}
###Open your match table
with open("matching_table.txt", 'r') as match_f:
###For each line, you create a key/value item in a dictionnary
for line in match_f:
gi_notation = line.rstrip().split("\t")[0]
chr_notation = line.rstrip().split("\t")[1]
###Check if the key doesn't exist in the dictionnary
if gi_notation not in match_dict:
match_dict[gi_notation] = chr_notation
else:
print("Care, duplicate in matching_table.txt, on : "+str(gi_notation))
###Open your vcf file
new_vcf_file = open("your_new_vcf_file.vcf", "a")
with open("your_vcf_file.vcf", 'r') as vcf_f:
###Read it line by line
headers_chromosome = ""
for line in vcf_f:
###Change VCF dictionnary headers
if line.startswith('##contig'):
###Get chromosome name
headers_chromosome = line.split("=")[2].split(",")[0]
###If your chromosome exist in your dictionnary
if headers_chromosome in match_dict:
###Replace in chromosome name in line
line = line.replace(headers_chromosome, match_dict[headers_chromosome])
###Skip metadata informations
if line[0] != '#':
###Retrieve your chromosome for each line
current_chromosome = line.split("\t")[0]
###If your chromosome exist in your dictionnary
if current_chromosome in match_dict:
###Change the value of your chromosome
new_vcf_file.write(match_dict[current_chromosome]+"\t"+'\t'.join(line.split("\t")[1:]))
###Your chromosome is not in your dictionnary (I write it as it is but you can do something else...)
else:
print("This chromosome is not in my matching_table.txt : "+str(current_chromosome))
new_vcf_file.write(line)
###Write unchanged metadata
else:
new_vcf_file.write(line)
new_vcf_file.close()
this code doesn't change the VCF dictionary headers ("## contig=....")
Right, i'll add it ! Thanks
I don’t think we should be providing code - how would OP learn if the solution is ready-made for them?
The way I think is a little bit different. I personally needed to rely on some code to understand how it works and how to do my own solution. In that way, I let OP some comments on key lines. Moreover, if OP just want a solution without getting involved in bioinformatics concept, here is one. Depends on how and what OP wants to learn, if OP want to.
Fair enough, I'm opposed to that concept a little as I have seen way too many people take advantage of people trying to help such as yourself and outright make demands for code or even subtly wait until code is provided.
Thanks Bastien for this help
Dear I am trying with the script you have generously created for me, But I am receiving the following error IndentationError: unexpected indent
Chr1 = line.rstrip().split("\t")[1]although I am trying to maintain the correct indentation in command
for line in match_f:
gi_notation = line.rstrip().split("\t")[0]
chr_notation = line.rstrip().split("\t")[1]
I found it, it works well for me now. Note that, it is easier to find an error by providing the error message and the line involved. As Ram said, if you want to go deeper in file manipulation, you have to try coding by your own. Thereby, I let you some comments in the code '###' to have a better understanding on the process. This code will only work for what you ask but you can reuse some line to write a new script, like opening a file (with open("your_vcf_file.vcf", 'r') as vcf_f) or loop on the different lines (for line in match_f)
Create a file with 2 columns - the gi| notation and the notation you want.
Read the above file to a dictionary. Then, read the VCF file in a streaming fashion and substitute each occurrence of the old notation with the new one.
This is the strategy. I'd recommend you coding this yourself so you can use it as a learning experience.
Thank you Sir for this advice I will try to do it
Log in to answer this question.
see VCF files: Change Chromosome Notation