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