ó
X¦çIc           @   s<  d  Z  d Z d d l m Z d d l m Z m Z d d l m Z d d l	 m
 Z
 d d l m Z d Z d	 Z d
 „  Z d „  Z d „  Z d „  Z d „  Z e d d „ Z e d d „ Z e d d „ Z d e
 f d „  ƒ  YZ d e
 f d „  ƒ  YZ d e
 f d „  ƒ  YZ e d d „ Z d „  Z e d k r8e ƒ  n  d S(   s½  Bio.SeqIO support for the "fastq" and "qual" file formats.

Note that you are expected to use this code via the Bio.SeqIO interface, as
shown below.

The FASTQ file format is used frequently at the Wellcome Trust Sanger Institute
to bundle a FASTA sequence and its PHRED quality data (integers between 0 and
90).  Rather than using a single FASTQ file, often paired FASTA and QUAL files
are used containing the sequence and the quality information separately.

The PHRED software reads DNA sequencing trace files, calls bases, and
assigns a quality value between 0 and 90 to each called base using a logged
transformation of the error probability, Q = -10 log10( Pe ), for example::

    Pe = 0.0,         Q =  0
    Pe = 0.1,         Q = 10
    Pe = 0.01,        Q = 20
    ...
    Pe = 0.00000001,  Q = 80
    Pe = 0.000000001, Q = 90

In the QUAL format these quality values are held as space separated text in
a FASTA like file format.  In the FASTQ format, each quality values is encoded
with a single ASCI character using chr(Q+33), meaning zero maps to the
character "!" and for example 80 maps to "q".  The sequences and quality are
then stored in pairs in a FASTA like format.

Unfortunately there is no official document describing the FASTQ file format,
and worse, several related but different variants exist.  Reasonable
documentation exists at: http://maq.sourceforge.net/fastq.shtml

Solexa/Illumina quality scores use Q = - 10 log10 ( Pe / (1-Pe) ), which can
be negative or easily exceed 90.  PHRED scores and Solexa scores are NOT
interchangeable (but a reasonable mapping can be achieved between them).
Confusingly Solexa produces a FASTQ like file but using their own score
mapping instead.

Also note that Roche 454 sequencers can output files in the QUAL format, and
thankfully they use PHREP style scores like Sanger.  To extract QUAL files from
a Roche 454 SFF binary file, use the Roche off instrument command line tool
"sffinfo" with the -q or -qual argument.  You can extract a matching FASTA file
using the -s or -seq argument instead.

You are expected to use this module via the Bio.SeqIO functions, with the
following format names:
 - "fastq" means Sanger style FASTQ files using PHRED scores.
 - "fastq-solexa" means Solexa/Illumina style FASTQ files.
 - "qual" means simple quality files using PHRED scores.

For example, consider the following short FASTQ file (extracted from a real
NCBI dataset)::

    @EAS54_6_R1_2_1_413_324
    CCCTTCTTGTCTTCAGCGTTTCTCC
    +
    ;;3;;;;;;;;;;;;7;;;;;;;88
    @EAS54_6_R1_2_1_540_792
    TTGGCAGGCCAAGGCCGATGGATCA
    +
    ;;;;;;;;;;;7;;;;;-;;;3;83
    @EAS54_6_R1_2_1_443_348
    GTTGCTTCTGGCGTGGGTGGGGGGG
    +
    ;;;;;;;;;;;9;7;;.7;393333

This contains three reads of length 25.  From the read length these were
probably originally from an early Solexa/Illumina sequencer but NCBI have
followed the Sanger FASTQ convention and this actually uses PHRED style
qualities.  This means we can parse this file using Bio.SeqIO using "fastq"
as the format name:

    >>> from Bio import SeqIO
    >>> for record in SeqIO.parse(open("Quality/example.fastq"), "fastq") :
    ...     print record.id, record.seq
    EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC
    EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA
    EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG

The qualities are held as a list of integers in each record's annotation:

    >>> print record
    ID: EAS54_6_R1_2_1_443_348
    Name: EAS54_6_R1_2_1_443_348
    Description: EAS54_6_R1_2_1_443_348
    Number of features: 0
    Per letter annotation for: phred_quality
    Seq('GTTGCTTCTGGCGTGGGTGGGGGGG', SingleLetterAlphabet())
    >>> print record.letter_annotations["phred_quality"]
    [26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]

You can use the SeqRecord format method you can show this in the QUAL format:

    >>> print record.format("qual")
    >EAS54_6_R1_2_1_443_348
    26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18
    24 18 18 18 18
    <BLANKLINE>

Or go back to the FASTQ format,

    >>> print record.format("fastq")
    @EAS54_6_R1_2_1_443_348
    GTTGCTTCTGGCGTGGGTGGGGGGG
    +
    ;;;;;;;;;;;9;7;;.7;393333
    <BLANKLINE>

You can also get Biopython to convert the scores and show a Solexa style
FASTQ file:

    >>> print record.format("fastq-solexa")
    @EAS54_6_R1_2_1_443_348
    GTTGCTTCTGGCGTGGGTGGGGGGG
    +
    ZZZZZZZZZZZXZVZZMVZRXRRRR
    <BLANKLINE>

If you wanted to trim your sequences (perhaps to remove low quality regions,
or to remove a primer sequence), try slicing the SeqRecord objects.  e.g.

    >>> sub_rec = record[5:15]
    >>> print sub_rec
    ID: EAS54_6_R1_2_1_443_348
    Name: EAS54_6_R1_2_1_443_348
    Description: EAS54_6_R1_2_1_443_348
    Number of features: 0
    Per letter annotation for: phred_quality
    Seq('TTCTGGCGTG', SingleLetterAlphabet())
    >>> print sub_rec.letter_annotations["phred_quality"]
    [26, 26, 26, 26, 26, 26, 24, 26, 22, 26]
    >>> print sub_rec.format("fastq")
    @EAS54_6_R1_2_1_443_348
    TTCTGGCGTG
    +
    ;;;;;;9;7;
    <BLANKLINE>
    
If you wanted to, you could read in this FASTQ file, and save it as a QUAL file:

    >>> from Bio import SeqIO
    >>> record_iterator = SeqIO.parse(open("Quality/example.fastq"), "fastq")
    >>> out_handle = open("Quality/temp.qual", "w")
    >>> SeqIO.write(record_iterator, out_handle, "qual")
    3
    >>> out_handle.close()

You can of course read in a QUAL file, such as the one we just created:

    >>> from Bio import SeqIO
    >>> for record in SeqIO.parse(open("Quality/temp.qual"), "qual") :
    ...     print record.id, record.seq
    EAS54_6_R1_2_1_413_324 ?????????????????????????
    EAS54_6_R1_2_1_540_792 ?????????????????????????
    EAS54_6_R1_2_1_443_348 ?????????????????????????

Notice that QUAL files don't have a proper sequence present!  But the quality
information is there:

    >>> print record
    ID: EAS54_6_R1_2_1_443_348
    Name: EAS54_6_R1_2_1_443_348
    Description: EAS54_6_R1_2_1_443_348
    Number of features: 0
    Per letter annotation for: phred_quality
    UnknownSeq(25, alphabet = SingleLetterAlphabet(), character = '?')
    >>> print record.letter_annotations["phred_quality"]
    [26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]

Just to keep things tidy, if you are following this example yourself, you can
delete this temporary file now:

    >>> import os
    >>> os.remove("Quality/temp.qual")

Sometimes you won't have a FASTQ file, but rather just a pair of FASTA and QUAL
files.  Because the Bio.SeqIO system is designed for reading single files, you
would have to read the two in separately and then combine the data.  However,
since this is such a common thing to want to do, there is a helper iterator
defined in this module that does this for you - PairedFastaQualIterator.

Alternatively, if you have enough RAM to hold all the records in memory at once,
then a simple dictionary approach would work:

    >>> from Bio import SeqIO
    >>> reads = SeqIO.to_dict(SeqIO.parse(open("Quality/example.fasta"), "fasta"))
    >>> for rec in SeqIO.parse(open("Quality/example.qual"), "qual") :
    ...     reads[rec.id].letter_annotations["phred_quality"]=rec.letter_annotations["phred_quality"]

You can then access any record by its key, and get both the sequence and the
quality scores.

    >>> print reads["EAS54_6_R1_2_1_540_792"].format("fastq")
    @EAS54_6_R1_2_1_540_792
    TTGGCAGGCCAAGGCCGATGGATCA
    +
    ;;;;;;;;;;;7;;;;;-;;;3;83
    <BLANKLINE>

It is important that you explicitly tell Bio.SeqIO which FASTQ variant you are
using ("fastq" for the Sanger standard using PHRED values, or "fastq-solexa"
for the Solexa/Illumina variant), as this cannot be detected reliably
automatically.
s
   epytext eniÿÿÿÿ(   t   single_letter_alphabet(   t   Seqt
   UnknownSeq(   t	   SeqRecord(   t   SequentialSequenceWriter(   t   logi!   i@   c         C   s   d t  d |  d d d ƒ S(   s0  Covert a PHRED quality (range 0 to about 90) to a Solexa quality.

    This will return a floating point number, it is up to you to round this to
    the nearest integer if appropriate.  e.g.

    >>> print "%0.2f" % round(solexa_quality_from_phred(80),2)
    80.00
    >>> print "%0.2f" % round(solexa_quality_from_phred(50),2)
    50.00
    >>> print "%0.2f" % round(solexa_quality_from_phred(20),2)
    19.96
    >>> print "%0.2f" % round(solexa_quality_from_phred(10),2)
    9.54
    >>> print "%0.2f" % round(solexa_quality_from_phred(1),2)
    -5.87
    i
   g      $@i   (   R   (   t   phred_quality(    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   solexa_quality_from_phredã   s    c         C   s   d t  d |  d d d ƒ S(   s3  Convert a Solexa quality (which can be negative) to a PHRED quality.

    This will return a floating point number, it is up to you to round this to
    the nearest integer if appropriate.  e.g.

    >>> print "%0.2f" % round(phred_quality_from_solexa(80),2)
    80.00
    >>> print "%0.2f" % round(phred_quality_from_solexa(20),2)
    20.04
    >>> print "%0.2f" % round(phred_quality_from_solexa(10),2)
    10.41
    >>> print "%0.2f" % round(phred_quality_from_solexa(0),2)
    3.01
    >>> print "%0.2f" % round(phred_quality_from_solexa(-10),2)
    0.41
    i
   g      $@i   (   R   (   t   solexa_quality(    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   phred_quality_from_solexaö   s    c         C   sv   y |  j  d SWn t k
 r" n Xy( g  |  j  d D] } t | ƒ ^ q4 SWn$ t k
 rq t d |  j ƒ ‚ n Xd S(   sÀ   Extract PHRED qualities from a SeqRecord's letter_annotations (PRIVATE).

    If there are no PHRED qualities, but there are Solexa qualities, those are
    used instead after conversion.
    R   R   sL   No suitable quality scores found in letter_annotations of SeqRecord (id=%s).N(   t   letter_annotationst   KeyErrorR	   t
   ValueErrort   id(   t   recordt   q(    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   _get_phred_quality	  s    %c         C   sv   y |  j  d SWn t k
 r" n Xy( g  |  j  d D] } t | ƒ ^ q4 SWn$ t k
 rq t d |  j ƒ ‚ n Xd S(   sÁ   Extract Solexa qualities from a SeqRecord's letter_annotations (PRIVATE).

    If there are no Solexa qualities, but there are PHRED qualities, those are
    used instead after conversion.
    R   R   sK   No suitable quality scores found in letter_annotation of SeqRecord (id=%s).N(   R
   R   R   R   R   (   R   R   (    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   _get_solexa_quality  s    %c         c   s'  x: t  r< |  j ƒ  } | d k r% d S| d d k r Pq q WxÑt  r| d d k re t d ƒ ‚ n  | d j ƒ  } g  } |  j ƒ  } x‡ t  r| s¥ t d ƒ ‚ n  | d d k rî | d j ƒ  rê | d j ƒ  | k rê t d	 ƒ ‚ n  Pn  | j | j ƒ  ƒ |  j ƒ  } qŠ Wd j | ƒ } ~ g  } |  j ƒ  } xk t  r¢| sHPn  | d d k r€t d j | ƒ ƒ t | ƒ k r€Pq€n  | j | j ƒ  ƒ |  j ƒ  } q8Wd j | ƒ } ~ t | ƒ t | ƒ k rõt d
 | t | ƒ t | ƒ f ƒ ‚ n  | | | f V| s@ d Sq@ Wt s#t	 d ƒ ‚ d S(   sF  Iterate over Fastq records as string tuples (not as SeqRecord objects).

    This code does not try to interpret the quality string numerically.  It
    just returns tuples of the title, sequence and quality as strings.  For
    the sequence and quality, any whitespace (such as new lines) is removed.

    Our SeqRecord based FASTQ iterators call this function internally, and then
    turn the strings into a SeqRecord objects, mapping the quality string into
    a list of numerical scores.  If you want to do a custom quality mapping,
    then you might consider calling this function directly.

    For parsing FASTQ files, the title string from the "@" line at the start
    of each record can optionally be omitted on the "+" lines.  If it is
    repeated, it must be identical.

    The sequence string and the quality string can optionally be split over
    multiple lines, although several sources discourage this.  In comparison,
    for the FASTA file format line breaks between 60 and 80 characters are
    the norm.

    WARNING - Because the "@" character can appear in the quality string,
    this can cause problems as this is also the marker for the start of
    a new sequence.  In fact, the "+" sign can also appear as well.  Some
    sources recommended having no line breaks in the  quality to avoid this,
    but even that is not enough, consider this example::

        @071113_EAS56_0053:1:1:998:236
        TTTCTTGCCCCCATAGACTGAGACCTTCCCTAAATA
        +071113_EAS56_0053:1:1:998:236
        IIIIIIIIIIIIIIIIIIIIIIIIIIIIICII+III
        @071113_EAS56_0053:1:1:182:712
        ACCCAGCTAATTTTTGTATTTTTGTTAGAGACAGTG
        +
        @IIIIIIIIIIIIIIICDIIIII<%<6&-*).(*%+
        @071113_EAS56_0053:1:1:153:10
        TGTTCTGAAGGAAGGTGTGCGTGCGTGTGTGTGTGT
        +
        IIIIIIIIIIIICIIGIIIII>IAIIIE65I=II:6
        @071113_EAS56_0053:1:3:990:501
        TGGGAGGTTTTATGTGGA
        AAGCAGCAATGTACAAGA
        +
        IIIIIII.IIIIII1@44
        @-7.%<&+/$/%4(++(%

    This is four PHRED encoded FASTQ entries originally from an NCBI source
    (given the read length of 36, these are probably Solexa Illumna reads where
    the quality has been mapped onto the PHRED values).

    This example has been edited to illustrate some of the nasty things allowed
    in the FASTQ format.  Firstly, on the "+" lines most but not all of the
    (redundant) identifiers are ommited.  In real files it is likely that all or
    none of these extra identifiers will be present.

    Secondly, while the first three sequences have been shown without line
    breaks, the last has been split over multiple lines.  In real files any line
    breaks are likely to be consistent.

    Thirdly, some of the quality string lines start with an "@" character.  For
    the second record this is unavoidable.  However for the fourth sequence this
    only happens because its quality string is split over two lines.  A naive
    parser could wrongly treat any line starting with an "@" as the beginning of
    a new sequence!  This code copes with this possible ambiguity by keeping track
    of the length of the sequence which gives the expected length of the quality
    string.

    Using this tricky example file as input, this short bit of code demonstrates
    what this parsing function would return:

    >>> handle = open("Quality/tricky.fastq", "rU")
    >>> for (title, sequence, quality) in FastqGeneralIterator(handle) :
    ...     print title
    ...     print sequence, quality
    071113_EAS56_0053:1:1:998:236
    TTTCTTGCCCCCATAGACTGAGACCTTCCCTAAATA IIIIIIIIIIIIIIIIIIIIIIIIIIIIICII+III
    071113_EAS56_0053:1:1:182:712
    ACCCAGCTAATTTTTGTATTTTTGTTAGAGACAGTG @IIIIIIIIIIIIIIICDIIIII<%<6&-*).(*%+
    071113_EAS56_0053:1:1:153:10
    TGTTCTGAAGGAAGGTGTGCGTGCGTGTGTGTGTGT IIIIIIIIIIIICIIGIIIII>IAIIIE65I=II:6
    071113_EAS56_0053:1:3:990:501
    TGGGAGGTTTTATGTGGAAAGCAGCAATGTACAAGA IIIIIII.IIIIII1@44@-7.%<&+/$/%4(++(%
    >>> handle.close()

    Finally we note that some sources state that the quality string should
    start with "!" (which using the PHRED mapping means the first letter always
    has a quality score of zero).  This rather restrictive rule is not widely
    observed, so is therefore ignored here.  One plus point about this "!" rule
    is that (provided there are no line breaks in the quality sequence) it
    would prevent the above problem with the "@" character.
    t    Ni    t   @s6   Records in Fastq files should start with '@' characteri   s(   End of file without quality information.t   +s%   Sequence and quality captions differ.sC   Lengths of sequence and quality values differs  for %s (%i and %i).s   Should not reach this line(
   t   Truet   readlineR   t   rstript   extendt   splitt   joint   lent   Falset   AssertionError(   t   handlet   linet
   title_linet	   seq_linest
   seq_stringt   quality_linest   quality_string(    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   FastqGeneralIterator-  sT    \	 		&	 !" c      	   c   s  xt  |  ƒ D]ù \ } } } | r: | | ƒ \ } } } n | } | j ƒ  d } | } t t | | ƒ d | d | d | ƒ}	 t t d ƒ k s• t ‚ g  | D] }
 t |
 ƒ t ^ qœ } | rô t | ƒ d k  sâ t | ƒ d k rô t	 d ƒ ‚ qô n  | |	 j
 d <|	 Vq Wd	 S(
   sú	  Generator function to iterate over FASTQ records (as SeqRecord objects).

     - handle - input file
     - alphabet - optional alphabet
     - title2ids - A function that, when given the title line from the FASTQ
                   file (without the beginning >), will return the id, name and
                   description (in that order) for the record as a tuple of
                   strings.  If this is not given, then the entire title line
                   will be used as the description, and the first word as the
                   id and name.

    Note that use of title2ids matches that of Bio.SeqIO.FastaIO.

    For each sequence in a (Sanger style) FASTQ file there is a matching string
    encoding the PHRED qualities (integers between 0 and about 90) using ASCII
    values with an offset of 33.

    For example, consider a file containing three short reads::

        @EAS54_6_R1_2_1_413_324
        CCCTTCTTGTCTTCAGCGTTTCTCC
        +
        ;;3;;;;;;;;;;;;7;;;;;;;88
        @EAS54_6_R1_2_1_540_792
        TTGGCAGGCCAAGGCCGATGGATCA
        +
        ;;;;;;;;;;;7;;;;;-;;;3;83
        @EAS54_6_R1_2_1_443_348
        GTTGCTTCTGGCGTGGGTGGGGGGG
        +
        ;;;;;;;;;;;9;7;;.7;393333

    For each sequence (e.g. "CCCTTCTTGTCTTCAGCGTTTCTCC") there is a matching
    string encoding the PHRED qualities using a ASCI values with an offset of
    33 (e.g. ";;3;;;;;;;;;;;;7;;;;;;;88").

    Using this module directly you might run:

    >>> handle = open("Quality/example.fastq", "rU")
    >>> for record in FastqPhredIterator(handle) :
    ...     print record.id, record.seq
    EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC
    EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA
    EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG
    >>> handle.close()

    Typically however, you would call this via Bio.SeqIO instead with "fastq" as
    the format:

    >>> from Bio import SeqIO
    >>> handle = open("Quality/example.fastq", "rU")
    >>> for record in SeqIO.parse(handle, "fastq") :
    ...     print record.id, record.seq
    EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC
    EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA
    EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG
    >>> handle.close()

    If you want to look at the qualities, they are record in each record's
    per-letter-annotation dictionary as a simple list of integers:

    >>> print record.letter_annotations["phred_quality"]
    [26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]
    i    R   t   namet   descriptiont   !iZ   s‰   Quality score outside 0 to 90 found - these are perhaps in a Solexa/Illumina format, not the Sanger FASTQ format which uses PHRED scores.R   N(   R%   R   R   R   t   SANGER_SCORE_OFFSETt   ordR   t   mint   maxR   R
   (   R   t   alphabett	   title2idsR    R"   R$   R   R&   t   descrR   t   lettert	   qualities(    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   FastqPhredIteratorÄ  s    A
#$c      	   c   s´   x­ t  |  ƒ D]Ÿ \ } } } | r4 | \ } } } n | } | j ƒ  d } | } t t | | ƒ d | d | d | ƒ}	 g  | D] }
 t |
 ƒ t ^ q~ } | |	 j d <|	 Vq Wd S(   så  Parsing the Solexa/Illumina FASTQ like files (which differ in the quality mapping).

    The optional arguments are the same as those for the FastqPhredIterator.

    For each sequence in Solexa/Illumina FASTQ files there is a matching string
    encoding the Solexa integer qualities using ASCII values with an offset
    of 64.  Solexa scores are scaled differently to PHRED scores, and Biopython
    will NOT perform any automatic conversion when loading.

    For example, consider a file containing these five records::
        
        @SLXA-B3_649_FC8437_R1_1_1_610_79
        GATGTGCAATACCTTTGTAGAGGAA
        +SLXA-B3_649_FC8437_R1_1_1_610_79
        YYYYYYYYYYYYYYYYYYWYWYYSU
        @SLXA-B3_649_FC8437_R1_1_1_397_389
        GGTTTGAGAAAGAGAAATGAGATAA
        +SLXA-B3_649_FC8437_R1_1_1_397_389
        YYYYYYYYYWYYYYWWYYYWYWYWW
        @SLXA-B3_649_FC8437_R1_1_1_850_123
        GAGGGTGTTGATCATGATGATGGCG
        +SLXA-B3_649_FC8437_R1_1_1_850_123
        YYYYYYYYYYYYYWYYWYYSYYYSY
        @SLXA-B3_649_FC8437_R1_1_1_362_549
        GGAAACAAAGTTTTTCTCAACATAG
        +SLXA-B3_649_FC8437_R1_1_1_362_549
        YYYYYYYYYYYYYYYYYYWWWWYWY
        @SLXA-B3_649_FC8437_R1_1_1_183_714
        GTATTATTTAATGGCATACACTCAA
        +SLXA-B3_649_FC8437_R1_1_1_183_714
        YYYYYYYYYYWYYYYWYWWUWWWQQ
        
    Using this module directly you might run:

    >>> handle = open("Quality/solexa_example.fastq", "rU")
    >>> for record in FastqSolexaIterator(handle) :
    ...     print record.id, record.seq
    SLXA-B3_649_FC8437_R1_1_1_610_79 GATGTGCAATACCTTTGTAGAGGAA
    SLXA-B3_649_FC8437_R1_1_1_397_389 GGTTTGAGAAAGAGAAATGAGATAA
    SLXA-B3_649_FC8437_R1_1_1_850_123 GAGGGTGTTGATCATGATGATGGCG
    SLXA-B3_649_FC8437_R1_1_1_362_549 GGAAACAAAGTTTTTCTCAACATAG
    SLXA-B3_649_FC8437_R1_1_1_183_714 GTATTATTTAATGGCATACACTCAA
    >>> handle.close()

    Typically however, you would call this via Bio.SeqIO instead with "fastq" as
    the format:

    >>> from Bio import SeqIO
    >>> handle = open("Quality/solexa_example.fastq", "rU")
    >>> for record in SeqIO.parse(handle, "fastq-solexa") :
    ...     print record.id, record.seq
    SLXA-B3_649_FC8437_R1_1_1_610_79 GATGTGCAATACCTTTGTAGAGGAA
    SLXA-B3_649_FC8437_R1_1_1_397_389 GGTTTGAGAAAGAGAAATGAGATAA
    SLXA-B3_649_FC8437_R1_1_1_850_123 GAGGGTGTTGATCATGATGATGGCG
    SLXA-B3_649_FC8437_R1_1_1_362_549 GGAAACAAAGTTTTTCTCAACATAG
    SLXA-B3_649_FC8437_R1_1_1_183_714 GTATTATTTAATGGCATACACTCAA
    >>> handle.close()

    If you want to look at the qualities, they are recorded in each record's
    per-letter-annotation dictionary as a simple list of integers:

    >>> print record.letter_annotations["solexa_quality"]
    [25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 23, 25, 25, 25, 25, 23, 25, 23, 23, 21, 23, 23, 23, 17, 17]

    These scores aren't very good, but they are high enough that they map
    almost exactly onto PHRED scores:

    >>> print "%0.2f" % phred_quality_from_solexa(25)
    25.01

    Let's look at another example read which is even worse, where there are
    more noticeable differences between the Solexa and PHRED scores::

         @slxa_0013_1_0001_24
         ACAAAAATCACAAGCATTCTTATACACC
         +slxa_0013_1_0001_24
         ??????????????????:??<?<-6%.

    Again, you would typically use Bio.SeqIO to read this file in (rather than
    calling the Bio.SeqIO.QualtityIO module directly).  Most FASTQ files will
    contain thousands of reads, so you would normally use Bio.SeqIO.parse()
    as shown above.  This example has only as one entry, so instead we can
    use the Bio.SeqIO.read() function:

    >>> from Bio import SeqIO
    >>> handle = open("Quality/solexa.fastq", "rU")
    >>> record = SeqIO.read(handle, "fastq-solexa")
    >>> handle.close()
    >>> print record.id, record.seq
    slxa_0013_1_0001_24 ACAAAAATCACAAGCATTCTTATACACC
    >>> print record.letter_annotations["solexa_quality"]
    [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -6, -1, -1, -4, -1, -4, -19, -10, -27, -18]

    These quality scores are so low that when converted from the Solexa scheme
    into PHRED scores they look quite different:

    >>> print "%0.2f" % phred_quality_from_solexa(-1)
    2.54

    Note you can use the Bio.SeqIO.write() function or the SeqRecord's format
    method to output the record(s):

    >>> print record.format("fastq-solexa")
    @slxa_0013_1_0001_24
    ACAAAAATCACAAGCATTCTTATACACC
    +
    ??????????????????:??<?<-6%.
    <BLANKLINE>

    Note this output is slightly different from the input file as Biopython
    has left out the optional repetition of the sequence identifier on the "+"
    line.  If you want the to use PHRED scores, use "fastq" or "qual" as the
    output format instead, and Biopython will do the conversion for you:

    >>> print record.format("fastq")
    @slxa_0013_1_0001_24
    ACAAAAATCACAAGCATTCTTATACACC
    +
    $$$$$$$$$$$$$$$$$$"$$"$"!!!!
    <BLANKLINE>

    >>> print record.format("qual")
    >slxa_0013_1_0001_24
    3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 3 3 1 3 1 0 0 0 0
    <BLANKLINE>
    i    R   R&   R'   R   N(   R%   R   R   R   R*   t   SOLEXA_SCORE_OFFSETR
   (   R   R-   R.   R    R"   R$   R   R&   R/   R   R0   R1   (    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   FastqSolexaIterator#  s    #c   
      c   sè  x: t  r< |  j ƒ  } | d k r% d S| d d k r Pq q Wx’t  rÑ| d d k re t d ƒ ‚ n  | r | | d j ƒ  ƒ \ } } } n& | d j ƒ  } | j ƒ  d } | } g  } |  j ƒ  } x` t  r'| sØ Pn  | d d k rì Pn  | j g  | j ƒ  D] } t | ƒ ^ qÿ ƒ |  j ƒ  } qÈ W| r…t | ƒ d k  sRt | ƒ d k r…t d d	 d
 | t | ƒ t | ƒ f ƒ ‚ q…n  t	 t
 t | ƒ | ƒ d | d | d | ƒ}	 | |	 j d <|	 V| s@ d Sq@ Wt sät d ƒ ‚ d S(   s
  For QUAL files which include PHRED quality scores, but no sequence.

    For example, consider this short QUAL file::

        >EAS54_6_R1_2_1_413_324
        26 26 18 26 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26
        26 26 26 23 23
        >EAS54_6_R1_2_1_540_792
        26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 12 26 26
        26 18 26 23 18
        >EAS54_6_R1_2_1_443_348
        26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18
        24 18 18 18 18
    
    Using this module directly you might run:

    >>> handle = open("Quality/example.qual", "rU")
    >>> for record in QualPhredIterator(handle) :
    ...     print record.id, record.seq
    EAS54_6_R1_2_1_413_324 ?????????????????????????
    EAS54_6_R1_2_1_540_792 ?????????????????????????
    EAS54_6_R1_2_1_443_348 ?????????????????????????
    >>> handle.close()

    Typically however, you would call this via Bio.SeqIO instead with "qual"
    as the format:

    >>> from Bio import SeqIO
    >>> handle = open("Quality/example.qual", "rU")
    >>> for record in SeqIO.parse(handle, "qual") :
    ...     print record.id, record.seq
    EAS54_6_R1_2_1_413_324 ?????????????????????????
    EAS54_6_R1_2_1_540_792 ?????????????????????????
    EAS54_6_R1_2_1_443_348 ?????????????????????????
    >>> handle.close()

    Becase QUAL files don't contain the sequence string itself, the seq
    property is set to an UnknownSeq object.  As no alphabet was given, this
    has defaulted to a generic single letter alphabet and the character "?"
    used.

    By specifying a nucleotide alphabet, "N" is used instead:

    >>> from Bio import SeqIO
    >>> from Bio.Alphabet import generic_dna
    >>> handle = open("Quality/example.qual", "rU")
    >>> for record in SeqIO.parse(handle, "qual", alphabet=generic_dna) :
    ...     print record.id, record.seq
    EAS54_6_R1_2_1_413_324 NNNNNNNNNNNNNNNNNNNNNNNNN
    EAS54_6_R1_2_1_540_792 NNNNNNNNNNNNNNNNNNNNNNNNN
    EAS54_6_R1_2_1_443_348 NNNNNNNNNNNNNNNNNNNNNNNNN
    >>> handle.close()

    However, the quality scores themselves are available as a list of integers
    in each record's per-letter-annotation:

    >>> print record.letter_annotations["phred_quality"]
    [26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]

    You can still slice one of these SeqRecord objects with an UnknownSeq:

    >>> sub_record = record[5:10]
    >>> print sub_record.id, sub_record.letter_annotations["phred_quality"]
    EAS54_6_R1_2_1_443_348 [26, 26, 26, 26, 26]
    R   Ni    t   >s6   Records in Fasta files should start with '>' characteri   iZ   s4   Quality score range for %s is %i to %i, outside the s5   expected 0 to 90.  Perhaps these are Solexa/Illumina s   scores, and not PHRED scores?R   R&   R'   R   s   Should not reach this line(   R   R   R   R   R   R   t   intR+   R,   R   R   R   R
   R   R   (
   R   R-   R.   R   R   R&   R/   R1   t   wordR   (    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   QualPhredIterator°  sF    C	 	"	  ,$% t   FastqPhredWriterc           B   s   e  Z d  Z d „  Z RS(   së  Class to write FASTQ format files (using PHRED quality scores).

    Although you can use this class directly, you are strongly encouraged
    to use the Bio.SeqIO.write() function instead.  For example, this code
    reads in a FASTQ (PHRED) file and re-saves it as another FASTQ (PHRED)
    file:

    >>> from Bio import SeqIO
    >>> record_iterator = SeqIO.parse(open("Quality/example.fastq"), "fastq")
    >>> out_handle = open("Quality/temp.fastq", "w")
    >>> SeqIO.write(record_iterator, out_handle, "fastq")
    3
    >>> out_handle.close()

    You might want to do this if the original file included extra line breaks,
    which while valid may not be supported by all tools.  The output file from
    Biopython will have each sequence on a single line, and each quality
    string on a single line (which is considered desirable for maximum
    compatibility).

    In this next example, a Solexa FASTQ file is converted into a standard
    Sanger style FASTQ file using PHRED qualities:

    >>> from Bio import SeqIO
    >>> record_iterator = SeqIO.parse(open("Quality/solexa.fastq"), "fastq-solexa")
    >>> out_handle = open("Quality/temp.fastq", "w")
    >>> SeqIO.write(record_iterator, out_handle, "fastq")
    1
    >>> out_handle.close()

    This code is also called if you use the .format("fastq") method of a
    SeqRecord.

    P.S. To avoid cluttering up your working directory, you can delete this
    temporary file now:

    >>> import os
    >>> os.remove("Quality/temp.fastq")

    c         C   s  |  j  s t ‚ |  j s t ‚ t |  _ t t d ƒ k s@ t ‚ d j g  t | ƒ D]% } t	 t
 t | t d ƒ ƒ ƒ ^ qS ƒ } | j d k r¦ t d | j ƒ ‚ n  t | ƒ t | ƒ k ré t d | j t | ƒ t | ƒ f ƒ ‚ n  |  j | j ƒ } |  j j d | | j | f ƒ d S(   s(   Write a single FASTQ record to the file.R(   R   i    s   No sequence for record %ss6   Record %s has sequence length %i but %i quality scoress   @%s
%s
+
%s
N(   t   _header_writtenR   t   _footer_writtenR   t   _record_writtenR)   R*   R   R   t   chrR6   t   roundt   seqt   NoneR   R   R   t   cleanR   t   write(   t   selfR   R   R1   t   title(    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   write_recordD  s    		8%(   t   __name__t
   __module__t   __doc__RE   (    (    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyR9     s   (t   QualPhredWriterc           B   s&   e  Z d  Z d d d „ Z d „  Z RS(   sþ  Class to write QUAL format files (using PHRED quality scores).

    Although you can use this class directly, you are strongly encouraged
    to use the Bio.SeqIO.write() function instead.  For example, this code
    reads in a FASTQ file and saves the quality scores into a QUAL file:

    >>> from Bio import SeqIO
    >>> record_iterator = SeqIO.parse(open("Quality/example.fastq"), "fastq")
    >>> out_handle = open("Quality/temp.qual", "w")
    >>> SeqIO.write(record_iterator, out_handle, "qual")
    3
    >>> out_handle.close()

    This code is also called if you use the .format("qual") method of a
    SeqRecord.

    P.S. Don't forget to clean up the temp file if you don't need it anymore:

    >>> import os
    >>> os.remove("Quality/temp.qual")
    i<   c         C   sM   t  j |  | ƒ d |  _ | r7 | d k  r7 t ‚ q7 n  | |  _ | |  _ d S(   sx  Create a QUAL writer.

        Arguments:
         - handle - Handle to an output file, e.g. as returned
                    by open(filename, "w")
         - wrap   - Optional line length used to wrap sequence lines.
                    Defaults to wrapping the sequence at 60 characters
                    Use zero (or None) for no wrapping, giving a single
                    long line for the sequence.
         - record2title - Optional function to return the text to be
                    used for the title line of each record.  By default
                    a combination of the record.id and record.description
                    is used.  If the record.description starts with the
                    record.id, then just the record.description is used.

        The record2title argument is present for consistency with the
        Bio.SeqIO.FastaIO writer class.
        i   N(   R   t   __init__R@   t   wrapR   t   record2title(   RC   R   RK   RL   (    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyRJ   n  s    		c   	      C   s½  |  j  s t ‚ |  j s t ‚ t |  _ |  j rL |  j |  j | ƒ ƒ } n_ |  j | j ƒ } |  j | j ƒ } | r› | j	 d	 d ƒ d | k r› | } n d | | f } d | k s½ t ‚ d | k sÏ t ‚ |  j j d | ƒ g  t | ƒ D] } d t | d ƒ ^ qð } |  j r–xž | r’| j d ƒ } xH | rzt | ƒ d t | d ƒ |  j k  rz| d | j d ƒ 7} q3W|  j j | d ƒ qWn# d j | ƒ } |  j j | d ƒ d	 S(
   s'   Write a single QUAL record to the file.i   i    s   %s %ss   
s   s   >%s
s   %it    N(   R:   R   R;   R   R<   RL   RA   R   R'   R   R@   R   RB   R   R>   RK   t   popR   R   (	   RC   R   RD   R   R'   R   R1   R   t   data(    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyRE   Š  s.    		"	,			'N(   RF   RG   RH   R@   RJ   RE   (    (    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyRI   X  s   t   FastqSolexaWriterc           B   s   e  Z d  Z d „  Z RS(   sS  Class to write FASTQ format files (using Solexa quality scores).

    Although you can use this class directly, you are strongly encouraged
    to use the Bio.SeqIO.write() function instead.  For example, this code
    reads in a FASTQ file and re-saves it as another FASTQ file:

    >>> from Bio import SeqIO
    >>> record_iterator = SeqIO.parse(open("Quality/solexa.fastq"), "fastq-solexa")
    >>> out_handle = open("Quality/temp.fastq", "w")
    >>> SeqIO.write(record_iterator, out_handle, "fastq-solexa")
    1
    >>> out_handle.close()

    You might want to do this if the original file included extra line
    breaks, which (while valid) may not be supported by all tools.  The
    output file from Biopython will have each sequence on a single line, and
    each quality string on a single line (which is considered desirable for
    maximum compatibility).

    This code is also called if you use the .format("fastq-solexa") method of
    a SeqRecord.

    P.S. Don't forget to delete the temp file if you don't need it anymore:

    >>> import os
    >>> os.remove("Quality/temp.fastq")
    c         C   s  |  j  s t ‚ |  j s t ‚ t |  _ d j g  t | ƒ D]% } t t t	 | t
 d ƒ ƒ ƒ ^ q; ƒ } | j d k rŽ t d | j ƒ ‚ n  t | ƒ t | ƒ k rÑ t d | j t | ƒ t | ƒ f ƒ ‚ n  |  j | j ƒ } |  j j d | | j | f ƒ d S(   s(   Write a single FASTQ record to the file.R   i    s   No sequence for record %ss6   Record %s has sequence length %i but %i quality scoress   @%s
%s
+
%s
N(   R:   R   R;   R   R<   R   R   R=   R6   R>   R3   R?   R@   R   R   R   RA   R   RB   (   RC   R   R   R1   RD   (    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyRE   Ì  s    		8%(   RF   RG   RH   RE   (    (    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyRP   °  s   c   	      c   sv  d d l  m } | |  d | d | ƒ} t | d | d | ƒ} x/t rqy | j ƒ  } Wn t k
 rr d
 } n Xy | j ƒ  } Wn t k
 rœ d
 } n X| d
 k r¹ | d
 k r¹ Pn  | d
 k rÔ t d ƒ ‚ n  | d
 k rï t d ƒ ‚ n  | j | j k r t d | j | j f ƒ ‚ n  t	 | ƒ t	 | j
 d ƒ k rUt d	 | j ƒ ‚ n  | j
 d | j
 d <| VqC Wd
 S(   sé	  Iterate over matched FASTA and QUAL files as SeqRecord objects.

    For example, consider this short QUAL file::

        >EAS54_6_R1_2_1_413_324
        26 26 18 26 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26
        26 26 26 23 23
        >EAS54_6_R1_2_1_540_792
        26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 12 26 26
        26 18 26 23 18
        >EAS54_6_R1_2_1_443_348
        26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18
        24 18 18 18 18
    
    And a matching FASTA file::

        >EAS54_6_R1_2_1_413_324
        CCCTTCTTGTCTTCAGCGTTTCTCC
        >EAS54_6_R1_2_1_540_792
        TTGGCAGGCCAAGGCCGATGGATCA
        >EAS54_6_R1_2_1_443_348
        GTTGCTTCTGGCGTGGGTGGGGGGG

    You can parse these separately using Bio.SeqIO with the "qual" and
    "fasta" formats, but then you'll get a group of SeqRecord objects with
    no sequence, and a matching group with the sequence but not the
    qualities.  Because it only deals with one input file handle, Bio.SeqIO
    can't be used to read the two files together - but this function can!
    For example,
    
    >>> rec_iter = PairedFastaQualIterator(open("Quality/example.fasta", "rU"),
    ...                                    open("Quality/example.qual", "rU"))
    >>> for record in rec_iter :
    ...     print record.id, record.seq
    EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC
    EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA
    EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG

    As with the FASTQ or QUAL parsers, if you want to look at the qualities,
    they are in each record's per-letter-annotation dictionary as a simple
    list of integers:

    >>> print record.letter_annotations["phred_quality"]
    [26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]

    If you have access to data as a FASTQ format file, using that directly
    would be simpler and more straight forward.  Note that you can easily use
    this function to convert paired FASTA and QUAL files into FASTQ files:

    >>> from Bio import SeqIO
    >>> rec_iter = PairedFastaQualIterator(open("Quality/example.fasta", "rU"),
    ...                                    open("Quality/example.qual", "rU"))
    >>> out_handle = open("Quality/temp.fastq", "w")
    >>> SeqIO.write(rec_iter, out_handle, "fastq")
    3
    >>> out_handle.close()

    And don't forget to clean up the temp file if you don't need it anymore:

    >>> import os
    >>> os.remove("Quality/temp.fastq")    
    iÿÿÿÿ(   t   FastaIteratorR-   R.   s/   FASTA file has more entries than the QUAL file.s/   QUAL file has more entries than the FASTA file.s/   FASTA and QUAL entries do not match (%s vs %s).R   s<   Sequence length and number of quality scores disagree for %sN(   t   Bio.SeqIO.FastaIORQ   R8   R   t   nextt   StopIterationR@   R   R   R   R
   (	   t   fasta_handlet   qual_handleR-   R.   RQ   t
   fasta_itert	   qual_itert   f_rect   q_rec(    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   PairedFastaQualIteratorÞ  s8    ?			

c          C   s  d d l  }  d d l } | j j | j j d d d ƒ ƒ rd GH| j j | j ƒ } | j | j j d d d ƒ ƒ | j j d ƒ s t	 ‚ | j j d ƒ s¥ t	 ‚ | j j d ƒ s½ t	 ‚ | j j d	 ƒ sÕ t	 ‚ | j j d
 ƒ sí t	 ‚ |  j
 ƒ  | j | ƒ ~ d GHn  d S(   sÁ   Run the Bio.SeqIO module's doctests.

    This will try and locate the unit tests directory, and run the doctests
    from there in order that the relative paths used in the examples work.
    iÿÿÿÿNs   ..t   Testss   Runing doctests...s   Quality/example.fastqs   Quality/example.fastas   Quality/example.quals   Quality/tricky.fastqs   Quality/solexa.fastqt   Done(   t   doctestt   ost   patht   isdirR   t   abspatht   curdirt   chdirt   isfileR   t   testmod(   R^   R_   t   cur_dir(    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   _testA  s    $
t   __main__N(   RH   t   __docformat__t   Bio.AlphabetR    t   Bio.SeqR   R   t   Bio.SeqRecordR   t
   InterfacesR   t   mathR   R)   R3   R   R	   R   R   R%   R@   R2   R4   R8   R9   RI   RP   R[   Rh   RF   (    (    (    s†   /oak/stanford/groups/akundaje/marinovg/programs/biopython-1.50.tar.gz/biopython-1.50/build/lib.linux-x86_64-2.7/Bio/SeqIO/QualityIO.pyt   <module>Ó   s.   					—_k=X.c	