Hello,
I have an alignment bam file, and I would like to retrieve all the bases aligning to that specific location in the genome. For example, if I am interested in position X in the genome, I would like to extract all the nucleotides that are aligning to that position X from my alignment file
•
written
by
rsafavi •
If you mean, you want to extract read base which is aligning to particular position of your genome, then you can try samtools mpileup or gatk Pileup. You can explore more about mpileup format here. Samtools mpileup will give you output in form of 5 columns :
From this you can idea about the base present in reference genome and nucleotide of the reads mapping to that particular position.
If you have just handful of positions to check then you can also view your alignment in genome viewer like IGV or Tablet to visualise the bases at that particular position.
I am working with Nanopore cDNA data. I tried samtool mpileup before like this:
samtools mpileup intersectedReads.bam -r chr1:13625035-13625035 > intersectedReads.pileup
as a test case, but I got:
chr1 13625035 N 0 * *
chr1 13625036 N 0 * *
which is not correct, since when I look at the genome browser I see two reads aligning to that region
Posted an answer here: A: Extract Base Based On Position From Bam File
Kevin
Log in to answer this question.
Hello rsafavi,
why do you like to do this and how should your desired output look like?
The procedure would be to get all reads that overlap the desired position using
samtools view. Due to the mapping position given for each read, one can than find the base on the location with respect to the CIGAR value.fin swimmer