comparison tools/fastq_paired_unpaired/fastq_paired_unpaired.py @ 6:f396701fbf32 draft

v0.1.3 Depends on Biopython 1.67 via Tool Shed package or bioconda.

5:b38bbcbd458d

This script is copyright 2010-2013 by Peter Cock, The James Hutton Institute
(formerly SCRI), Scotland, UK. All rights reserved.

See accompanying text file for licence details (MIT license).
"""
import os
import sys
import re

if "-v" in sys.argv or "--version" in sys.argv:
    print("Version 0.1.0")
    sys.exit(0)

def sys_exit(msg, err=1):
   sys.stderr.write(msg.rstrip() + "\n")
   sys.exit(err)

try:
    from Bio.SeqIO.QualityIO import FastqGeneralIterator
except ImportError:
    sys_exit("Biopython missing")

msg = """Expect either 3 or 4 arguments, all FASTQ filenames.

If you want two output files, use four arguments:
 - FASTQ variant (e.g. sanger, solexa, illumina or cssanger)

convention (see http://staden.sourceforge.net/manual/pregap4_unix_50.html
for details), and the new Illumina convention where the reads have the
same identifier with the fragment at the start of the description, e.g.

@HWI-ST916:79:D04M5ACXX:1:1101:10000:100326 1:N:0:TGNCCA
@HWI-ST916:79:D04M5ACXX:1:1101:10000:100326 2:N:0:TGNCCA 

Note that this does support multiple forward and reverse reads per template
(which is quite common with Sanger sequencing), e.g. this which is sorted
alphabetically:

WTSI_1055_4p17.p1kapIBF paired with WTSI_1055_4p17.q1kapIBR
WTSI_1055_4p17.p1kpIBF paired with WTSI_1055_4p17.q1kpIBR
"""

if len(sys.argv) == 5:
    format, input_fastq, pairs_fastq, singles_fastq = sys.argv[1:]
elif len(sys.argv) == 6:
    pairs_fastq = None
    format, input_fastq, pairs_f_fastq, pairs_r_fastq, singles_fastq = sys.argv[1:]
else:
    sys_exit(msg)

format = format.replace("fastq", "").lower()
if not format:
    format="sanger" #safe default
elif format not in ["sanger","solexa","illumina","cssanger"]:
    sys_exit("Unrecognised format %s" % format)

#Cope with three widely used suffix naming convensions,
#Illumina: /1 or /2
#Forward/revered: .f or .r
#Sanger, e.g. .p1k and .q1k
#See http://staden.sourceforge.net/manual/pregap4_unix_50.html
re_f = re.compile(r"(/1|\.f|\.[sfp]\d\w*)$")
re_r = re.compile(r"(/2|\.r|\.[rq]\d\w*)$")

#assert re_f.match("demo/1")
assert re_f.search("demo.f")
assert re_f.search("demo.s1")
assert re_f.search("demo.f1k")
assert re_f.search("demo.p1")
assert re_f.search("demo.p1k")

singles_handle = open(singles_fastq, "w")
last_template, buffered_reads = None, []

for title, seq, qual in FastqGeneralIterator(in_handle):
    count += 1
    name = title.split(None,1)[0]
    is_forward = False
    suffix = re_f.search(name)
    if suffix:
        # ============
        # Forward read

if last_template:
    # Left over singles...
    for old in buffered_reads:
        singles_handle.write(FASTQ_TEMPLATE % old)
        singles += 1
in_handle.close
singles_handle.close()
if pairs_fastq:
    pairs_f_handle.close()
    assert pairs_r_handle.closed
else:

    print("%i reads (%i forward, %i reverse), %i in pairs, %i as singles"
          % (count, forward, reverse, pairs, singles))

assert count == pairs + singles == forward + reverse + neither, \
    "%i vs %i+%i=%i vs %i+%i+%i=%i" \
    % (count,pairs,singles,pairs+singles,forward,reverse,neither,forward+reverse+neither)

6:f396701fbf32

This script is copyright 2010-2013 by Peter Cock, The James Hutton Institute
(formerly SCRI), Scotland, UK. All rights reserved.

See accompanying text file for licence details (MIT license).
"""

import re
import sys

if "-v" in sys.argv or "--version" in sys.argv:
    print("Version 0.1.3")
    sys.exit(0)

try:
    from Bio.SeqIO.QualityIO import FastqGeneralIterator
except ImportError:
    sys.exit("Biopython missing")

msg = """Expect either 3 or 4 arguments, all FASTQ filenames.

If you want two output files, use four arguments:
 - FASTQ variant (e.g. sanger, solexa, illumina or cssanger)

convention (see http://staden.sourceforge.net/manual/pregap4_unix_50.html
for details), and the new Illumina convention where the reads have the
same identifier with the fragment at the start of the description, e.g.

@HWI-ST916:79:D04M5ACXX:1:1101:10000:100326 1:N:0:TGNCCA
@HWI-ST916:79:D04M5ACXX:1:1101:10000:100326 2:N:0:TGNCCA

Note that this does support multiple forward and reverse reads per template
(which is quite common with Sanger sequencing), e.g. this which is sorted
alphabetically:

WTSI_1055_4p17.p1kapIBF paired with WTSI_1055_4p17.q1kapIBR
WTSI_1055_4p17.p1kpIBF paired with WTSI_1055_4p17.q1kpIBR
"""

if len(sys.argv) == 5:
    seq_format, input_fastq, pairs_fastq, singles_fastq = sys.argv[1:]
elif len(sys.argv) == 6:
    pairs_fastq = None
    seq_format, input_fastq, pairs_f_fastq, pairs_r_fastq, singles_fastq = sys.argv[1:]
else:
    sys.exit(msg)

seq_format = seq_format.replace("fastq", "").lower()
if not seq_format:
    seq_format = "sanger"  # safe default
elif seq_format not in ["sanger", "solexa", "illumina", "cssanger"]:
    sys.exit("Unrecognised format %s" % seq_format)

# Cope with three widely used suffix naming convensions,
# Illumina: /1 or /2
# Forward/revered: .f or .r
# Sanger, e.g. .p1k and .q1k
# See http://staden.sourceforge.net/manual/pregap4_unix_50.html
re_f = re.compile(r"(/1|\.f|\.[sfp]\d\w*)$")
re_r = re.compile(r"(/2|\.r|\.[rq]\d\w*)$")

# assert re_f.match("demo/1")
assert re_f.search("demo.f")
assert re_f.search("demo.s1")
assert re_f.search("demo.f1k")
assert re_f.search("demo.p1")
assert re_f.search("demo.p1k")

singles_handle = open(singles_fastq, "w")
last_template, buffered_reads = None, []

for title, seq, qual in FastqGeneralIterator(in_handle):
    count += 1
    name = title.split(None, 1)[0]
    is_forward = False
    suffix = re_f.search(name)
    if suffix:
        # ============
        # Forward read

if last_template:
    # Left over singles...
    for old in buffered_reads:
        singles_handle.write(FASTQ_TEMPLATE % old)
        singles += 1
in_handle.close()
singles_handle.close()
if pairs_fastq:
    pairs_f_handle.close()
    assert pairs_r_handle.closed
else:

    print("%i reads (%i forward, %i reverse), %i in pairs, %i as singles"
          % (count, forward, reverse, pairs, singles))

assert count == pairs + singles == forward + reverse + neither, \
    "%i vs %i+%i=%i vs %i+%i+%i=%i" \
    % (count, pairs, singles, pairs + singles, forward, reverse, neither, forward + reverse + neither)