This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Is there a way to save a specific Gene in IGV of a whole exome sample?

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

bam igv sequence
• 2,567 views • link
updated by sacha sacha • written by Clark_BioMorgan • Clark_BioMorgan

2 answers

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.

• 1 views • link
written by cpad0112 cpad0112

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.

• 1 views • link
written by Clark_BioMorgan • Clark_BioMorgan

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.

• 1 views • link
written by GenoMax GenoMax

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
• 1 views • link
written by cpad0112 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.

• 1 views • link
written by cpad0112 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

• 1 views • link
written by cpad0112 cpad0112

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

• 1 views • link
written by ATpoint ATpoint

@ ATpoint. Thanks and updated the code.

• 1 views • link
written by cpad0112 cpad0112

Something like this ? 

bedtools intersect -a exome.bam -b mygene.bed  > mygene.bam
• 1 views • link
written by sacha sacha

Log in to answer this question.

More posts like this