Hello everyone,
I am having problems extracting the mapped reads at a certain position from a BAM alignment file. I hope that you can help me.
I have paired end RNA-seq reads which I mapped to the human reference genome using TopHat. What I want to do now, is for a predefined list of SNP positions, extract the reads that mapped to each position. I am currently using samtools view for this, for instance:
For most of the positions this works fine, however for some it goes wrong. It happens when the SNP position is in the insert region between two paired reads, then samtools returns the reads of those pairs, although at the exact location no reads were mapped. What I want is just an empty file for these locations.
Is there any filtering that I can apply, so that it only returns reads that acually map at the given location? Or is there another tool that I can use for this? I tried bedtools (intersectBed), but it acts the same as samtools view. Any help is appreciated.
Thanx in advance!
I am having problems extracting the mapped reads at a certain position from a BAM alignment file. I hope that you can help me.
I have paired end RNA-seq reads which I mapped to the human reference genome using TopHat. What I want to do now, is for a predefined list of SNP positions, extract the reads that mapped to each position. I am currently using samtools view for this, for instance:
Code:
samtools view reads_unique_hits.sorted.bam chr2:96858914-96858914 > extracted_reads.sam
Is there any filtering that I can apply, so that it only returns reads that acually map at the given location? Or is there another tool that I can use for this? I tried bedtools (intersectBed), but it acts the same as samtools view. Any help is appreciated.
Thanx in advance!
Comment