I have RNAseq data and aligned the fastq files to the genome. so now I have bam file. I want to make the bam file only for one gene for example gapdh. to do so at first I have converted the bam file to sam file to be able to grep it. in the sam file we only have coordinates but when I tried to look for the coordinates of gapdh I did not find them in the sam file. do you know how I can make a sam or bam file only for the gapdh?
•
written
by
alireza346 •
Get the coordinates of your gene of interest and use samtools view yourfile.bam chrN:xxxxxxx-xxxxxxxx -o GAPDH.bam
Hello,
there is no need to convert bam to sam to use grep. If you use samtools view samtools is doing this for you.
But at all you don't need grep here. All you have to know is, against which reference genome you have mapped your fastq data and how your coordinates are for the desired gene. Let's assume you used hg38. Than the coordinates for GAPDH are 12:6533927-6538374.
With this informations you can create a subset of your bam file:
$ samtools view -b -o output.bam input.bam 12:6533927-6538374
fin swimmer
AFAIK in recent versions the -b is not necessary because samtools figures out the output format based on the extension specified in -o ;-)
Log in to answer this question.
