Hi,
I'm creating a fasta alignment of concatenated SNPs from a combined vcf file, but I'm having some trouble with the gaps.
I'm using the useful script vcf_tab_to_fasta_alignment.pl (Christina Bergey, 2012, url = "http://code.google.com/p/vcf-tab-to-fasta" ). From an input.vcf.gz it I obtain all SNPs aligned to the reference in fasta format (see below). However, the SNPs coordinates are different between isolates so the loci with no SNPs info are filled with gaps (-).
I get:
>REF ATCCTTGCA
>ID1 -CT-A-CT-
>ID2 CG---CAA-
Is it possible to fill in those gaps using the nucleotide in the reference genome? I imagine that it can be achieved by using a simple script but I'm not a programmer.
I wish:
>REF ATCCTTGCA
>ID1 ACTCATCTA
>ID2 CGCCTCAAA
Any help would be very much appreciated! :D
Many thanks,
A
Since you are already having fasta file, I would recommend you to simply have your reference + your fasta sequences, do a codon aware alignment or any alignment of your choice.
After that, use a manual editing in Aliview, this will be curated and manually checked. Sure, one can write script to do this, but that would need time and also would require to check if it really does as it should.
Its just a suggestion, if you find script, you can use it anyway.
https://ormbunkar.se/aliview/
Hi Ana,
When merging single vcf files into a combined vcf you can set the genotype of a sample with no genotype to ref using VCFtools vcf-merge '-R' flag:
vcf-merge -R {0 for haploid} input1.vcf.gz input2.vcf.gz > merged_output.vcf
I have used this bash script to create a fasta from this output vcf:
for samp in input1 input2
do
printf '>'$samp'\n'
bcftools query -s $samp -f '[%TGT]' merged_output.vcf
printf '\n'
done > alignment.fa
I got this script from:
https://github.com/samtools/bcftools/issues/693
However, the output fasta will contain all nucleotides of the reference only if your vcf files have all positions, otherwise it will be an alignment of only variable sites. You'll need to make a vcf with all sites if you want it to be the same length and line up to the reference, but this may not be necessary depending on what you are doing.
1. if you've got Python interpreter installed, go to command line and type:
pip install sequman
2. using command line navigate to the directory where your fasta files are and type:
python -c "from sequman.sequman import fill_gaps; fill_gaps()"
this should create new files with nucleotides inserted from the reference sequence. these files will have trailing '_gaps_filled' sentence added to their original names.
the original files will remain intact.
note that the reference sequence for each file should be placed at the top of the file as the first record.
Many thanks Yuriy,
I'll try this last option since it seems to perfectly fit my data. I'll give feedback with my results.
Cheers!
This topic is pretty old, but just for those who may have same question and google to this topic:
I found this script by Stephan Kamrad, to be very useful:
https://github.com/Bahler-Lab/alignment-from-vcf
It's an old script, I believe from 2014, so it requires Python 2.x, in addition to pysam (I tried this release pysam-0.16.0.1-py27ha863e18_1), and biopython.
The good thing is that is deals very well with indels, which are included in the fasta output (so not just the SNPs, like some other scripts). The bad thing is that you have to specify a contig (you have to do it contig by contig) and you will get the whole contig sequence alignment not just the variable sites (but this latter could be preferable in some situations). The other thing that I wish it could do is to include the Ref sequence, but you can eventually make a fake VCF file for the Ref, which will have all 0, and merge it along with the samples.
As to the missing data (-) in the above question by Ana, the solution suggested by some participants above is enough to correct that. For example bcftools you can use --missing-to-ref option to fill those gaps, like what was described for vcftools above. Good luck...
- Badges
- Science topic
- Similar topics
- Phytochemistry
- Extracts