1
$\begingroup$

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?

$\endgroup$
0

1 Answer 1

4
$\begingroup$

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.)

$\endgroup$
1
  • $\begingroup$ htslib.org/doc/samtools-fasta.html -- the documentation for samtools fastait should work as @terdon wrote it but just in case you need more information $\endgroup$
    – rimo
    Nov 28, 2023 at 18:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.