Pileup to consensus sequence aligned to reference - Please help

Hi everyone,

I am new to bioinformatics and have made some really good progress in the last few weeks however I have found myself a little stuck and would greatly appreciate any help from the galaxy community.

I have taken some illumina MiSeq Fastq files, trimmed and quality checked them before aligning to multiple custom reference genomes. I have then created a pileup from the BAM file using generate pileup in SAMTools and the same set of reference genomes. I then used pileup to interval to try to get a consensus sequence however because my sequence is incomplete I get multiple chains of sequence. Is there a tool which will allow me to get a consensus fasta file from this pileup with Ns aligned to the area of the reference where there is no sequence data available?

Thanks in advance


Hi, I am trying to solve the same issue. In my case, I created a pileup from BAM and used pileup-to-interval tool, but this gives me the reference sequence, not the consensus sequence. As I understood, I need to have a 10 column pileup (http://samtools.sourceforge.net/cns0.shtml), but I couldn’t figure out how to do this. Advanced pileup option gave me 8 columns, but this resulted an error when used for pileup-to-interval. Any help will be highly appreciated.


1 Like

Hi @kdcs & @emmad76!

Please see the choices in the BFCtools tool suite.

These tools will call variants (pileup or VCF), fill in reference bases where they are not represented in your data (a few different ways), and generate new consensus sequences given the 1) original reference sequence the variants were called against and the 2) variation output VCF. Plus many other manipulations. Please give these a try and see if it produces the output you each want – these are flexible tools with many options.

You will probably need to make use of a custom reference genome/transcriptome/exome fasta dataset. These tools do not have built-in indexes like mapping tools. If you are not sure where to find the fasta version of a pre-indexed reference genome you mapped against, please write back and we can help. No matter where you get it, it must be an exact match (genome build/source/version) for what you originally mapped against – plus the fasta should be in a very simple format – meaning, no “>” identifier line description content. The tool NormalizeFasta can be used in most cases to standardize the format of fasta datasets.

The “consensus sequence” that used to be generated by older versions of Mpileup were encoded and probably not what you are both wanting as a final result (is NOT a fasta “consensus sequence” result based on the variation in your data – what you might think of as a type of “assembly” result).

Also, using coordinates of regions in a pileup result (or VCF result, or gtf/bed/interval result) to Extract sequences from the genomic sequence will only result in fasta sequence based on that original reference genomic sequence again. It won’t including any base-level variation your read data may have had.

There isn’t a Galaxy Training Network tutorial that covers using these tools in detail, but looking at other workflows variant calling tutorials would probably help. Plus, you might want to compare tools/methods and compare. If interested, please see:


1 Like