Is there a way to save a specific Gene in IGV of a whole exome sample?

Question

[mhmtgenc85]

Dear all, Is there a way to save a specific Gene in IGV of a whole exome sample? For example I have the whole exome sequencing of a sample but I would like to view and save only one gene out of that whole BAM as only one gene specific BAM?

Thanks in advance

Answer 1

[cpad0112]

Make a bam for gene of interest (using gene coordinates), index the bam and load the indexed bam in IGV.

samtools view -h test.bam gene_coordinates
samtools view -h aln.sorted.bam chr2:20,100,000-20,200,000

Extract Alignment From Very Large Bam File: Extract Alignment From Very Large Bam File

or you can write a session file in IGV, so that IGV always zooms in to that gene of interest.

[mhmtgenc85]

Thank you for the tip. I tried it and get the BAM file. But now I need to get the index file of the bam. When I use

samtools index test.bam > test.bai I get an error like this

EOF marker is absent. invalid BAM binary header (this is not a BAM file) Invalid BAM header. Fail to index..

I have tried this command and it worked.

samtools view -b test.bam > output.bam

Changed -h to -b ...

Thank you.

[GenoMax]

Are you using latest version of samtools? If so, command to index is samtools index test.bam. You should also sort the test.bam file to be safe by doing samtools sort -o test_sorted.bam test.bam.

[cpad0112]

Probably bam headers are not extracted.

@OP: Try this instead:

samtools view -b aln.sorted.bam chr2:20,100,000-20,200,000 > out.bam
samtools index out.bam

[cpad0112]

@OP: Btw, you do not have to output bai file like samtools index test.bam > test.bai. Run samtools index test.bam. That would create .bai file automatically.

[cpad0112]

if you are on ubuntu, could you please post the output from test.bam:

$ od -x test.bam | head -1

It should be 0000000 8b1f 0408 0000 0000 ff00 0006 4342 0002

[ATpoint]

Command for bam files is -b (implies -h).

[cpad0112]

@ ATpoint. Thanks and updated the code.

Answer 2

[sacha]

Something like this ?

bedtools intersect -a exome.bam -b mygene.bed  > mygene.bam