Question: Is there a way to save a specific Gene in IGV of a whole exome sample?
0
gravatar for mhmtgenc85
2.0 years ago by
mhmtgenc8540
Turkey
mhmtgenc8540 wrote:

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

igv bam sequence • 546 views
ADD COMMENTlink modified 2.0 years ago by sacha1.9k • written 2.0 years ago by mhmtgenc8540
0
gravatar for cpad0112
2.0 years ago by
cpad011213k
India
cpad011213k wrote:

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 COMMENTlink modified 2.0 years ago • written 2.0 years ago by cpad011213k

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 REPLYlink modified 2.0 years ago • written 2.0 years ago by mhmtgenc8540

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 REPLYlink written 2.0 years ago by genomax83k

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 REPLYlink modified 2.0 years ago • written 2.0 years ago by cpad011213k

@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 REPLYlink written 2.0 years ago by cpad011213k

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 REPLYlink written 2.0 years ago by cpad011213k

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

ADD REPLYlink written 2.0 years ago by ATpoint34k

@ ATpoint. Thanks and updated the code.

ADD REPLYlink written 2.0 years ago by cpad011213k
0
gravatar for sacha
2.0 years ago by
sacha1.9k
France
sacha1.9k wrote:

Something like this ?

bedtools intersect -a exome.bam -b mygene.bed  > mygene.bam
ADD COMMENTlink written 2.0 years ago by sacha1.9k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1591 users visited in the last hour