Recent changes to 15: read.get_sam_line() does not reflect strand of read intervalhttps://sourceforge.net/p/htseq/bugs/15/2013-08-02T21:49:40.779000ZRecent changes to 15: read.get_sam_line() does not reflect strand of read intervalread.get_sam_line() does not reflect strand of read interval2013-08-02T21:49:40.779000Z2013-08-02T21:49:40.779000ZDavid Koppsteinhttps://sourceforge.net/u/koppology/https://sourceforge.netd12c10666ec570ba7e431c73aba0d1edde968c08<div class="markdown_content"><div class="codehilite"><pre><span class="o">>>></span> <span class="n">import</span> <span class="n">HTSeq</span>
<span class="o">>>></span> <span class="n">reader</span> <span class="o">=</span> <span class="n">HTSeq</span><span class="p">.</span><span class="n">BAM_Reader</span><span class="p">(</span><span class="err">'</span><span class="n">test</span><span class="p">.</span><span class="n">bam</span><span class="err">'</span><span class="p">)</span>
<span class="o">>>></span> <span class="k">for</span> <span class="n">read</span> <span class="n">in</span> <span class="n">reader</span><span class="o">:</span>
<span class="p">...</span> <span class="n">pass</span>
<span class="p">...</span>
<span class="o">>>></span> <span class="n">read</span>
<span class="o"><</span><span class="n">SAM_Alignment</span> <span class="n">object</span><span class="o">:</span> <span class="n">Paired</span><span class="o">-</span><span class="n">end</span> <span class="n">Read</span> <span class="err">'</span><span class="n">WIGTC</span><span class="o">-</span><span class="n">HISEQ</span><span class="o">:</span><span class="mi">8</span><span class="o">:</span><span class="mi">2204</span><span class="o">:</span><span class="mi">1251</span><span class="o">:</span><span class="mi">1976</span><span class="err">#</span><span class="n">TATAAT</span><span class="o">/</span><span class="mi">1</span><span class="p">;</span><span class="mi">1</span><span class="o">:</span><span class="n">count</span><span class="o">=</span><span class="mi">1</span><span class="err">'</span> <span class="n">aligned</span> <span class="n">to</span> <span class="n">chr7</span><span class="o">:</span><span class="p">[</span><span class="mi">5570231</span><span class="p">,</span><span class="mi">5570247</span><span class="p">)</span><span class="o">/+></span>
</pre></div>
<h1 id="note-that-the-strand-is">Note that the strand is +</h1>
<div class="codehilite"><pre><span class="o">>>></span> <span class="n">read</span><span class="p">.</span><span class="n">iv</span>
<span class="o"><</span><span class="n">GenomicInterval</span> <span class="n">object</span> <span class="err">'</span><span class="n">chr7</span><span class="err">'</span><span class="p">,</span> <span class="p">[</span><span class="mi">5570231</span><span class="p">,</span><span class="mi">5570247</span><span class="p">),</span> <span class="n">strand</span> <span class="sc">'+'</span><span class="o">></span>
</pre></div>
<h1 id="correspondingly-the-sam-code-is-0">Correspondingly, the SAM code is 0</h1>
<div class="codehilite"><pre><span class="o">>>></span> <span class="n">read</span><span class="p">.</span><span class="n">get_sam_line</span><span class="p">()</span>
</pre></div>
<p>'WIGTC-HISEQ:8:2204:1251:1976#TATAAT/1;1:count=1\t<strong><em>0</em></strong>\tchr7\t<strong><em>5570232</em></strong>\t255\t16M\t*\t0\t0\tACCGCCGAGACCGCGT\tgggggiiiiegiihii\tXA:i:0'</p>
<h1 id="giving-a-new-genomic-interval-changes-the-output-of-the-coordinates">Giving a new Genomic Interval changes the output of the coordinates</h1>
<div class="codehilite"><pre><span class="o">>>></span> <span class="n">read</span><span class="p">.</span><span class="n">iv</span> <span class="o">=</span> <span class="n">HTSeq</span><span class="p">.</span><span class="n">GenomicInterval</span><span class="p">(</span><span class="err">'</span><span class="n">chr7</span><span class="err">'</span><span class="p">,</span> <span class="mi">5570232</span><span class="p">,</span> <span class="mi">5570248</span><span class="p">,</span> <span class="sc">'+'</span><span class="p">)</span>
<span class="o">>>></span> <span class="n">read</span><span class="p">.</span><span class="n">get_sam_line</span><span class="p">()</span>
</pre></div>
<p>'WIGTC-HISEQ:8:2204:1251:1976#TATAAT/1;1:count=1\t<strong><em>0</em></strong>\tchr7\t<strong><em>5570233</em></strong>\t255\t16M\t*\t0\t0\tACCGCCGAGACCGCGT\tgggggiiiiegiihii\tXA:i:0'</p>
<h1 id="but-changing-the-strand-of-the-genomic-interval-does-not-change-the-sam-code-to-16-the-reverse-complement-sam-code-which-i-think-it-should">But changing the strand of the genomic interval does not change the SAM code to 16 (the "reverse complement" SAM code, which I think it should)</h1>
<div class="codehilite"><pre><span class="o">>>></span> <span class="n">read</span><span class="p">.</span><span class="n">iv</span> <span class="o">=</span> <span class="n">HTSeq</span><span class="p">.</span><span class="n">GenomicInterval</span><span class="p">(</span><span class="err">'</span><span class="n">chr7</span><span class="err">'</span><span class="p">,</span> <span class="mi">5570232</span><span class="p">,</span> <span class="mi">5570248</span><span class="p">,</span> <span class="sc">'-'</span><span class="p">)</span>
<span class="o">>>></span> <span class="n">read</span><span class="p">.</span><span class="n">get_sam_line</span><span class="p">()</span>
</pre></div>
<p>'WIGTC-HISEQ:8:2204:1251:1976#TATAAT/1;1:count=1\t<strong><em>0</em></strong>\tchr7\t5570233\t255\t16M\t*\t0\t0\tACCGCCGAGACCGCGT\tgggggiiiiegiihii\tXA:i:0'</p></div>