Body: <p>Let's assume I have a BAM file and several positions that I would like to examine more closely in this alignment. My goal is to find out whether these positions are on the same reads and which nucleotides each read has in these two positions.
  So my output should be a table of all reads that occur in the positions and list the associated nucleotides.</p>
  <p>The output should look like this (with 980, 1350 and 1400 the first three positions; I made up the read names):</p>
  <div class="s-table-container">
  <table class="s-table">
  <thead>
  <tr>
  <th>read name</th>
  <th>980</th>
  <th>1350</th>
  <th>1400</th>
  <th>...</th>
  </tr>
  </thead>
  <tbody>
  <tr>
  <td>ReAdN4m1L03</td>
  <td>A</td>
  <td>T</td>
  <td></td>
  <td>...</td>
  </tr>
  <tr>
  <td>ReAdN4m2L02</td>
  <td>A</td>
  <td>T</td>
  <td>G</td>
  <td>...</td>
  </tr>
  <tr>
  <td>ReAdN4m3L01</td>
  <td></td>
  <td>T</td>
  <td>G</td>
  <td>...</td>
  </tr>
  <tr>
  <td>...</td>
  <td>...</td>
  <td>...</td>
  <td>...</td>
  <td>...</td>
  </tr>
  </tbody>
  </table>
  </div>
  <p>At the end I would like to have a table that links the nucleotides via the reads.</p>
  <p>So far, I have tried to read in the Bam file at the desired positions:</p>
  <pre><code>#Positions of interest
  pos &lt;- c(980, 1350, 1400)
  
  #path to BAM file
  path &lt;- pathToBamFile
  
  #name of reference used for mapping
  refname &lt;- &quot;myreference&quot;
  
  #loading the bam file
  bamFile &lt;- BamFile(path)
  gr &lt;- GRanges(seqnames = refname, ranges = IRanges(start = pos, end = pos))
  params &lt;- ScanBamParam(which = gr, what = scanBamWhat())
  aln &lt;- scanBam(bamFile, param = params)
  </code></pre>
  <p>With the following code I can access the reads, but how do I get the corresponding nucleotides in the individual positions?</p>
  <pre><code>pos.no &lt;- c(1:length(pos))
  
  i &lt;- 1
  aln[[pos.no[i]]]$qname
  </code></pre>
  <p>I have already tried to work with the CIGAR strings <code>aln[[i]]$cigar</code> with the sequences <code>aln[[i]]$seq</code> in order to extract to the nucleotides manually, but that led to more and more problems.</p>
  <p>Is there any easy way to extract the nucleotides?</p>
  <p>I have searched for a long time now on the internet for an easy solution but without success. Hopefully someone will help me here! Many thanks for your time in advance.</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-08-21T07:41:41.543Z
Comments.Text: Do you have to use R? This is more convenient if you use the pileup interface from pysam (perhaps Rsamtools provides a similar functionality).
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-08-22T06:35:00.410Z
Comments.Text: @DevonRyan but does pileup give you an output where the nucleotides are linked to the reads?
Comments.UserDisplayName: wejt
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-08-22T08:39:32.973Z
Comments.Text: Yes, it does, at least in non-R interfaces (i.e., python and C).
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-08-29T22:18:24.600Z
Comments.Text: This answer does it in Python, perhaps there are some ideas there: https://bioinformatics.stackexchange.com/questions/7401/access-base-aligned-to-particular-reference-position
Comments.UserDisplayName: Maximilian Press
Tags: <r><sam><snp><reads><cigar>
Title: Extract read names and the associated nucleotides on specific positions from a BAM file (in R)
ThreadText: # Extract read names and the associated nucleotides on specific positions from a BAM file (in R)
  
  (Asked by: username_0 on 2021-08-20)
  
  Let's assume I have a BAM file and several positions that I would like to examine more closely in this alignment. My goal is to find out whether these positions are on the same reads and which nucleotides each read has in these two positions.
  So my output should be a table of all reads that occur in the positions and list the associated nucleotides.
  
  The output should look like this (with 980, 1350 and 1400 the first three positions; I made up the read names):
  
  | read name | 980 | 1350 | 1400 | ... |
  | --- | --- | --- | --- | --- |
  | ReAdN4m1L03 | A | T |  | ... |
  | ReAdN4m2L02 | A | T | G | ... |
  | ReAdN4m3L01 |  | T | G | ... |
  | ... | ... | ... | ... | ... |
  
  At the end I would like to have a table that links the nucleotides via the reads.
  
  So far, I have tried to read in the Bam file at the desired positions:
  
  ```
  #Positions of interest
  pos <- c(980, 1350, 1400)
  
  #path to BAM file
  path <- pathToBamFile
  
  #name of reference used for mapping
  refname <- "myreference"
  
  #loading the bam file
  bamFile <- BamFile(path)
  gr <- GRanges(seqnames = refname, ranges = IRanges(start = pos, end = pos))
  params <- ScanBamParam(which = gr, what = scanBamWhat())
  aln <- scanBam(bamFile, param = params)
  
  ```
  
  With the following code I can access the reads, but how do I get the corresponding nucleotides in the individual positions?
  
  ```
  pos.no <- c(1:length(pos))
  
  i <- 1
  aln[[pos.no[i]]]$qname
  
  ```
  
  I have already tried to work with the CIGAR strings `aln[[i]]$cigar` with the sequences `aln[[i]]$seq` in order to extract to the nucleotides manually, but that led to more and more problems.
  
  Is there any easy way to extract the nucleotides?
  
  I have searched for a long time now on the internet for an easy solution but without success. Hopefully someone will help me here! Many thanks for your time in advance.
  
  **Tags:** r, sam, snp, reads, cigar
  
  > **Comments:**
  > **username_1**: Do you have to use R? This is more convenient if you use the pileup interface from pysam (perhaps Rsamtools provides a similar functionality).  
  > **username_0**: @username_1 but does pileup give you an output where the nucleotides are linked to the reads?  
  > **username_1**: Yes, it does, at least in non-R interfaces (i.e., python and C).  
  > **username_2**: This answer does it in Python, perhaps there are some ideas there: https://bioinformatics.stackexchange.com/questions/7401/access-base-aligned-to-particular-reference-position  
  
  ---
  
  ## 0 answers
  
  