Answers.Body: <p>The first approach that comes to mind is to cut the bam file and then convert the already cut bam to fasta:</p>
  <pre><code>samtools view -hb chr1:10100500-10200500 &gt; small.bam
  samtools fasta small.bam &gt; small.fa
  </code></pre>
  <p>(note that I am not certain about the <code>samtools fasta</code> command, I have never used it and I don't know if it will produce many sequences, one for each haplotype supported by reads; please test this.)</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-11-28T18:11:43.147Z
Answers.Text: https://www.htslib.org/doc/samtools-fasta.html -- the documentation for `samtools fasta`it should work as @terdon wrote it but just in case you need more information
Answers.UserDisplayName: rimo
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-11-19T12:01:35.780Z
Answers.LastActivityDate: 2023-11-20T10:56:45.683Z
Answers.LastEditDate: 2023-11-20T10:56:45.683Z
Answers.UserDisplayName: terdon
Body: <p>I need to get a fasta sequence of a certain gene for a certain worm strain that is different from reference. I have a reference genome, BAM for the strain of interest, and coordinates of the gene. I know that vcftools can convert bam to fasta, but I do not see if I can do it only for a certain coordinates (like, for example, chr1:10100500-10200500). I converted the whole bam to fasta and then I tried using bedtools getfasta, but my system ran out of space. Is there a more rational way to do it?</p>
  
Tags: <fasta><samtools><bedtools>
Title: Get a certain gene sequence from bam/vcf and reference
ThreadText: # Get a certain gene sequence from bam/vcf and reference
  
  (Asked by: username_0 on 2023-11-17)
  
  I need to get a fasta sequence of a certain gene for a certain worm strain that is different from reference. I have a reference genome, BAM for the strain of interest, and coordinates of the gene. I know that vcftools can convert bam to fasta, but I do not see if I can do it only for a certain coordinates (like, for example, chr1:10100500-10200500). I converted the whole bam to fasta and then I tried using bedtools getfasta, but my system ran out of space. Is there a more rational way to do it?
  
  **Tags:** fasta, samtools, bedtools
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2023-11-19. Score: 4
  
  The first approach that comes to mind is to cut the bam file and then convert the already cut bam to fasta:
  
  ```
  samtools view -hb chr1:10100500-10200500 > small.bam
  samtools fasta small.bam > small.fa
  
  ```
  
  (note that I am not certain about the `samtools fasta` command, I have never used it and I don't know if it will produce many sequences, one for each haplotype supported by reads; please test this.)
  
  > **Comments:**
  > **username_2**: https://www.htslib.org/doc/samtools-fasta.html -- the documentation for `samtools fasta`it should work as @username_1 wrote it but just in case you need more information  
  
  ---
  
  