How to call variants from a whole genome alignment file

[baika]

Hi All
I have a quick question about calling SNP/InDel from a whole genome alignment file which is in FASTA format. I have aligned 45 assembled bacterial genomes (~4Mb each) using Mugsy tool and converted the MAF output to FASTA. One of the 45 genomes is the reference. Now I want a table (preferably in VCF format) showing SNP/InDels in 44 strains as compared to my reference strain. Please suggest how to do this?

Thanks

baika

[arkal]

I think if u use the samtools mpileup command using the reference genome of ur choice as the input reference fasta and the remaining 44 alignments in the bam format u shud get what u want! Correct me if im wrong!

[baika]

Is it possible to convert an aligned FASTA file into a BAM file? I would really appreciate if you could suggest me any tool to do that.

Thanks

[brofallon]

I don't really think its possible to convert a fasta file directly into BAM. BAM files are meant to store many short reads with associated quality scores, but fasta is just a listing of a single sequence. You may be able to use some trickery to force the fasta into BAM (or, more likely, BAMs uncompressed version, SAM), but I'm guessing this would be more trouble than its worth.
If you want to convert the aligned, fasta-formatted genomes into a vcf, I bet you'll end up writing a script (bash, perl, python, etc.) to do the job. VCF format isn't super complicated, and the script would simply look at each alignment column and see which samples differed from the reference. That would be my advice...
HTH

[arkal]

Is it possible to convert an aligned FASTA file into a BAM file? I would really appreciate if you could suggest me any tool to do that.

Thanks

Sorry there my mistake... when you said you aligned genomes i automatically thought in terms of fastq onto ref fast alignment stored in the form of bams.
So in my opinion, what you CAN do instead is
1. Identify your reference bacterial genome fasta and store it in a separate file.
2. Download a read simulator/generator (eg DWGsim, ART, Maq, etc)
3. generate PAIRED END reads for the remaining 44 genomes. ENSURE THAT SNP AND INDEL INTRODUCTION RATES ARE SET TO ZERO.
4. Align these reads to the reference separately using an aligner such as BWA, NOVOAlign or Stampy, etc.
5. Do a samtools mpileup to get your results in a bcf format!
6. Thank me later

[aniruddha.otago]

I have analysed my NGS data from Illumina.. (RRBS).. I have fragments (chromosome, genomic position etc).. I would like to search whether there is any common SNPs present in my fragments. I mean I would like to search against databases and see what are the chances that these fragments could contain a potential common SNPS. I presume I would like to DbSNPs (131 or 135 may be). But if anyone details the process or advices which will enable me to do this search quickly that will be much appreciated.

[DR.AYAH]

How to call variants from a whole genome alignment file

Hi All
I have a quick question about calling SNP/InDel from a whole genome alignment file which is in FASTA format. I have aligned 45 assembled bacterial genomes (~4Mb each) using Mugsy tool and converted the MAF output to FASTA. One of the 45 genomes is the reference. Now I want a table (preferably in VCF format) showing SNP/InDels in 44 strains as compared to my reference strain. Please suggest how to do this?

Thanks

baika

i am interested to call variants from a whole-genome alignment file, i am working with fungi have 6 whole genomes one of them is reference strain (~13Mb for each).

i read more about Mugsy software, i am new for using Linux and dealing with Terminal command.

could you help to know the command to run the alignment for 6 genomes?

[salarshaaf]

Have a look on TASSEL. It can convert your fasta file into VCF without writing commands.
here is the download link:

https://www.maizegenetics.net/tassel