How to create a consensus of a contig with samtools or bbmap?

Question

[marongiu.luigi]

I have assembled my reads using samtools and obtained a contig.fasta files that I use to remap with BBmap to generate an alignment following these indications. The point is to generate the genome of a viral isolate.

IGV shows that there are 3 main contigs in the obtained bam file:

So I have extracted these contigs with samtools as follows:

samtools view assembled.bam "NODE_1_length_152472_cov_5.695244_cutoff_5_type_linears" > contig_1.bam

I then wanted to generate a consensus of this alignment using again samtools but I got an error:

samtools consensus -f fasta contig_1.bam -o cont_1.fa 
[E::sam_parse1] no SQ lines present in the header
bam_next_seq() failure.
samtools consensus: failed

Why is there a SQ line missing? The head of the bam file shows:

head phiM/Mapped/contig_1.bam 
LH00452:285:232WNLLT4:7:2479:15335:23308 1:N:0:AACAGCCAGT+GAGTGCTGTT    99  NODE_1_length_152472_cov_5.695244_cutoff_5_type_linears 1   45  150=    =122    271 AGAATAACAACGCCATTAAAATAACCACAAACATTAATGTCACAATAATCATTAAAAATAACCAACAAACGATAGTCAAAATATCTTGATAAAAAAAGAACGCAATCACTTTTAGTATAATAGCCGCTACCTCTTTTCCTCTAAAGAATA  IIIII9IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII9IIIIIIIIIIIIIIIIIIII9IIIIIIIIIIIIIIIIIIIIIIII9IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII9IIIII  NM:i:0  AM:i:45
LH00452:285:232WNLLT4:7:2263:22236:23153 1:N:0:AACAGCCAGT+GAGTGCTGTT    163 NODE_1_length_152472_cov_5.695244_cutoff_5_type_linears 3   45  150=    =169    316 AATAACAACGCCATTAAAATAACCACAAACATTAATGTCACAATAATCATTAAAAATAACCAACAAACGATAGTCAAAATATCTTGATAAAAAAAGAACGCAATCACTTTTAGTATAATAGCCGCTACCTCTTTTCCTCTAAAGAATAAA  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII  NM:i:0  AM:i:45

...

Alternatively, BBmap also allows the creation of a consensus: consensus.sh in=infile.bam ref=ref.fa out=cons.fa. Is it possible to run it without a reference? If I run it without the ref argument, I of course get the error:

Exception in thread "main" java.lang.RuntimeException: Error - a reference file is required.

Thank you.

Answer 1

[lieven.sterck]

You assembled a genome using samtools ???? that sounds very weird, if not impossible ...

The header (SQ lines) is missing from your file because you did not ask to extract it. use 'view -h' (== the -h is the key parameter here) , or alternatively extract only the header (that's done with the -H parameter) and paste it on top of your sam file.

BBmap is an aligner tool so yes it will always need some kind of reference to map against (as the tool error message clearly says ;-) )