Is there a way to save a specific Gene in IGV of a whole exome sample?
2
0
Entering edit mode
6.4 years ago

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 • 1.7k views
ADD COMMENT
0
Entering edit mode
6.4 years ago

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.

ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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
ADD REPLY
0
Entering edit mode

@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.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

@ ATpoint. Thanks and updated the code.

ADD REPLY
0
Entering edit mode
6.4 years ago
sacha ★ 2.4k

Something like this ?

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

Login before adding your answer.

Traffic: 1295 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6