This is a test version of Biostars. For the public version, visit https://www.biostars.org.
BAM/SAM to FASTA conversion

Hi everone,

I searched Biostars for BAM/SAM to FASTA conversion method, and found the tools EMBOSS Picard could do this (Convert Bam File To Fasta File). I am wondering is there any perl script to accomplish this work? I should note that my seq data is strand-specific seq, so the command samtools view filename.bam | awk '{OFS="\t"; print ">"$1"\n"$10}' may not be very good.

I appreciate any of your solutions. Thanks a lot!

sam bam
• 89,294 views • link
updated by Ram Ram • written by biolab biolab

6 answers

For future record: samtools versions >1.3 can convert bam to fasta directly via samtools fasta

samtools --help

Program: samtools (Tools for alignments in the SAM format)
Version: 1.4 (using htslib 1.4)

Usage:   samtools <command> [options]

Commands:
  [...]
  -- File operations
     [...]
     fastq          converts a BAM to a FASTQ
     fasta          converts a BAM to a FASTA
  [...]
• 0 views • link
written by Carambakaracho Carambakaracho

The advice in the post you link to is greatly outdated. One of the reasons we should revisit these posts. There are far easier/faster ways to accomplish the same task.

The latest versions of samtools have the bam2fq command. Chain that to seqtk seq like so:

samtools bam2fq input.bam | seqtk seq -A > output.fa
• 0 views • link
updated by Ram Ram • written by Istvan Albert Istvan Albert

Hi Istvan, thank you very much for your answer. It's much helpful for my work. I just have two more minor questions: (1) my samtools version is Version: 0.1.8 (r613). I did not see samtools bam2fq function, so I need to update samtools? (2) Additionally is seqtk inside samtools? If not, where to get it? Thanks a lot!

• 0 views • link
updated by Ram Ram • written by biolab biolab

Hi Istvan, your answer is really helpful.

• 0 views • link
written by biolab biolab

You need to add '-' for seqtk in order to read from stdin :

samtools bam2fq input.bam | seqtk seq -A - > output.fa
• 0 views • link
updated by ATpoint ATpoint • written by Chadi Saad Chadi Saad

Please use ADD COMMENT or ADD REPLY to answer to previous reactions, as such this thread remains logically structured and easy to follow. I have now moved your post but as you can see it's not optimal. Adding an answer should only be used for providing a solution to the question asked.

• 0 views • link
written by WouterDeCoster WouterDeCoster

You can add - but it is not necessary.

• 0 views • link
written by ATpoint ATpoint

reformat.sh from BBMap suite will do the conversion efficiently.

reformat.sh in=your.bam out=seq.fa
reformat.sh in=your.sam out=seq.fa
• 0 views • link
written by GenoMax GenoMax

samtools fasta input.bam > output.fasta

• 0 views • link
written by kristina.mahan kristina.mahan

Via bam2bed and bed2faidxsta.pl:

$ bam2bed < reads.bam | bed2faidxsta.pl > reads.fa

Or sam2bed:

$ sam2bed < reads.sam | bed2faidxsta.pl > reads.fa
• 0 views • link
written by Alex Reynolds Alex Reynolds

I followed your instruction and even copied the unzipped bed2faidxstal.pl to the target folder, I got error "(base) bieniaszs-ipro:K20 bieniaszlab$ sam2bed < NL43_5562_filterredreads_header.sam | bed2faidxsta.pl > NL43_5562_filterredreads_header.fa -bash: bed2faidxsta.pl: command not found".

Then I add "perl" before the bed2faidxstal.pl, I still got error as below:

(base) bieniaszs-ipro:K20 bieniaszlab$ sam2bed < NL43_5562_filterredreads_header.sam | perl bed2faidxsta.pl > NL43_5562_filterredreads_header.fa
[E::fai_build3_core] Failed to open the file /Users/bieniaszlab/miniconda3/data/K_SRF11_U/K_20/K20/NL43_WT_CLIP.fa.gz
[faidx] Could not load fai index of /Users/bieniaszlab/miniconda3/data/K_SRF11_U/K_20/K20/NL43_WT_CLIP.fa.gz
[E::fai_build3_core] Failed to open the file /Users/bieniaszlab/miniconda3/data/K_SRF11_U/K_20/K20/NL43_WT_CLIP.fa.gz
[faidx] Could not load fai index of /Users/bieniaszlab/miniconda3/data/K_SRF11_U/K_20/K20/NL43_WT_CLIP.fa.gz
[E::fai_build3_core] Failed to open the file /Users/bieniaszlab/miniconda3/data/K_SRF11_U/K_20/K20/NL43_WT_CLIP.fa.gz
[faidx] Could not load fai index of /Users/bieniaszlab/miniconda3/data/K_SRF11_U/K_20/K20/NL43_WT_CLIP.fa.gz
[E::fai_build3_core] Failed to open the file /Users/bieniaszlab/miniconda3/data/K_SRF11_U/K_20/K20/NL43_WT_CLIP.fa.gz
• 0 views • link
updated by GenoMax GenoMax • written by xiaoleiusc xiaoleiusc

Put the directory containing the Perl script into your environment PATH, or run it with Perl. Make sure your FASTA is indexed with samtools.

• 0 views • link
written by Alex Reynolds Alex Reynolds

Thanks. I need FASTA output from SAM input, you mean to make sure my SAM file (not FASTA file) is indexed right?

• 0 views • link
written by xiaoleiusc xiaoleiusc

No, I mean that your reference genome FASTA must be indexed, in order to get out FASTA based on the BED coords that come out of the SAM conversion step. (Sorry for brevity, I was on a cell phone. Feel free to follow up if this isn't clear.)

• 0 views • link
written by Alex Reynolds Alex Reynolds

Log in to answer this question.

More posts like this