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?
0
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
1
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.)
-
$\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$– rimoCommented Nov 28, 2023 at 18:11