<div dir="ltr">Sorry, in plain English, my previous reply amounts to this:<div><br></div><div style>Your trimmed sequence only had 5 k-mers of length 32.</div><div style>Your full sequence had 43 k-mers of length 32.</div>
<div style>Since the median samples from the middle of the sorted counts, you will not see a count of 1 unless at least half of the counts are 1.</div><div style>To achieve this, you need a trimmed sequence of at least length 53.</div>
<div style><br></div><div style>Make sense?</div><div style><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Fri, Mar 15, 2013 at 6:09 PM, Eric McDonald <span dir="ltr"><<a href="mailto:emcd.msu@gmail.com" target="_blank">emcd.msu@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Yes.<div><br></div><div>Default k = 32.</div><div>Full sequence length is 75. 75 - 32 = 43.</div><div>Trimmed sequence length is 37. 37 -32 = 5.</div>
<div>To get a median of >0, you need a trimmed sequence length of at least 53, because 43 / 2 = 21, 21 - 5 = 16, 37 + 16 = 53.</div>
<div><br></div><div>For example, try this as your sequence in test1.fa:</div><div> AAGACACCCTCACCCCTAGCTGCGCGAGGCCCTCTCCCCTGGGTAGAGGGTCA</div><div><br></div><div>Nothing wrong with the code - just your math skills, Adina. :-)</div>
<div><br></div><div>Hope that helps,</div><div> Eric</div><div><br></div></div><div class="gmail_extra"><div><div class="h5"><br><br><div class="gmail_quote">On Wed, Mar 13, 2013 at 12:03 AM, Adina Chuang Howe <span dir="ltr"><<a href="mailto:adina.chuang@gmail.com" target="_blank">adina.chuang@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">It is definitely reproducible. And I have at least one sequence that can be identified that is causing this.<div><br>
</div>
<div>Here's what I'm seeing: (these are on the same HPC scratch space as before)</div><div>
<br></div><div>test1.fa contains only (37 bp):</div><div><br></div><div><div>>SRR172902.316476<span style="white-space:pre-wrap"> </span>42391</div><div>AAGACACCCTCACCCCTAGCTGCGCGAGGCCCTCTCC</div></div>
<div><br></div><div>SRR172902.316476.single.fq contains the untrimmed same:</div><div><br></div><div><div>@SRR172902.316476 USI-EAS376:1:5:311:233 length=75</div><div>AAGACACCCTCACCCCTAGCTGCGCGAGGCCCTCTCCCCTGGGTAGAGGGTCAAACAGCGCAAGGCAACAGATCG</div>
<div>+SRR172902.316476 USI-EAS376:1:5:311:233 length=75</div><div>BBBBBABBB@BBBBBB>B>B=AA>B?A=4ABA;B>>AAA79;;068799;2;====>3>9>7922=;739#####</div><div><br></div><div class="gmail_quote">python sweep-reads3.py test1.fa SRR172902.316476.single.fq </div>
<div class="gmail_quote"><br></div><div class="gmail_quote">results in an empty test1.fa.sweep</div><div class="gmail_quote"><br></div><div class="gmail_quote">But...</div><div class="gmail_quote"><br></div><div class="gmail_quote">
python sweep-reads3.py SRR172902.316476.single.fq SRR172902.316476.single.fq</div><div class="gmail_quote"><br></div><div class="gmail_quote">Results in</div><div class="gmail_quote"><br></div><div class="gmail_quote"><div class="gmail_quote">
==> SRR172902.316476.single.fq.sweep3 <==</div><div class="gmail_quote">>SRR172902.316476</div><div class="gmail_quote">AAGACACCCTCACCCCTAGCTGCGCGAGGCCCTCTCCCCTGGGTAGAGGGTCAAACAGCGCAAGGCAACAGATCG</div></div><div class="gmail_quote">
<br></div><div class="gmail_quote">Any clue?</div><div><div><div class="gmail_quote"><br></div><div class="gmail_quote"><br></div><div class="gmail_quote">On Sun, Mar 10, 2013 at 11:33 PM, C. Titus Brown <span dir="ltr"><<a href="mailto:ctb@msu.edu" target="_blank">ctb@msu.edu</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div>On Thu, Mar 07, 2013 at 02:10:17PM -0500, Adina Chuang Howe wrote:<br>
> Possible bug in sweep-reads...I'm not recovering the partitioned reads<br>
> from the original dataset.<br>
><br>
> First observed this when I looked at lotsa partitions and trying to<br>
> recover swept reads - some swept files would show up empty:<br>
><br>
> command:<br>
> python sweep-reads3.py -N 4 -k 32 -x 1e9<br>
> /mnt/research/gpgc/hmp-mock-partitions/001264-files/no-sweep-pids/pid*fa<br>
> /mnt/research/gpgc/hmp-mock-partitions/SRR-combined.fastq<br>
><br>
> troubleshooting:<br>
> Then I looked at just one partition:<br>
> on HPC: /mnt/scratch/howead/test<br>
> python sweep-reads3.py pid-42391.fa SRR-combined.fastq<br>
><br>
> And resulting sweepfile is empty.<br>
><br>
> If I run:<br>
> python sweep-reads3.py pid-42391.fa pid-42391.fa<br>
><br>
> Behavior is correct.<br>
<br>
</div></div>Very weird. I don't see anything obviously wrong in the script, which<br>
just means it's a subtle and deep bug.<br>
<br>
Two questions --<br>
<br>
- is it reproducible, i.e. do you get the same results every time you<br>
run it? (please yes)<br>
<br>
- can you break it down to a smaller failure point than with SRR-combined,<br>
e.g. maybe a few hundred k reads?<br>
<div><div><br>
thanks,<br>
--titus<br>
--<br>
C. Titus Brown, <a href="mailto:ctb@msu.edu" target="_blank">ctb@msu.edu</a><br>
</div></div></blockquote></div><br></div></div></div>
</blockquote></div><br><br clear="all"><div><br></div></div></div><div class="im">-- <br><div dir="ltr"><div>Eric McDonald</div><div>HPC/Cloud Software Engineer</div><div> for the Institute for Cyber-Enabled Research (iCER)</div>
<div> and the Laboratory for Genomics, Evolution, and Development (GED)</div>
<div>Michigan State University</div><div>P: <a href="tel:517-355-8733" value="+15173558733" target="_blank">517-355-8733</a></div></div>
</div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div dir="ltr"><div>Eric McDonald</div><div>HPC/Cloud Software Engineer</div><div> for the Institute for Cyber-Enabled Research (iCER)</div><div> and the Laboratory for Genomics, Evolution, and Development (GED)</div>
<div>Michigan State University</div><div>P: 517-355-8733</div></div>
</div>