Error when tabulating SAM file

[superbobry]

I'm using MethylCoder from GitHub with gsnap:

gsnap --quiet-if-excessive -A sam -k 15  --nofails --nthreads 24 -D /mnt/stripe -m 1 -B 5 -d human_g1k_v37 --mode=cmet-
nonstranded /mnt/stripe/GSM818003/SRR364089_1.fastq > /mnt/stripe/GSM818003/SRR364089_1_methylcoder/methylcoded.gsnap.s
am 2> /mnt/stripe/GSM818003/SRR364089_1_methylcoder/gsnap_run.log
^ executing gsnap. ^
tabulating methylation for /mnt/stripe/GSM818003/SRR364089_1_methylcoder/methylcoded.gsnap.sam
Traceback (most recent call last):
  File "/home/user/.virtualenvs/env/bin/methylcoder", line 9, in <module>
    load_entry_point('methylcoder==0.0', 'console_scripts', 'methylcoder')()
  File "/home/user/.virtualenvs/env/lib/python2.6/site-packages/methylcoder/__init__.py", line 732, in main
    opts.extra_args, opts.write_bin)
  File "/home/user/.virtualenvs/env/lib/python2.6/site-packages/methylcoder/gsnap.py", line 165, in main
    parse_gsnap_sam(gsnap_sam, ref_fasta, out_dir, paired_end, write_bin)
  File "/home/user/.virtualenvs/env/lib/python2.6/site-packages/methylcoder/gsnap.py", line 108, in parse_gsnap_sam
    sam_flag = int(line[1])
ValueError: invalid literal for int() with base 10: 'VN:1.0'

It seems that the parser isn't expecting SAM header. Have you considered using a separate library for reading and writing SAM files, i. e. pysam?

[brentp]

If you change that line to test .startswith('@') and open a pull request, I'll be glad to accept. I'd rather not depend on pysam.

[superbobry]

Yup, I've already tried this and ended up with a KeyError in pafasta.

gsnap-meth.py script seems to work however. Is there a reason for keeping it separately from the rest of the project?

[brentp]

Hmm, can you send the traceback for the KeyError after changing

re: gsnap-meth: I wrote that later and didn't test it fully. As far as I know it is working, but use with caution. I wanted to check it to make sure paired-end was working as expected. I had reports that it did but as I remember, there was a couple more cases I wanted to check.

-Brent

On Mon, Nov 11, 2013 at 7:30 AM, Sergei Lebedev notifications@github.comwrote:

Yup, I've already tried this and ended up with a KeyError in pafasta.

gsnap-meth.py script seems to work however. Is there a reason for keeping it separately from the rest of the project?

— Reply to this email directly or view it on GitHubhttps://github.com/brentp/methylcode/issues/7#issuecomment-28203809 .