diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,675 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc.
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+
+ Copyright (C)
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see .
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+ Copyright (C)
+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+.
+
+ The GNU General Public License does not permit incorporating your program
+into proprietary programs. If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License. But first, please read
+.
+
diff --git a/RNAlien.cabal b/RNAlien.cabal
new file mode 100644
--- /dev/null
+++ b/RNAlien.cabal
@@ -0,0 +1,75 @@
+-- Initial RNAlien.cabal generated by cabal init. For further documentation,
+-- see http://haskell.org/cabal/users-guide/
+
+name: RNAlien
+version: 1.0.0
+synopsis: Unsupervized construction of RNA family models
+description: RNAlien is a tool for automatic construction of RNAfamily models from a single sequence.
+ .
+ It is available as a commandline tool, for testing or construction of few sequences the webservice can be used.
+ .
+ The source code of RNAlien, as well as the webserver is open source and available via GitHub (License GPL-3):
+ .
+ *
+ .
+ *
+ .
+ TaxonomyTools which can be used to visualise the organisms included in a RNAlien result can be found here (License GPL-3):
+ .
+ *
+ .
+ *
+ .
+ For instruction how to use RNAlien please see the .
+ .
+ Dependencies:
+ .
+ *
+ .
+ *
+ .
+ *
+ .
+ *
+ .
+ *
+ .
+ Installation via cabal-install:
+ .
+ > cabal install RNAlien
+
+license: GPL-3
+license-file: LICENSE
+author: Florian Eggenhofer
+maintainer: egg@tbi.univie.ac.at
+-- copyright:
+category: Bioinformatics
+build-type: Simple
+cabal-version: >=1.8
+
+executable RNAlien
+ Hs-Source-Dirs: ./src/Bio/
+ main-is: RNAlien.hs
+ ghc-options: -Wall -O2
+ build-depends: base >=4.5 && <5, cmdargs, directory, biofasta, random, biocore, containers, RNAlien, time, either-unwrap, filepath
+
+executable RNAlienStatistics
+ Hs-Source-Dirs: ./src/Bio/
+ main-is: RNAlienStatistics.hs
+ ghc-options: -Wall -O2
+ build-depends: base >=4.5 && <5, cmdargs, cassava, vector, process, bytestring, either-unwrap, RNAlien, directory, biofasta, biocore, split, filepath, ViennaRNAParser
+
+executable SelectSequences
+ Hs-Source-Dirs: ./src/Bio/
+ main-is: SelectSequences.hs
+ ghc-options: -Wall -O2
+ build-depends: base >=4.5 && <5, cmdargs, either-unwrap, RNAlien
+
+Library
+ Hs-Source-Dirs: ./src/
+ ghc-options: -Wall -O2 -fno-warn-unused-do-bind
+ build-depends: base >=4.5 && <5, cmdargs, ViennaRNAParser>=1.2.6, process, directory, blastxml>=0.3.2, biofasta, parsec, random, BlastHTTP, biocore, bytestring, Taxonomy >= 1.0.1, either-unwrap, containers, ClustalParser>=1.1.0, EntrezHTTP>=1.0.0, split, vector, edit-distance, cassava, matrix, hierarchical-clustering, filepath, HTTP
+ Exposed-Modules: Bio.RNAlienData
+ Bio.RNAlienLibrary
+
+
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/src/Bio/RNAlien.hs b/src/Bio/RNAlien.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/RNAlien.hs
@@ -0,0 +1,100 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+
+-- | Unsupervized construction of RNA family models
+-- For more information on RNA family models consult
+-- Testcommand: dist/build/RNAlien/RNAlien -i ~egg/initialfasta/RybB.fa -c 3 -o /scr/kronos/egg/temp/ > ~egg/Desktop/alieninitialtest
+module Main where
+
+import System.Console.CmdArgs
+import System.Directory
+import Bio.Sequence.Fasta
+import Bio.RNAlienData
+import Bio.RNAlienLibrary
+import Data.Maybe
+import Data.Either.Unwrap
+import Data.Time
+import qualified System.FilePath as FP
+
+data Options = Options
+ { inputFastaFilePath :: String,
+ outputPath :: String,
+ inputTaxId :: Maybe Int,
+ inputnSCICutoff :: Maybe Double,
+ inputEvalueCutoff :: Maybe Double,
+ inputBlastDatabase :: Maybe String,
+ lengthFilter :: Bool,
+ coverageFilter :: Bool,
+ singleHitperTax :: Bool,
+ threads :: Int,
+ taxonomyRestriction :: Maybe String,
+ sessionIdentificator :: Maybe String
+ } deriving (Show,Data,Typeable)
+
+options :: Options
+options = Options
+ { inputFastaFilePath = def &= name "i" &= help "Path to input fasta file",
+ outputPath = def &= name "o" &= help "Path to output directory",
+ inputTaxId = Nothing &= name "t" &= help "NCBI taxonomy ID number of input RNA organism",
+ inputnSCICutoff = (Just (1 :: Double)) &= name "z" &= help "Only candidate sequences with a normalized structure conservation index (nSCI) higher than this value are accepted. Default: 1",
+ inputEvalueCutoff = (Just (0.001 :: Double)) &= name "e" &= help "Evalue cutoff for cmsearch filtering. Default: 0.001",
+ inputBlastDatabase = Just "nt" &= name "b" &= help "Specify name of blast database to use. Default: nt",
+ lengthFilter = True &= name "l" &= help "Filter blast hits per genomic length. Default: True",
+ coverageFilter = True &= name "a" &= help "Filter blast hits by coverage of at least 80%. Default: True",
+ singleHitperTax = True &= name "s" &= help "Only the best blast hit per taxonomic entry is considered. Default: True",
+ threads = 1 &= name "c" &= help "Number of available cpu slots/cores. Default: 1",
+ taxonomyRestriction = Nothing &= name "r" &= help "Restrict search space to taxonomic kingdom (bacteria,archea,eukaryia). Default: not set",
+ sessionIdentificator = Nothing &= name "d" &= help "Optional session id that is used instead of automatically generated one."
+ } &= summary "RNAlien version 1.0.0" &= help "Florian Eggenhofer, Ivo L. Hofacker, Christian Höner zu Siederdissen - 2013 - 2015" &= verbosity
+
+main :: IO ()
+main = do
+ Options{..} <- cmdArgs options
+ verboseLevel <- getVerbosity
+ -- Generate SessionID
+ sessionId <- createSessionID sessionIdentificator
+ timestamp <- getCurrentTime
+ let temporaryDirectoryPath = FP.addTrailingPathSeparator outputPath ++ sessionId ++ "/"
+ createDirectoryIfMissing False temporaryDirectoryPath
+ createDirectoryIfMissing False (temporaryDirectoryPath ++ "log")
+ -- Create Log files
+ writeFile (temporaryDirectoryPath ++ "Log") ("RNAlien 1.0.0" ++ "\n")
+ writeFile (temporaryDirectoryPath ++ "log/warnings") ("")
+ logMessage ("Timestamp: " ++ (show timestamp) ++ "\n") temporaryDirectoryPath
+ logMessage ("Temporary Directory: " ++ temporaryDirectoryPath ++ "\n") temporaryDirectoryPath
+ inputFasta <- readFasta inputFastaFilePath
+ networkCheck <- checkNCBIConnection
+ if isLeft networkCheck
+ then do
+ putStrLn ("Error - Could not contact NCBI server: " ++ fromLeft networkCheck ++ "\n")
+ logMessage ("Error - Could not contact NCBI server: " ++ fromLeft networkCheck ++ "\n") temporaryDirectoryPath
+ else do
+ if null inputFasta
+ then do
+ putStrLn "Error: Input fasta file is empty."
+ logMessage "Error: Input fasta file is empty.\n" temporaryDirectoryPath
+ else do
+ let iterationNumber = 0
+ let tools = ["clustalo","mlocarna","RNAfold","RNAalifold","cmcalibrate","cmstat","cmbuild","RNAz"]
+ toolsCheck <- checkTools tools temporaryDirectoryPath
+ -- Check required commandline tools
+ if isLeft toolsCheck
+ then do
+ putStrLn ("Error - Not all required tools could be found in $PATH: " ++ fromLeft toolsCheck ++ "\n")
+ logMessage ("Error - Not all required tools could be found in $PATH: " ++ fromLeft toolsCheck ++ "\n") temporaryDirectoryPath
+ else do
+ logToolVersions temporaryDirectoryPath
+ let inputSequence = reformatFasta (head inputFasta)
+ initialTaxId <- setInitialTaxId inputBlastDatabase temporaryDirectoryPath inputTaxId inputSequence
+ let checkedTaxonomyRestriction = checkTaxonomyRestriction taxonomyRestriction
+ let staticOptions = StaticOptions temporaryDirectoryPath sessionId (fromJust inputnSCICutoff) inputTaxId singleHitperTax lengthFilter coverageFilter threads inputBlastDatabase checkedTaxonomyRestriction (setVerbose verboseLevel)
+ let initialization = ModelConstruction iterationNumber inputSequence [] initialTaxId Nothing (fromJust inputEvalueCutoff) False [] []
+ logMessage (show initialization) temporaryDirectoryPath
+ modelConstructionResults <- modelConstructer staticOptions initialization
+ let resultTaxonomyRecordsCSVTable = constructTaxonomyRecordsCSVTable modelConstructionResults
+ let resultNumber = (length (concatMap sequenceRecords (taxRecords modelConstructionResults))) + 1
+ resultEvaluation <- evaluateConstructionResult staticOptions resultNumber
+ appendFile (temporaryDirectoryPath ++ "Log") resultEvaluation
+ writeFile (temporaryDirectoryPath ++ "result.csv") resultTaxonomyRecordsCSVTable
+ resultSummary modelConstructionResults staticOptions
+ writeFile (temporaryDirectoryPath ++ "done") ""
diff --git a/src/Bio/RNAlienData.hs b/src/Bio/RNAlienData.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/RNAlienData.hs
@@ -0,0 +1,144 @@
+-- | This module contains data structures for RNAlien
+
+module Bio.RNAlienData where
+
+import qualified Data.ByteString.Lazy.Char8 as L
+import Bio.Core.Sequence
+import Bio.Sequence.Fasta
+import Bio.Taxonomy
+
+-- | Static construction options
+data StaticOptions = StaticOptions
+ { tempDirPath :: String,
+ sessionID :: String,
+ nSCICutoff :: Double,
+ userTaxId :: Maybe Int,
+ singleHitperTaxToggle :: Bool,
+ lengthFilterToggle :: Bool,
+ coverageFilterToggle :: Bool,
+ cpuThreads :: Int,
+ blastDatabase :: Maybe String,
+ taxRestriction :: Maybe String,
+ verbositySwitch :: Bool
+ } deriving (Show)
+
+-- | Keeps track of model construction
+data ModelConstruction = ModelConstruction
+ { iterationNumber :: Int,
+ inputFasta :: Sequence,
+ taxRecords :: [TaxonomyRecord],
+ --Taxonomy ID of the highest node in taxonomic subtree used in search
+ upperTaxonomyLimit :: Maybe Int,
+ taxonomicContext :: Maybe Taxon,
+ evalueThreshold :: Double,
+ alignmentModeInfernal :: Bool,
+ selectedQueries :: [String],
+ potentialMembers :: [SearchResult]
+ }
+
+instance Show ModelConstruction where
+ show (ModelConstruction _iterationNumber _inputFasta _taxRecords _upperTaxonomyLimit _taxonomicContext _evalueThreshold _alignmentModeInfernal _selectedQueries _potentialMembers) = a ++ b ++ c ++ d ++ e ++ g ++ h ++ i
+ where a = "Modelconstruction iteration: " ++ show _iterationNumber ++ "\n"
+ b = "Input fasta:\n" ++ L.unpack (unSL (seqheader _inputFasta)) ++ "\n" ++ L.unpack (unSD (seqdata _inputFasta)) ++ "\n"
+ c = show _taxRecords
+ d = "Upper taxonomy limit: " ++ maybe "not set" show _upperTaxonomyLimit ++ "\n"
+ e = "Taxonomic Context: " ++ maybe "not set" show _taxonomicContext ++ "\n"
+ g = "Evalue cutoff: " ++ show _evalueThreshold ++ "\n"
+ h = "Selected queries: \n" ++ concatMap (\x -> x ++ "\n") _selectedQueries
+ i = "Potential Members: \n" ++ concatMap show _potentialMembers
+
+data TaxonomyRecord = TaxonomyRecord
+ { recordTaxonomyId :: Int,
+ sequenceRecords :: [SequenceRecord]
+ }
+
+instance Show TaxonomyRecord where
+ show (TaxonomyRecord _recordTaxonomyId _sequenceRecords) = a ++ b
+ where a = "TaxonomyRecord TaxonomyId: " ++ show _recordTaxonomyId ++ "\n"
+ b = show _sequenceRecords
+
+data SequenceRecord = SequenceRecord
+ { --Sequence consisting of SeqLabel, and SeqData
+ nucleotideSequence :: Sequence,
+ -- 0 is unaligned, number is the iteration the sequence has been included into the alignment
+ aligned :: Int,
+ recordDescription :: String,
+ -- Is the sequence derived from the blast hit coordinates (B) or from a corresponding genbank feature (G)
+ sequenceOrigin :: Char
+ }
+
+instance Show SequenceRecord where
+ show (SequenceRecord _nucleotideSequence _aligned _recordDescription _sequenceOrigin) = a ++ b ++ c ++ d
+ where a = "SequenceRecord TaxonomyId: " ++ show _recordDescription ++ "\n"
+ b = "Sequence Origin: " ++ _recordDescription ++ "\n"
+ c = "Aligned in iteration: " ++ show _aligned ++ "\n"
+ d = "Sequence Origin: " ++ show _nucleotideSequence ++ "\n"
+-- |
+data CMsearch = CMsearch
+ { queryCMfile :: String,
+ targetSequenceDatabase :: String,
+ numberOfWorkerThreads :: String,
+ cmsearchHits :: [CMsearchHit]
+-- hitAlignments :: [CMsearchHitAlignment]
+-- internalCMPipelineStatisticsSummary
+ } deriving (Show, Eq, Read)
+
+-- |
+data CMsearchHit = CMsearchHit
+ { hitRank :: Int,
+ hitSignificance :: Char,
+ hitEvalue :: Double,
+ hitScore :: Double,
+ hitBias :: Double,
+ hitSequenceHeader :: L.ByteString,
+ hitStart :: Int,
+ hitEnd :: Int,
+ hitStrand :: Char,
+ hitModel :: L.ByteString,
+ hitTruncation :: L.ByteString,
+ hitGCContent :: Double,
+ hitDescription :: L.ByteString
+ } deriving (Show, Eq, Read)
+
+data SearchResult = SearchResult
+ { candidates :: [(Sequence,Int,String,Char)],
+ blastDatabaseSize :: Maybe Double
+ }
+
+instance Show SearchResult where
+ show (SearchResult _candidates _blastDatabaseSize) = a ++ b
+ where a = "SearchResults :\n " ++ concatMap show _candidates ++ "\n"
+ b = "BlastDb Size: " ++ show _blastDatabaseSize ++ "\n"
+
+-- |
+data CMstat = CMstat
+ { statIndex :: Int,
+ statName :: String,
+ statAccession :: String,
+ statSequenceNumber :: Int,
+ statEffectiveSequences :: Double,
+ statConsensusLength :: Int,
+ -- W The expected maximum length of a hit to the model.
+ statW :: Int,
+ statBasepairs :: Int,
+ statBifurcations :: Int,
+ statModel :: String,
+ relativeEntropyCM :: Double,
+ relativeEntropyHMM :: Double
+ } deriving (Eq, Read)
+
+instance Show CMstat where
+ show (CMstat _statIndex _statName _statAccession _statSequenceNumber _statEffectiveSequences _statConsensusLength _statW _statBasepairs _statBifurcations _statModel _relativeEntropyCM _relativeEntropyHMM) = a ++ b ++ c ++ d ++ e ++ f ++ g ++ h ++ i ++ j ++ k ++ l
+ where a = "CMstat - covariance model statistics:\nIndex: " ++ show _statIndex ++ "\n"
+ b = "Name: " ++ show _statName ++ "\n"
+ c = "Accession: " ++ show _statAccession ++ "\n"
+ d = "Sequence Number: " ++ show _statSequenceNumber ++ "\n"
+ e = "Effective Sequences: " ++ show _statEffectiveSequences ++ "\n"
+ f = "Consensus length: " ++ show _statConsensusLength ++ "\n"
+ g = "Expected maximum hit-length: " ++ show _statW ++ "\n"
+ h = "Basepairs: " ++ show _statBasepairs ++ "\n"
+ i = "Bifurcations: " ++ show _statBifurcations ++ "\n"
+ j = "Modeltype: " ++ show _statModel ++ "\n"
+ k = "Relative Entropy CM: " ++ show _relativeEntropyCM ++ "\n"
+ l = "Relative Entropy HMM: " ++ show _relativeEntropyHMM ++ "\n"
+
diff --git a/src/Bio/RNAlienLibrary.hs b/src/Bio/RNAlienLibrary.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/RNAlienLibrary.hs
@@ -0,0 +1,1819 @@
+-- | This module contains functions for RNAlien
+
+module Bio.RNAlienLibrary (
+ module Bio.RNAlienData,
+ createSessionID,
+ logMessage,
+ logEither,
+ modelConstructer,
+ constructTaxonomyRecordsCSVTable,
+ resultSummary,
+ setVerbose,
+ logToolVersions,
+ checkTools,
+ systemCMsearch,
+ readCMSearch,
+ compareCM,
+ parseCMSearch,
+ cmSearchsubString,
+ setInitialTaxId,
+ evaluateConstructionResult,
+ readCMstat,
+ parseCMstat,
+ checkNCBIConnection,
+ preprocessClustalForRNAz,
+ preprocessClustalForRNAzExternal,
+ rnaZEvalOutput,
+ reformatFasta,
+ checkTaxonomyRestriction,
+ evaluePartitionTrimCMsearchHits
+ )
+where
+
+import System.Process
+import qualified System.FilePath as FP
+import Text.ParserCombinators.Parsec
+import Data.List
+import Data.Char
+import Bio.Core.Sequence
+import Bio.Sequence.Fasta
+import Bio.BlastXML
+import Bio.ClustalParser
+import Data.Int (Int16)
+import Bio.RNAlienData
+import qualified Data.ByteString.Lazy.Char8 as L
+import Bio.Taxonomy
+import Data.Either.Unwrap
+import Data.Maybe
+import Bio.EntrezHTTP
+import qualified Data.List.Split as DS
+import System.Exit
+import Data.Either (lefts,rights)
+import qualified Text.EditDistance as ED
+import qualified Data.Vector as V
+import Control.Concurrent
+import System.Random
+import Data.Csv
+import Data.Matrix
+import Bio.BlastHTTP
+import Data.Clustering.Hierarchical
+import System.Directory
+import System.Console.CmdArgs
+import qualified Control.Exception.Base as CE
+import Bio.RNAfoldParser
+import Bio.RNAalifoldParser
+import Bio.RNAzParser
+import Network.HTTP
+
+-- | Initial RNA family model construction - generates iteration number, seed alignment and model
+modelConstructer :: StaticOptions -> ModelConstruction -> IO ModelConstruction
+modelConstructer staticOptions modelConstruction = do
+ logMessage ("Iteration: " ++ show (iterationNumber modelConstruction) ++ "\n") (tempDirPath staticOptions)
+ iterationSummary modelConstruction staticOptions
+ let currentIterationNumber = (iterationNumber modelConstruction)
+ let foundSequenceNumber = length (concatMap sequenceRecords (taxRecords modelConstruction))
+ --extract queries
+ let queries = extractQueries foundSequenceNumber modelConstruction
+ logVerboseMessage (verbositySwitch staticOptions) ("Queries:" ++ show queries ++ "\n") (tempDirPath staticOptions)
+ let iterationDirectory = (tempDirPath staticOptions) ++ (show currentIterationNumber) ++ "/"
+ let maybeLastTaxId = extractLastTaxId (taxonomicContext modelConstruction)
+ if (isNothing maybeLastTaxId) then logMessage ("Lineage: Could not extract last tax id \n") (tempDirPath staticOptions) else (return ())
+ --If highest node in linage was used as upper taxonomy limit, taxonomic tree is exhausted
+ if (maybe True (\uppertaxlimit -> maybe True (\lastTaxId -> uppertaxlimit /= lastTaxId) maybeLastTaxId) (upperTaxonomyLimit modelConstruction))
+ then do
+ createDirectory (iterationDirectory)
+ let (upperTaxLimit,lowerTaxLimit) = setTaxonomicContextEntrez currentIterationNumber (taxonomicContext modelConstruction) (upperTaxonomyLimit modelConstruction)
+ logVerboseMessage (verbositySwitch staticOptions) ("Upper taxonomy limit: " ++ (show upperTaxLimit) ++ "\n " ++ "Lower taxonomy limit: "++ show lowerTaxLimit ++ "\n") (tempDirPath staticOptions)
+ --search queries
+ let expectThreshold = setBlastExpectThreshold modelConstruction
+ searchResults <- catchAll (searchCandidates staticOptions Nothing currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold queries)
+ (\e -> do logWarning ("Warning: Search results iteration" ++ show (iterationNumber modelConstruction) ++ " - exception: " ++ show e) (tempDirPath staticOptions)
+ return (SearchResult [] Nothing))
+ currentTaxonomicContext <- getTaxonomicContextEntrez upperTaxLimit (taxonomicContext modelConstruction)
+ if null (candidates searchResults)
+ then do
+ alignmentConstructionWithoutCandidates currentTaxonomicContext upperTaxLimit staticOptions modelConstruction
+ else do
+ alignmentConstructionWithCandidates currentTaxonomicContext upperTaxLimit searchResults staticOptions modelConstruction
+ else do
+ logMessage ("Message: Modelconstruction complete: Out of queries or taxonomic tree exhausted\n") (tempDirPath staticOptions)
+ modelConstructionResult staticOptions modelConstruction
+
+catchAll :: IO a -> (CE.SomeException -> IO a) -> IO a
+catchAll = CE.catch
+
+setInitialTaxId :: Maybe String -> String -> Maybe Int -> Sequence -> IO (Maybe Int)
+setInitialTaxId inputBlastDatabase tempdir inputTaxId inputSequence = do
+ if (isNothing inputTaxId)
+ then do
+ initialTaxId <- findTaxonomyStart inputBlastDatabase tempdir inputSequence
+ return (Just initialTaxId)
+ else do
+ return inputTaxId
+
+extractLastTaxId :: Maybe Taxon -> Maybe Int
+extractLastTaxId taxon
+ | isNothing taxon = Nothing
+ | V.null lineageExVector = Nothing
+ | otherwise = Just (lineageTaxId (V.head lineageExVector))
+ where lineageExVector = V.fromList (lineageEx (fromJust taxon))
+
+modelConstructionResult :: StaticOptions -> ModelConstruction -> IO ModelConstruction
+modelConstructionResult staticOptions modelConstruction = do
+ let currentIterationNumber = iterationNumber modelConstruction
+ let outputDirectory = tempDirPath staticOptions
+ logMessage ("Global search iteration: " ++ show currentIterationNumber ++ "\n") outputDirectory
+ iterationSummary modelConstruction staticOptions
+ let foundSequenceNumber = length (concatMap sequenceRecords (taxRecords modelConstruction))
+ --extract queries
+ --let querySeqIds = selectedQueries modelConstruction ---
+ let queries = extractQueries foundSequenceNumber modelConstruction ---
+ --let alignedSequences' = map nucleotideSequence (concatMap sequenceRecords (taxRecords modelConstruction)) ---
+ logVerboseMessage (verbositySwitch staticOptions) ("Queries:" ++ show queries ++ "\n") outputDirectory
+ let iterationDirectory = outputDirectory ++ (show currentIterationNumber) ++ "/"
+ createDirectory (iterationDirectory)
+ let logFileDirectoryPath = iterationDirectory ++ "log"
+ createDirectoryIfMissing False logFileDirectoryPath
+ let expectThreshold = setBlastExpectThreshold modelConstruction
+ (alignmentResults,currentPotentialMembers) <- if isJust (taxRestriction staticOptions)
+ then do
+ --taxonomic restriction
+ let (upperTaxLimit,lowerTaxLimit) = setRestrictedTaxonomyLimits (fromJust (taxRestriction staticOptions))
+ restrictedCandidates <- catchAll (searchCandidates staticOptions (taxRestriction staticOptions) currentIterationNumber upperTaxLimit lowerTaxLimit expectThreshold queries)
+ (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory
+ return (SearchResult [] Nothing))
+ let uniqueCandidates = filterDuplicates modelConstruction restrictedCandidates
+ (restrictedAlignmentResults,restrictedPotentialMembers) <- catchAll (alignCandidates staticOptions modelConstruction (fromJust (taxRestriction staticOptions)) uniqueCandidates)
+ (\e -> do logWarning ("Warning: Alignment results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory
+ return ([],[]))
+ let currentPotentialMembers = [SearchResult restrictedPotentialMembers (blastDatabaseSize restrictedCandidates)]
+ return (restrictedAlignmentResults,currentPotentialMembers)
+ else do
+ --taxonomic context archea
+ let (upperTaxLimit1,lowerTaxLimit1) = (Just (2157 :: Int), Nothing)
+ candidates1 <- catchAll (searchCandidates staticOptions (Just "archea") currentIterationNumber upperTaxLimit1 lowerTaxLimit1 expectThreshold queries)
+ (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory
+ return (SearchResult [] Nothing))
+ let uniqueCandidates1 = filterDuplicates modelConstruction candidates1
+ (alignmentResults1,potentialMembers1) <- catchAll (alignCandidates staticOptions modelConstruction "archea" uniqueCandidates1)
+ (\e -> do logWarning ("Warning: Alignment results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory
+ return ([],[]))
+ --taxonomic context bacteria
+ let (upperTaxLimit2,lowerTaxLimit2) = (Just (2 :: Int), Nothing)
+ candidates2 <- catchAll (searchCandidates staticOptions (Just "bacteria") currentIterationNumber upperTaxLimit2 lowerTaxLimit2 expectThreshold queries)
+ (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory
+ return (SearchResult [] Nothing))
+ let uniqueCandidates2 = filterDuplicates modelConstruction candidates2
+ (alignmentResults2,potentialMembers2) <- catchAll (alignCandidates staticOptions modelConstruction "bacteria" uniqueCandidates2)
+ (\e -> do logWarning ("Warning: Alignment results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory
+ return ([],[]))
+ --taxonomic context eukaryia
+ let (upperTaxLimit3,lowerTaxLimit3) = (Just (2759 :: Int), Nothing)
+ candidates3 <- catchAll (searchCandidates staticOptions (Just "eukaryia") currentIterationNumber upperTaxLimit3 lowerTaxLimit3 expectThreshold queries)
+ (\e -> do logWarning ("Warning: Search results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory
+ return (SearchResult [] Nothing))
+ let uniqueCandidates3 = filterDuplicates modelConstruction candidates3
+ (alignmentResults3,potentialMembers3) <- catchAll (alignCandidates staticOptions modelConstruction "eukaryia" uniqueCandidates3)
+ (\e -> do logWarning ("Warning: Alignment results iteration" ++ show currentIterationNumber ++ " - exception: " ++ show e) outputDirectory
+ return ([],[]))
+ let alignmentResults = alignmentResults1 ++ alignmentResults2 ++ alignmentResults3
+ let currentPotentialMembers = [SearchResult potentialMembers1 (blastDatabaseSize candidates1), SearchResult potentialMembers2 (blastDatabaseSize candidates2), SearchResult potentialMembers3 (blastDatabaseSize candidates3)]
+ return (alignmentResults,currentPotentialMembers)
+ let preliminaryFastaPath = iterationDirectory ++ "model.fa"
+ let preliminaryCMPath = iterationDirectory ++ "model.cm"
+ let preliminaryAlignmentPath = iterationDirectory ++ "model.stockholm"
+ let preliminaryCMLogPath = iterationDirectory ++ "model.cm.log"
+ let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults Nothing Nothing [] currentPotentialMembers (alignmentModeInfernal modelConstruction)
+ if (length alignmentResults == 0) && (not (alignmentModeInfernal modelConstruction))
+ then do
+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment result initial mode\n") outputDirectory
+ logMessage ("Message: No sequences found that statisfy filters. Try to reconstruct model with less strict cutoff parameters.") outputDirectory
+ let alignedSequences = extractAlignedSequences (iterationNumber modelConstruction) modelConstruction
+ let alignmentSequences = map snd (V.toList (V.concat [alignedSequences]))
+ writeFasta preliminaryFastaPath alignmentSequences
+ let cmBuildFilepath = iterationDirectory ++ "model" ++ ".cmbuild"
+ let foldFilepath = iterationDirectory ++ "model" ++ ".fold"
+ _ <- systemRNAfold preliminaryFastaPath foldFilepath
+ foldoutput <- readRNAfold foldFilepath
+ let seqStructure = foldSecondaryStructure (fromRight foldoutput)
+ let stockholAlignment = convertFastaFoldStockholm (head alignmentSequences) seqStructure
+ writeFile preliminaryAlignmentPath stockholAlignment
+ _ <- systemCMbuild preliminaryAlignmentPath preliminaryCMPath cmBuildFilepath
+ _ <- systemCMcalibrate "fast" (cpuThreads staticOptions) preliminaryCMPath preliminaryCMLogPath
+ resultModelConstruction <- reevaluatePotentialMembers staticOptions nextModelConstructionInput
+ return resultModelConstruction
+ else do
+ if (alignmentModeInfernal modelConstruction)
+ then do
+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - infernal mode\n") outputDirectory
+ constructModel nextModelConstructionInput staticOptions
+ writeFile (iterationDirectory ++ "done") ""
+ logMessage (iterationSummaryLog nextModelConstructionInput) outputDirectory
+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) outputDirectory
+ resultModelConstruction <- reevaluatePotentialMembers staticOptions nextModelConstructionInput
+ return resultModelConstruction
+ else do
+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - initial mode\n") outputDirectory
+ constructModel nextModelConstructionInput staticOptions
+ let nextModelConstructionInputInfernalMode = nextModelConstructionInput {alignmentModeInfernal = True}
+ logMessage (iterationSummaryLog nextModelConstructionInputInfernalMode) outputDirectory
+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputInfernalMode) outputDirectory
+ writeFile (iterationDirectory ++ "done") ""
+ resultModelConstruction <- reevaluatePotentialMembers staticOptions nextModelConstructionInputInfernalMode
+ return resultModelConstruction
+
+-- | Reevaluate collected potential members for inclusion in the result model
+reevaluatePotentialMembers :: StaticOptions -> ModelConstruction -> IO ModelConstruction
+reevaluatePotentialMembers staticOptions modelConstruction = do
+ let currentIterationNumber = iterationNumber modelConstruction
+ let outputDirectory = tempDirPath staticOptions
+ iterationSummary modelConstruction staticOptions
+ logMessage ("Reevaluation of potential members iteration: " ++ show currentIterationNumber ++ "\n") outputDirectory
+ let iterationDirectory = outputDirectory ++ (show currentIterationNumber) ++ "/"
+ createDirectory (iterationDirectory)
+ let indexedPotentialMembers = V.indexed (V.fromList (potentialMembers modelConstruction))
+ potentialMembersAlignmentResultVector <- V.mapM (\(i,searchresult) -> (alignCandidates staticOptions modelConstruction (show i ++ "_") searchresult)) indexedPotentialMembers
+ let potentialMembersAlignmentResults = V.toList potentialMembersAlignmentResultVector
+ let alignmentResults = concatMap fst potentialMembersAlignmentResults
+ let discardedMembers = concatMap snd potentialMembersAlignmentResults
+ writeFile (outputDirectory ++ "log/discarded") (concatMap show discardedMembers)
+ let resultFastaPath = outputDirectory ++ "result.fa"
+ let resultCMPath = outputDirectory ++ "result.cm"
+ let resultAlignmentPath = outputDirectory ++ "result.stockholm"
+ let resultCMLogPath = outputDirectory ++ "log/result.cm.log"
+ if null alignmentResults
+ then do
+ let lastIterationFastaPath = outputDirectory ++ show (currentIterationNumber - 1)++ "/model.fa"
+ let lastIterationAlignmentPath = outputDirectory ++ show (currentIterationNumber - 1) ++ "/model.stockholm"
+ let lastIterationCMPath = outputDirectory ++ show (currentIterationNumber - 1)++ "/model.cm"
+ copyFile lastIterationCMPath resultCMPath
+ copyFile lastIterationFastaPath resultFastaPath
+ copyFile lastIterationAlignmentPath resultAlignmentPath
+ _ <- systemCMcalibrate "standard" (cpuThreads staticOptions) resultCMPath resultCMLogPath
+ writeFile (iterationDirectory ++ "done") ""
+ return modelConstruction
+ else do
+ let lastIterationFastaPath = outputDirectory ++ show currentIterationNumber ++ "/model.fa"
+ let lastIterationAlignmentPath = outputDirectory ++ show currentIterationNumber ++ "/model.stockholm"
+ let lastIterationCMPath = outputDirectory ++ show currentIterationNumber ++ "/model.cm"
+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - infernal mode\n") outputDirectory
+ let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults Nothing Nothing [] [] (alignmentModeInfernal modelConstruction)
+ constructModel nextModelConstructionInput staticOptions
+ copyFile lastIterationCMPath resultCMPath
+ copyFile lastIterationFastaPath resultFastaPath
+ copyFile lastIterationAlignmentPath resultAlignmentPath
+ logMessage (iterationSummaryLog nextModelConstructionInput) outputDirectory
+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) outputDirectory
+ _ <- systemCMcalibrate "standard" (cpuThreads staticOptions) resultCMPath resultCMLogPath
+ writeFile (iterationDirectory ++ "done") ""
+ return nextModelConstructionInput
+
+---------------------------------------------------------
+
+alignmentConstructionWithCandidates :: Maybe Taxon -> Maybe Int -> SearchResult -> StaticOptions -> ModelConstruction -> IO ModelConstruction
+alignmentConstructionWithCandidates currentTaxonomicContext currentUpperTaxonomyLimit searchResults staticOptions modelConstruction = do
+ --candidates usedUpperTaxonomyLimit blastDatabaseSize
+ let currentIterationNumber = (iterationNumber modelConstruction)
+ let iterationDirectory = (tempDirPath staticOptions) ++ (show currentIterationNumber) ++ "/"
+ --let usedUpperTaxonomyLimit = (snd (head candidates))
+ --align search result
+ (alignmentResults,potentialMemberEntries) <- catchAll (alignCandidates staticOptions modelConstruction "" searchResults)
+ (\e -> do logWarning ("Warning: Alignment results iteration" ++ show (iterationNumber modelConstruction) ++ " - exception: " ++ show e) (tempDirPath staticOptions)
+ return ([],[]))
+ let currentPotentialMembers = [SearchResult potentialMemberEntries (blastDatabaseSize searchResults)]
+ if (length alignmentResults == 0) && (not (alignmentModeInfernal modelConstruction))
+ then do
+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - length 1 - inital mode" ++ "\n") (tempDirPath staticOptions)
+ --too few sequences for alignment. because of lack in sequences no cm was constructed before
+ --reusing previous modelconstruction with increased upperTaxonomyLimit but include found sequence
+ --prepare next iteration
+ let newTaxEntries = (taxRecords modelConstruction) ++ (buildTaxRecords alignmentResults currentIterationNumber)
+ let nextModelConstructionInputWithThreshold = modelConstruction {iterationNumber = (currentIterationNumber + 1),upperTaxonomyLimit = currentUpperTaxonomyLimit, taxRecords = newTaxEntries,taxonomicContext = currentTaxonomicContext}
+ logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)
+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) ----
+ writeFile (iterationDirectory ++ "done") ""
+ nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithThreshold
+ return nextModelConstruction
+ else do
+ --select queries
+ currentSelectedQueries <- selectQueries staticOptions modelConstruction alignmentResults
+ if (alignmentModeInfernal modelConstruction)
+ then do
+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - infernal mode\n") (tempDirPath staticOptions)
+ --prepare next iteration
+ let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext currentSelectedQueries currentPotentialMembers True
+ constructModel nextModelConstructionInput staticOptions
+ writeFile (iterationDirectory ++ "done") ""
+ logMessage (iterationSummaryLog nextModelConstructionInput) (tempDirPath staticOptions)
+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInput) (tempDirPath staticOptions) ----
+ nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInput
+ return nextModelConstruction
+ else do
+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction with candidates - initial mode\n") (tempDirPath staticOptions)
+ --First round enough candidates are available for modelconstruction, alignmentModeInfernal is set to true after this iteration
+ --prepare next iteration
+ let nextModelConstructionInput = constructNext currentIterationNumber modelConstruction alignmentResults currentUpperTaxonomyLimit currentTaxonomicContext currentSelectedQueries currentPotentialMembers False
+ constructModel nextModelConstructionInput staticOptions
+ let nextModelConstructionInputWithInfernalMode = nextModelConstructionInput {alignmentModeInfernal = True}
+ logMessage (iterationSummaryLog nextModelConstructionInputWithInfernalMode) (tempDirPath staticOptions)
+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithInfernalMode) (tempDirPath staticOptions) ----
+ writeFile (iterationDirectory ++ "done") ""
+ nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithInfernalMode
+ return nextModelConstruction
+
+alignmentConstructionWithoutCandidates :: Maybe Taxon -> Maybe Int -> StaticOptions -> ModelConstruction -> IO ModelConstruction
+alignmentConstructionWithoutCandidates currentTaxonomicContext upperTaxLimit staticOptions modelConstruction = do
+ let currentIterationNumber = (iterationNumber modelConstruction)
+ let iterationDirectory = (tempDirPath staticOptions) ++ (show currentIterationNumber) ++ "/"
+ --Found no new candidates in this iteration, reusing previous modelconstruction with increased upperTaxonomyLimit
+ let nextModelConstructionInputWithThreshold = modelConstruction {iterationNumber = (currentIterationNumber + 1),upperTaxonomyLimit = upperTaxLimit,taxonomicContext = currentTaxonomicContext}
+ --copy model and alignment from last iteration in place if present
+ let previousIterationCMPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber - 1)) ++ "/model.cm"
+ previousIterationCMexists <- doesFileExist previousIterationCMPath
+ if previousIterationCMexists
+ then do
+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction no candidates - previous cm\n") (tempDirPath staticOptions)
+ let previousIterationFastaPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber - 1)) ++ "/model.fa"
+ let previousIterationAlignmentPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber - 1)) ++ "/model.stockholm"
+ let thisIterationFastaPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber)) ++ "/model.fa"
+ let thisIterationAlignmentPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber)) ++ "/model.stockholm"
+ let thisIterationCMPath = (tempDirPath staticOptions) ++ (show (currentIterationNumber)) ++ "/model.cm"
+ copyFile previousIterationFastaPath thisIterationFastaPath
+ copyFile previousIterationAlignmentPath thisIterationAlignmentPath
+ copyFile previousIterationCMPath thisIterationCMPath
+ logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)
+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)
+ writeFile (iterationDirectory ++ "done") ""
+ nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithThreshold
+ return nextModelConstruction
+ else do
+ logVerboseMessage (verbositySwitch staticOptions) ("Alignment construction no candidates - no previous iteration cm\n") (tempDirPath staticOptions)
+ logMessage (iterationSummaryLog nextModelConstructionInputWithThreshold) (tempDirPath staticOptions)
+ logVerboseMessage (verbositySwitch staticOptions) (show nextModelConstructionInputWithThreshold) (tempDirPath staticOptions) ----
+ writeFile (iterationDirectory ++ "done") ""
+ nextModelConstruction <- modelConstructer staticOptions nextModelConstructionInputWithThreshold
+ return nextModelConstruction
+
+findTaxonomyStart :: Maybe String -> String -> Sequence -> IO Int
+findTaxonomyStart inputBlastDatabase temporaryDirectory querySequence = do
+ let queryIndexString = "1"
+ let hitNumberQuery = buildHitNumberQuery "&HITLIST_SIZE=10"
+ let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at"
+ let blastQuery = BlastHTTPQuery (Just "ncbi") (Just "blastn") inputBlastDatabase [querySequence] (Just (hitNumberQuery ++ registrationInfo)) (Just (5400000000 :: Int))
+ logMessage ("No tax id provided - Sending find taxonomy start blast query \n") temporaryDirectory
+ blastOutput <- CE.catch (blastHTTP blastQuery)
+ (\e -> do let err = show (e :: CE.IOException)
+ logWarning ("Warning: Blast attempt failed:" ++ " " ++ err) temporaryDirectory
+ error "findTaxonomyStart: Blast attempt failed"
+ return (Left ""))
+ let logFileDirectoryPath = temporaryDirectory ++ "taxonomystart" ++ "/"
+ createDirectory logFileDirectoryPath
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_1blastOutput") (show blastOutput)
+ logEither blastOutput temporaryDirectory
+ let blastHitsArePresent = either (\_ -> False) blastMatchesPresent blastOutput
+ if (blastHitsArePresent)
+ then do
+ let rightBlast = fromRight blastOutput
+ let bestHit = getBestHit rightBlast
+ bestBlastHitTaxIdOutput <- retrieveBlastHitTaxIdEntrez [bestHit]
+ let taxIdFromEntrySummaries = extractTaxIdFromEntrySummaries (snd bestBlastHitTaxIdOutput)
+ if (null taxIdFromEntrySummaries) then (error "findTaxonomyStart: - head: empty list of taxonomy entry summary for best hit") else return ()
+ let rightBestTaxIdResult = head taxIdFromEntrySummaries
+ logMessage ("Initial TaxId: " ++ (show rightBestTaxIdResult) ++ "\n") temporaryDirectory
+ CE.evaluate rightBestTaxIdResult
+ else error "Find taxonomy start: Could not find blast hits to use as a taxonomic starting point"
+
+searchCandidates :: StaticOptions -> Maybe String -> Int -> Maybe Int -> Maybe Int -> Double -> [Sequence] -> IO SearchResult
+searchCandidates staticOptions finaliterationprefix iterationnumber upperTaxLimit lowerTaxLimit expectThreshold querySequences' = do
+ --let fastaSeqData = seqdata _querySequence
+ if (null querySequences') then error "searchCandidates: - head: empty list of query sequences" else return ()
+ let queryLength = fromIntegral (seqlength (head querySequences'))
+ let queryIndexString = "1"
+ let entrezTaxFilter = buildTaxFilterQuery upperTaxLimit lowerTaxLimit
+ logVerboseMessage (verbositySwitch staticOptions) ("entrezTaxFilter" ++ show entrezTaxFilter ++ "\n") (tempDirPath staticOptions)
+ let hitNumberQuery = buildHitNumberQuery "&HITLIST_SIZE=5000&EXPECT=" ++ show expectThreshold
+ let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at"
+ let blastQuery = BlastHTTPQuery (Just "ncbi") (Just "blastn") (blastDatabase staticOptions) querySequences' (Just (hitNumberQuery ++ entrezTaxFilter ++ registrationInfo)) (Just (5400000000 :: Int))
+ logVerboseMessage (verbositySwitch staticOptions) ("Sending blast query " ++ (show iterationnumber) ++ "\n") (tempDirPath staticOptions)
+ blastOutput <- CE.catch (blastHTTP blastQuery)
+ (\e -> do let err = show (e :: CE.IOException)
+ logWarning ("Warning: Blast attempt failed:" ++ " " ++ err) (tempDirPath staticOptions)
+ return (Left ""))
+ let logFileDirectoryPath = (tempDirPath staticOptions) ++ (show iterationnumber) ++ "/" ++ (fromMaybe "" finaliterationprefix) ++ "log"
+ logDirectoryPresent <- doesDirectoryExist logFileDirectoryPath
+ if (not logDirectoryPresent)
+ then createDirectory (logFileDirectoryPath) else return ()
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_1blastOutput") (show blastOutput)
+ logEither blastOutput (tempDirPath staticOptions)
+ let blastHitsArePresent = either (const False) blastMatchesPresent blastOutput
+ if blastHitsArePresent
+ then do
+ let rightBlast = fromRight blastOutput
+ -- bestBlastHitTaxIdOutput <- retrieveBlastHitTaxIdEntrez [bestHit]
+ -- let taxIdFromEntrySummaries = extractTaxIdFromEntrySummaries bestBlastHitTaxIdOutput
+ -- if (null taxIdFromEntrySummaries) then (error "searchCandidates: - head: empty list of taxonomy entry summary for best hit") else return ()
+ -- let rightBestTaxIdResult = head taxIdFromEntrySummaries
+ -- logVerboseMessage (verbositySwitch staticOptions) ("rightbestTaxIdResult: " ++ (show rightBestTaxIdResult) ++ "\n") (tempDirPath staticOptions)
+ let blastHits = concatMap hits (results rightBlast)
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_2blastHits") (showlines blastHits)
+ --filter by length
+ let blastHitsFilteredByLength = filterByHitLength blastHits queryLength (lengthFilterToggle staticOptions)
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_3blastHitsFilteredByLength") (showlines blastHitsFilteredByLength)
+ let blastHitsFilteredByCoverage = filterByCoverage blastHitsFilteredByLength queryLength (coverageFilterToggle staticOptions)
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_3ablastHitsFilteredByLength") (showlines blastHitsFilteredByCoverage)
+ --tag BlastHits with TaxId
+ blastHitsWithTaxIdOutput <- retrieveBlastHitsTaxIdEntrez blastHitsFilteredByLength
+ let uncheckedBlastHitsWithTaxIdList = map (\(blasthits,taxIdout) -> (blasthits,extractTaxIdFromEntrySummaries taxIdout)) blastHitsWithTaxIdOutput
+ let checkedBlastHitsWithTaxId = filter (\(_,taxids) -> not (null taxids)) uncheckedBlastHitsWithTaxIdList
+ --todo checked blasthittaxidswithblasthits need to be merged as taxid blasthit pairs
+ let blastHitsWithTaxId = zip (concatMap fst checkedBlastHitsWithTaxId) (concatMap snd checkedBlastHitsWithTaxId)
+ blastHitsWithParentTaxIdOutput <- retrieveParentTaxIdsEntrez blastHitsWithTaxId
+ --let blastHitsWithParentTaxId = concat blastHitsWithParentTaxIdOutput
+ -- filter by ParentTaxId (only one hit per TaxId)
+ let blastHitsFilteredByParentTaxIdWithParentTaxId = filterByParentTaxId blastHitsWithParentTaxIdOutput True
+ --let blastHitsFilteredByParentTaxId = map fst blastHitsFilteredByParentTaxIdWithParentTaxId
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_4blastHitsFilteredByParentTaxId") (showlines blastHitsFilteredByParentTaxIdWithParentTaxId)
+ -- Filtering with TaxTree (only hits from the same subtree as besthit)
+ --let blastHitsWithTaxId = zip blastHitsFilteredByParentTaxId blastHittaxIdList
+ --let (_, filteredBlastResults) = filterByNeighborhoodTreeConditional alignndmentModeInfernalToggle upperTaxLimit blastHitsWithTaxId (inputTaxNodes staticOptions) (fromJust upperTaxLimit) (singleHitperTaxToggle staticOptions)
+ --writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_5filteredBlastResults") (showlines filteredBlastResults)
+ -- Coordinate generation
+ let nonEmptyfilteredBlastResults = filter (\(blasthit,_) -> not (null (matches blasthit))) blastHitsFilteredByParentTaxIdWithParentTaxId
+ let requestedSequenceElements = map (getRequestedSequenceElement queryLength) nonEmptyfilteredBlastResults
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_6requestedSequenceElements") (showlines requestedSequenceElements)
+ -- Retrieval of full sequences from entrez
+ --fullSequencesWithSimilars <- retrieveFullSequences requestedSequenceElements
+ fullSequencesWithSimilars <- retrieveFullSequences staticOptions requestedSequenceElements
+ if null fullSequencesWithSimilars
+ then do
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_10afullSequencesWithSimilars") "No sequences retrieved"
+ CE.evaluate (SearchResult [] Nothing)
+ else do
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_10afullSequencesWithSimilars") (showlines fullSequencesWithSimilars)
+ let fullSequences = filterIdenticalSequences fullSequencesWithSimilars 100
+ let fullSequencesWithOrigin = map (\(parsedFasta,taxid,seqSubject) -> (parsedFasta,taxid,seqSubject,'B')) fullSequences
+ writeFile (logFileDirectoryPath ++ "/" ++ queryIndexString ++ "_10fullSequences") (showlines fullSequences)
+ let maybeFractionEvalueMatch = getHitWithFractionEvalue rightBlast
+ if isNothing maybeFractionEvalueMatch
+ then CE.evaluate (SearchResult [] Nothing)
+ else do
+ let fractionEvalueMatch = fromJust maybeFractionEvalueMatch
+ let dbSize = computeDataBaseSize (e_val fractionEvalueMatch) (bits fractionEvalueMatch) (fromIntegral queryLength ::Double)
+ CE.evaluate (SearchResult fullSequencesWithOrigin (Just dbSize))
+ else CE.evaluate (SearchResult [] Nothing)
+
+-- |Computes size of blast db in Mb
+computeDataBaseSize :: Double -> Double -> Double -> Double
+computeDataBaseSize evalue bitscore querylength = ((evalue * 2 ** bitscore) / querylength)/10^(6 :: Integer)
+
+alignCandidates :: StaticOptions -> ModelConstruction -> String -> SearchResult -> IO ([(Sequence,Int,String,Char)],[(Sequence,Int,String,Char)])
+alignCandidates staticOptions modelConstruction multipleSearchResultPrefix searchResults = do
+ let iterationDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/" ++ multipleSearchResultPrefix
+ createDirectoryIfMissing False (iterationDirectory ++ "log")
+ if null (candidates searchResults)
+ then do
+ writeFile (iterationDirectory ++ "log" ++ "/11candidates") ("No candidates to align")
+ return ([],[])
+ else do
+ --refilter for similarity
+ writeFile (iterationDirectory ++ "log" ++ "/11candidates") (showlines (candidates searchResults))
+ let alignedSequences = map snd (V.toList (extractAlignedSequences (iterationNumber modelConstruction) modelConstruction))
+ let filteredCandidates = filterWithCollectedSequences (candidates searchResults) alignedSequences 99
+ writeFile (iterationDirectory ++ "log" ++ "/12candidatesFilteredByCollected") (showlines filteredCandidates)
+ if alignmentModeInfernal modelConstruction
+ then alignCandidatesInfernalMode staticOptions modelConstruction multipleSearchResultPrefix (blastDatabaseSize searchResults) filteredCandidates
+ else alignCandidatesInitialMode staticOptions modelConstruction multipleSearchResultPrefix filteredCandidates
+
+alignCandidatesInfernalMode :: StaticOptions -> ModelConstruction -> String -> Maybe Double -> [(Sequence,Int,String,Char)] -> IO ([(Sequence,Int,String,Char)],[(Sequence,Int,String,Char)])
+alignCandidatesInfernalMode staticOptions modelConstruction multipleSearchResultPrefix blastDbSize filteredCandidates = do
+ let iterationDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/" ++ multipleSearchResultPrefix
+ let candidateSequences = extractCandidateSequences filteredCandidates
+ logVerboseMessage (verbositySwitch staticOptions) "Alignment Mode Infernal\n" (tempDirPath staticOptions)
+ let indexedCandidateSequenceList = (V.toList candidateSequences)
+ let cmSearchFastaFilePaths = map (constructFastaFilePaths iterationDirectory) indexedCandidateSequenceList
+ let cmSearchFilePaths = map (constructCMsearchFilePaths iterationDirectory) indexedCandidateSequenceList
+ let covarianceModelPath = tempDirPath staticOptions ++ show (iterationNumber modelConstruction - 1) ++ "/" ++ "model.cm"
+ mapM_ (\(number,_nucleotideSequence) -> writeFasta (iterationDirectory ++ show number ++ ".fa") [_nucleotideSequence]) indexedCandidateSequenceList
+ let zippedFastaCMSearchResultPaths = zip cmSearchFastaFilePaths cmSearchFilePaths
+ --check with cmSearch
+ mapM_ (uncurry (systemCMsearch (cpuThreads staticOptions) ("-Z " ++ show (fromJust blastDbSize)) covarianceModelPath)) zippedFastaCMSearchResultPaths
+ cmSearchResults <- mapM readCMSearch cmSearchFilePaths
+ writeFile (iterationDirectory ++ "cm_error") (concatMap show (lefts cmSearchResults))
+ let rightCMSearchResults = rights cmSearchResults
+ let cmSearchCandidatesWithSequences = zip rightCMSearchResults filteredCandidates
+ let (trimmedSelectedCandidates,potentialCandidates,rejectedCandidates) = evaluePartitionTrimCMsearchHits (evalueThreshold modelConstruction) cmSearchCandidatesWithSequences
+ createDirectoryIfMissing False (iterationDirectory ++ "log")
+ writeFile (iterationDirectory ++ "log" ++ "/13selectedCandidates'") (showlines trimmedSelectedCandidates)
+ writeFile (iterationDirectory ++ "log" ++ "/14rejectedCandidates'") (showlines rejectedCandidates)
+ writeFile (iterationDirectory ++ "log" ++ "/15potentialCandidates'") (showlines potentialCandidates)
+ CE.evaluate (map snd trimmedSelectedCandidates,map snd potentialCandidates)
+
+alignCandidatesInitialMode :: StaticOptions -> ModelConstruction -> String -> [(Sequence,Int,String,Char)] -> IO ([(Sequence,Int,String,Char)],[(Sequence,Int,String,Char)])
+alignCandidatesInitialMode staticOptions modelConstruction multipleSearchResultPrefix filteredCandidates = do
+ let iterationDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/" ++ multipleSearchResultPrefix
+ --writeFile (iterationDirectory ++ "log" ++ "/11bcandidates") (showlines filteredCandidates)
+ createDirectoryIfMissing False (iterationDirectory ++ "log")
+ let candidateSequences = extractCandidateSequences filteredCandidates
+ --writeFile (iterationDirectory ++ "log" ++ "/11ccandidates") (showlines (V.toList candidateSequences))
+ logVerboseMessage (verbositySwitch staticOptions) "Alignment Mode Initial\n" (tempDirPath staticOptions)
+ --write Fasta sequences
+ let inputFastaFilepath = iterationDirectory ++ "input.fa"
+ let inputFoldFilepath = iterationDirectory ++ "input.fold"
+ writeFasta (iterationDirectory ++ "input.fa") ([inputFasta modelConstruction])
+ logMessage (showlines (V.toList candidateSequences)) (tempDirPath staticOptions)
+ V.mapM_ (\(number,nucleotideSequence') -> writeFasta (iterationDirectory ++ show number ++ ".fa") [nucleotideSequence']) candidateSequences
+ let candidateFastaFilepath = V.toList (V.map (\(number,_) -> iterationDirectory ++ show number ++ ".fa") candidateSequences)
+ let candidateFoldFilepath = V.toList (V.map (\(number,_) -> iterationDirectory ++ show number ++ ".fold") candidateSequences)
+ let locarnainClustalw2FormatFilepath = V.toList (V.map (\(number,_) -> iterationDirectory ++ show number ++ "." ++ "clustalmlocarna") candidateSequences)
+ let candidateAliFoldFilepath = V.toList (V.map (\(number,_) -> iterationDirectory ++ show number ++ ".alifold") candidateSequences)
+ let locarnaFilepath = V.toList (V.map (\(number,_) -> iterationDirectory ++ show number ++ "." ++ "mlocarna") candidateSequences)
+ alignSequences "locarna" " --write-structure --free-endgaps=++-- " (replicate (V.length candidateSequences) inputFastaFilepath) candidateFastaFilepath locarnainClustalw2FormatFilepath locarnaFilepath
+ --compute SequenceIdentities
+ let sequenceIdentities = V.map (\(_,s) -> sequenceIdentity (inputFasta modelConstruction) s/(100 :: Double)) candidateSequences
+ --compute SCI
+ systemRNAfold inputFastaFilepath inputFoldFilepath
+ inputfoldResult <- readRNAfold inputFoldFilepath
+ let inputFoldMFE = foldingEnergy (fromRight inputfoldResult)
+ mapM_ (uncurry systemRNAfold) (zip candidateFastaFilepath candidateFoldFilepath)
+ foldResults <- mapM readRNAfold candidateFoldFilepath
+ let candidateMFEs = map (foldingEnergy . fromRight) foldResults
+ let averageMFEs = map (\candidateMFE -> (candidateMFE + inputFoldMFE)/2) candidateMFEs
+ mapM_ (uncurry (systemRNAalifold "")) (zip locarnainClustalw2FormatFilepath candidateAliFoldFilepath)
+ alifoldResults <- mapM readRNAalifold candidateAliFoldFilepath
+ let consensusMFE = map (alignmentConsensusMinimumFreeEnergy . fromRight) alifoldResults
+ let sciidfraction = map (\(consMFE,averMFE,seqId) -> (consMFE/averMFE)/seqId) (zip3 consensusMFE averageMFEs (V.toList sequenceIdentities))
+ let idlog = concatMap (\(sciidfraction',consMFE,averMFE,seqId) -> show sciidfraction' ++ "," ++ show consMFE ++ "," ++ show averMFE ++ "," ++ show seqId ++ "\n")(zip4 sciidfraction consensusMFE averageMFEs (V.toList sequenceIdentities))
+ writeFile (iterationDirectory ++ "log" ++ "/idlog") (idlog)
+ let alignedCandidates = zip sciidfraction filteredCandidates
+ writeFile (iterationDirectory ++ "log" ++ "/zscores") (showlines alignedCandidates)
+ let (selectedCandidates,rejectedCandidates) = partition (\(sciidfraction',_) -> sciidfraction' > nSCICutoff staticOptions) alignedCandidates
+ mapM_ print (zip3 consensusMFE averageMFEs (V.toList sequenceIdentities))
+ writeFile (iterationDirectory ++ "log" ++ "/13selectedCandidates") (showlines selectedCandidates)
+ writeFile (iterationDirectory ++ "log" ++ "/14rejectedCandidates") (showlines rejectedCandidates)
+ CE.evaluate (map snd selectedCandidates,[])
+
+setClusterNumber :: Int -> Int
+setClusterNumber x
+ | x <= 5 = x
+ | otherwise = 5
+
+findCutoffforClusterNumber :: Dendrogram a -> Int -> Distance -> Distance
+findCutoffforClusterNumber clustaloDendrogram numberOfClusters currentCutoff
+ | currentClusterNumber >= numberOfClusters = currentCutoff
+ | otherwise = findCutoffforClusterNumber clustaloDendrogram numberOfClusters (currentCutoff-0.01)
+ where currentClusterNumber = length (cutAt clustaloDendrogram currentCutoff)
+
+selectQueries :: StaticOptions -> ModelConstruction -> [(Sequence,Int,String,Char)] -> IO [String]
+selectQueries staticOptions modelConstruction selectedCandidates = do
+ logVerboseMessage (verbositySwitch staticOptions) "SelectQueries\n" (tempDirPath staticOptions)
+ --Extract sequences from modelconstruction
+ let alignedSequences = extractAlignedSequences (iterationNumber modelConstruction) modelConstruction
+ let candidateSequences = extractQueryCandidates selectedCandidates
+ let iterationDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/"
+ let alignmentSequences = map snd (V.toList (V.concat [candidateSequences,alignedSequences]))
+ if length alignmentSequences > 3
+ then do
+ --write Fasta sequences
+ writeFasta (iterationDirectory ++ "query" ++ ".fa") alignmentSequences
+ let fastaFilepath = iterationDirectory ++ "query" ++ ".fa"
+ let clustaloFilepath = iterationDirectory ++ "query" ++ ".clustalo"
+ let clustaloDistMatrixPath = iterationDirectory ++ "query" ++ ".matrix"
+ alignSequences "clustalo" ("--full --distmat-out=" ++ clustaloDistMatrixPath ++ " ") [fastaFilepath] [] [clustaloFilepath] []
+ idsDistancematrix <- readClustaloDistMatrix clustaloDistMatrixPath
+ logEither idsDistancematrix (tempDirPath staticOptions)
+ let (clustaloIds,clustaloDistMatrix) = fromRight idsDistancematrix
+ logVerboseMessage (verbositySwitch staticOptions) ("Clustalid: " ++ intercalate "," clustaloIds ++ "\n") (tempDirPath staticOptions)
+ logVerboseMessage (verbositySwitch staticOptions) ("Distmatrix: " ++ show clustaloDistMatrix ++ "\n") (tempDirPath staticOptions)
+ let clustaloDendrogram = dendrogram UPGMA clustaloIds (getDistanceMatrixElements clustaloIds clustaloDistMatrix)
+ logVerboseMessage (verbositySwitch staticOptions) ("ClustaloDendrogram: " ++ show clustaloDendrogram ++ "\n") (tempDirPath staticOptions)
+ logVerboseMessage (verbositySwitch staticOptions) ("ClustaloDendrogram: " ++ show clustaloDistMatrix ++ "\n") (tempDirPath staticOptions)
+ let numberOfClusters = setClusterNumber (length alignmentSequences)
+ logVerboseMessage (verbositySwitch staticOptions) ("numberOfClusters: " ++ show numberOfClusters ++ "\n") (tempDirPath staticOptions)
+ let dendrogramStartCutDistance = 1 :: Double
+ let dendrogramCutDistance' = findCutoffforClusterNumber clustaloDendrogram numberOfClusters dendrogramStartCutDistance
+ logVerboseMessage (verbositySwitch staticOptions) ("dendrogramCutDistance': " ++ show dendrogramCutDistance' ++ "\n") (tempDirPath staticOptions)
+ let cutDendrogram = cutAt clustaloDendrogram dendrogramCutDistance'
+ --putStrLn "cutDendrogram: "
+ --print cutDendrogram
+ let currentSelectedQueries = take 5 (concatMap (take 1 . elements) cutDendrogram)
+ logVerboseMessage (verbositySwitch staticOptions) ("SelectedQueries: " ++ show currentSelectedQueries ++ "\n") (tempDirPath staticOptions)
+ writeFile (tempDirPath staticOptions ++ show (iterationNumber modelConstruction) ++ "/log" ++ "/13selectedQueries") (showlines currentSelectedQueries)
+ CE.evaluate currentSelectedQueries
+ else return []
+
+constructModel :: ModelConstruction -> StaticOptions -> IO String
+constructModel modelConstruction staticOptions = do
+ --Extract sequences from modelconstruction
+ let alignedSequences = extractAlignedSequences (iterationNumber modelConstruction) modelConstruction
+ --The CM resides in the iteration directory where its input alignment originates from
+ let outputDirectory = tempDirPath staticOptions ++ show (iterationNumber modelConstruction - 1) ++ "/"
+ let alignmentSequences = map snd (V.toList (V.concat [alignedSequences]))
+ --write Fasta sequences
+ writeFasta (outputDirectory ++ "model" ++ ".fa") alignmentSequences
+ let fastaFilepath = outputDirectory ++ "model" ++ ".fa"
+ let locarnaFilepath = outputDirectory ++ "model" ++ ".mlocarna"
+ let stockholmFilepath = outputDirectory ++ "model" ++ ".stockholm"
+ --- let reformatedClustalFilepath = outputDirectory ++ "model" ++ ".clustal.reformated"
+ let updatedStructureStockholmFilepath = outputDirectory ++ "newstructuremodel" ++ ".stockholm"
+ let cmalignCMFilepath = (tempDirPath staticOptions) ++ show (iterationNumber modelConstruction - 2) ++ "/" ++ "model" ++ ".cm"
+ let cmFilepath = outputDirectory ++ "model" ++ ".cm"
+ let cmCalibrateFilepath = outputDirectory ++ "model" ++ ".cmcalibrate"
+ let cmBuildFilepath = outputDirectory ++ "model" ++ ".cmbuild"
+ let alifoldFilepath = outputDirectory ++ "model" ++ ".alifold"
+ if alignmentModeInfernal modelConstruction
+ then do
+ logVerboseMessage (verbositySwitch staticOptions) "Construct Model - infernal mode\n" (tempDirPath staticOptions)
+ systemCMalign ("--cpu " ++ show (cpuThreads staticOptions)) cmalignCMFilepath fastaFilepath stockholmFilepath
+ systemRNAalifold "-r --cfactor 0.6 --nfactor 0.5" stockholmFilepath alifoldFilepath
+ replaceStatus <- replaceStockholmStructure stockholmFilepath alifoldFilepath updatedStructureStockholmFilepath
+ if null replaceStatus
+ then do
+ systemCMbuild updatedStructureStockholmFilepath cmFilepath cmBuildFilepath
+ systemCMcalibrate "fast" (cpuThreads staticOptions) cmFilepath cmCalibrateFilepath
+ return cmFilepath
+ else do
+ logWarning ("Warning: A problem occured updating the secondary structure of iteration " ++ show (iterationNumber modelConstruction) ++ " stockholm alignment: " ++ replaceStatus) (tempDirPath staticOptions)
+ systemCMbuild stockholmFilepath cmFilepath cmBuildFilepath
+ systemCMcalibrate "fast" (cpuThreads staticOptions) cmFilepath cmCalibrateFilepath
+ return cmFilepath
+ else do
+ logVerboseMessage (verbositySwitch staticOptions) "Construct Model - initial mode\n" (tempDirPath staticOptions)
+ alignSequences "mlocarna" ("--threads=" ++ show (cpuThreads staticOptions) ++ " ") [fastaFilepath] [] [locarnaFilepath] []
+ mlocarnaAlignment <- readStructuralClustalAlignment locarnaFilepath
+ logEither mlocarnaAlignment (tempDirPath staticOptions)
+ let stockholAlignment = convertClustaltoStockholm (fromRight mlocarnaAlignment)
+ writeFile stockholmFilepath stockholAlignment
+ _ <- systemCMbuild stockholmFilepath cmFilepath cmBuildFilepath
+ _ <- systemCMcalibrate "fast" (cpuThreads staticOptions) cmFilepath cmCalibrateFilepath
+ return cmFilepath
+
+-- | Replaces structure of input stockholm file with the consensus structure of alifoldFilepath and outputs updated stockholmfile
+replaceStockholmStructure :: String -> String -> String -> IO String
+replaceStockholmStructure stockholmFilepath alifoldFilepath updatedStructureStockholmFilepath = do
+ inputAln <- readFile stockholmFilepath
+ inputRNAalifold <- readRNAalifold alifoldFilepath
+ if isLeft inputRNAalifold
+ then do
+ return (show (fromLeft inputRNAalifold))
+ else do
+ let alifoldstructure = alignmentConsensusDotBracket (fromRight inputRNAalifold)
+ let seedLinesVector = V.fromList (lines inputAln)
+ let structureIndices = V.toList (V.findIndices isStructureLine seedLinesVector)
+ let updatedStructureElements = updateStructureElements seedLinesVector alifoldstructure structureIndices
+ let newVector = seedLinesVector V.// updatedStructureElements
+ let newVectorString = concatMap (++ "\n") (V.toList newVector)
+ writeFile updatedStructureStockholmFilepath newVectorString
+ return []
+
+updateStructureElements :: V.Vector String -> String -> [Int] -> [(Int,String)]
+updateStructureElements inputVector structureString indices
+ | null indices = []
+ | otherwise = newElement ++ updateStructureElements inputVector (drop structureLength structureString) (tail indices)
+ where currentIndex = head indices
+ currentElement = inputVector V.! currentIndex
+ elementLength = length currentElement
+ structureStartIndex = maximum (elemIndices ' ' currentElement) + 1
+ structureLength = elementLength - structureStartIndex
+ newElementHeader = take structureStartIndex currentElement
+ newElementStructure = take structureLength structureString
+ newElement = [(currentIndex,newElementHeader ++ newElementStructure)]
+
+isStructureLine :: String -> Bool
+isStructureLine input = "#=GC SS_cons" `isInfixOf` input
+
+-- Generates iteration string for Log
+iterationSummaryLog :: ModelConstruction -> String
+iterationSummaryLog mC = output
+ where upperTaxonomyLimitOutput = maybe "not set" show (upperTaxonomyLimit mC)
+ output = "Upper taxonomy id limit: " ++ upperTaxonomyLimitOutput ++ ", Collected members: " ++ show (length (concatMap sequenceRecords (taxRecords mC))) ++ "\n"
+
+-- | Used for passing progress to Alien server
+iterationSummary :: ModelConstruction -> StaticOptions -> IO()
+iterationSummary mC sO = do
+ --iteration -- tax limit -- bitscore cutoff -- blastresult -- aligned seqs --queries --fa link --aln link --cm link
+ let upperTaxonomyLimitOutput = maybe "not set" show (upperTaxonomyLimit mC)
+ let output = show (iterationNumber mC) ++ "," ++ upperTaxonomyLimitOutput ++ "," ++ show (length (concatMap sequenceRecords (taxRecords mC)))
+ writeFile ((tempDirPath sO) ++ "/log/" ++ show (iterationNumber mC) ++ ".log") output
+
+-- | Used for passing progress to Alien server
+resultSummary :: ModelConstruction -> StaticOptions -> IO()
+resultSummary mC sO = do
+ --iteration -- tax limit -- bitscore cutoff -- blastresult -- aligned seqs --queries --fa link --aln link --cm link
+ let upperTaxonomyLimitOutput = maybe "not set" show (upperTaxonomyLimit mC)
+ let output = show (iterationNumber mC) ++ "," ++ upperTaxonomyLimitOutput ++ "," ++ show (length (concatMap sequenceRecords (taxRecords mC)))
+ writeFile (tempDirPath sO ++ "/log/result" ++ ".log") output
+
+readClustaloDistMatrix :: String -> IO (Either ParseError ([String],Matrix Double))
+readClustaloDistMatrix = parseFromFile genParserClustaloDistMatrix
+
+genParserClustaloDistMatrix :: GenParser Char st ([String],Matrix Double)
+genParserClustaloDistMatrix = do
+ _ <- many1 digit
+ newline
+ clustaloDistRow <- many1 (try genParserClustaloDistRow)
+ eof
+ return (map fst clustaloDistRow,fromLists (map snd clustaloDistRow))
+
+genParserClustaloDistRow :: GenParser Char st (String,[Double])
+genParserClustaloDistRow = do
+ entryId <- many1 (noneOf " ")
+ many1 space
+ distances <- many1 (try genParserClustaloDistance)
+ newline
+ return (entryId,distances)
+
+genParserClustaloDistance :: GenParser Char st Double
+genParserClustaloDistance = do
+ distance <- many1 (oneOf "1234567890.")
+ optional (try (char ' ' ))
+ return (readDouble distance)
+
+getDistanceMatrixElements :: [String] -> Matrix Double -> String -> String -> Double
+getDistanceMatrixElements ids distMatrix id1 id2 = distance
+ -- Data.Matrix is indexed starting with 1
+ where indexid1 = fromJust (elemIndex id1 ids) + 1
+ indexid2 = fromJust (elemIndex id2 ids) + 1
+ distance = getElem indexid1 indexid2 distMatrix
+
+-- | Filter duplicates removes hits in sequences that were already collected. This happens during revisiting the starting subtree.
+filterDuplicates :: ModelConstruction -> SearchResult -> SearchResult
+filterDuplicates modelConstruction inputSearchResult = uniqueSearchResult
+ where alignedSequences = map snd (V.toList (extractAlignedSequences (iterationNumber modelConstruction) modelConstruction))
+ collectedIdentifiers = map seqid alignedSequences
+ uniques = filter (\(s,_,_,_) -> notElem (seqid s) collectedIdentifiers) (candidates inputSearchResult)
+ uniqueSearchResult = SearchResult uniques (blastDatabaseSize inputSearchResult)
+
+-- | Filter a list of similar extended blast hits
+--filterIdenticalSequencesWithOrigin :: [(Sequence,Int,String,Char)] -> Double -> [(Sequence,Int,String,Char)]
+--filterIdenticalSequencesWithOrigin (headSequence:rest) identitycutoff = result
+-- where filteredSequences = filter (\x -> (sequenceIdentity (firstOfQuadruple headSequence) (firstOfQuadruple x)) < identitycutoff) rest
+-- result = headSequence:(filterIdenticalSequencesWithOrigin filteredSequences identitycutoff)
+--filterIdenticalSequencesWithOrigin [] _ = []
+
+-- | Filter a list of similar extended blast hits
+filterIdenticalSequences :: [(Sequence,Int,String)] -> Double -> [(Sequence,Int,String)]
+filterIdenticalSequences (headSequence:rest) identitycutoff = result
+ where filteredSequences = filter (\x -> sequenceIdentity (firstOfTriple headSequence) (firstOfTriple x) < identitycutoff) rest
+ result = headSequence:filterIdenticalSequences filteredSequences identitycutoff
+filterIdenticalSequences [] _ = []
+
+-- | Filter sequences too similar to already aligned sequences
+filterWithCollectedSequences :: [(Sequence,Int,String,Char)] -> [Sequence] -> Double -> [(Sequence,Int,String,Char)]
+filterWithCollectedSequences inputCandidates collectedSequences identitycutoff = filter (isUnSimilarSequence collectedSequences identitycutoff . firstOfQuadruple) inputCandidates
+--filterWithCollectedSequences [] [] _ = []
+
+-- | Filter alignment entries by similiarity
+filterIdenticalAlignmentEntry :: [ClustalAlignmentEntry] -> Double -> [ClustalAlignmentEntry]
+filterIdenticalAlignmentEntry (headEntry:rest) identitycutoff = result
+ where filteredEntries = filter (\x -> (stringIdentity (entryAlignedSequence headEntry) (entryAlignedSequence x)) < identitycutoff) rest
+ result = headEntry:filterIdenticalAlignmentEntry filteredEntries identitycutoff
+filterIdenticalAlignmentEntry [] _ = []
+
+
+isUnSimilarSequence :: [Sequence] -> Double -> Sequence -> Bool
+isUnSimilarSequence collectedSequences identitycutoff checkSequence = any (\ x -> (sequenceIdentity checkSequence x) < identitycutoff) collectedSequences
+
+
+firstOfTriple :: (t, t1, t2) -> t
+firstOfTriple (a,_,_) = a
+
+firstOfQuadruple :: (t, t1, t2, t3) -> t
+firstOfQuadruple (a,_,_,_) = a
+
+-- | Check if the result field of BlastResult is filled and if hits are present
+blastMatchesPresent :: BlastResult -> Bool
+blastMatchesPresent blastResult
+ | null resultList = False
+ | otherwise = True
+ where resultList = concatMap matches (concatMap hits (results blastResult))
+
+-- | Compute identity of sequences
+stringIdentity :: String -> String -> Double
+stringIdentity string1 string2 = identityPercent
+ where distance = ED.levenshteinDistance ED.defaultEditCosts string1 string2
+ maximumDistance = maximum [length string1,length string2]
+ identityPercent = 100 - ((fromIntegral distance/fromIntegral maximumDistance) * (read "100" ::Double))
+
+-- | Compute identity of sequences
+sequenceIdentity :: Sequence -> Sequence -> Double
+sequenceIdentity sequence1 sequence2 = identityPercent
+ where distance = ED.levenshteinDistance ED.defaultEditCosts sequence1string sequence2string
+ sequence1string = L.unpack (unSD (seqdata sequence1))
+ sequence2string = L.unpack (unSD (seqdata sequence2))
+ maximumDistance = maximum [length sequence1string,length sequence2string]
+ identityPercent = 100 - ((fromIntegral distance/fromIntegral maximumDistance) * (read "100" ::Double))
+
+getTaxonomicContextEntrez :: Maybe Int -> Maybe Taxon -> IO (Maybe Taxon)
+getTaxonomicContextEntrez upperTaxLimit currentTaxonomicContext = do
+ if isJust upperTaxLimit
+ then if isJust currentTaxonomicContext
+ then return currentTaxonomicContext
+ else retrieveTaxonomicContextEntrez (fromJust upperTaxLimit)
+ --return retrievedTaxonomicContext
+ else return Nothing
+
+setTaxonomicContextEntrez :: Int -> Maybe Taxon -> Maybe Int -> (Maybe Int, Maybe Int)
+setTaxonomicContextEntrez currentIterationNumber currentTaxonomicContext subTreeTaxId
+ | currentIterationNumber == 0 = (subTreeTaxId, Nothing)
+ | otherwise = setUpperLowerTaxLimitEntrez (fromJust subTreeTaxId) (fromJust currentTaxonomicContext)
+
+-- setTaxonomic Context for next candidate search, the upper bound of the last search become the lower bound of the next
+setUpperLowerTaxLimitEntrez :: Int -> Taxon -> (Maybe Int, Maybe Int)
+setUpperLowerTaxLimitEntrez subTreeTaxId currentTaxonomicContext = (upperLimit,lowerLimit)
+ where upperLimit = raiseTaxIdLimitEntrez subTreeTaxId currentTaxonomicContext
+ lowerLimit = Just subTreeTaxId
+
+raiseTaxIdLimitEntrez :: Int -> Taxon -> Maybe Int
+raiseTaxIdLimitEntrez subTreeTaxId taxon = parentNodeTaxId
+ where lastUpperBoundNodeIndex = fromJust (V.findIndex (\node -> lineageTaxId node == subTreeTaxId) lineageExVector)
+ linageNodeTaxId = Just (lineageTaxId (lineageExVector V.! (lastUpperBoundNodeIndex -1)))
+ lineageExVector = V.fromList (lineageEx taxon)
+ --the input taxid is not part of the lineage, therefor we look for further taxids in the lineage after we used the parent tax id of the input node
+ parentNodeTaxId = if subTreeTaxId == taxonTaxId taxon then Just (taxonParentTaxId taxon) else linageNodeTaxId
+
+constructNext :: Int -> ModelConstruction -> [(Sequence, Int, String, Char)] -> Maybe Int -> Maybe Taxon -> [String] -> [SearchResult] -> Bool -> ModelConstruction
+constructNext currentIterationNumber modelconstruction alignmentResults upperTaxLimit inputTaxonomicContext inputSelectedQueries inputPotentialMembers toggleInfernalAlignmentModeTrue = nextModelConstruction
+ where newIterationNumber = currentIterationNumber + 1
+ taxEntries = taxRecords modelconstruction ++ buildTaxRecords alignmentResults currentIterationNumber
+ potMembers = potentialMembers modelconstruction ++ inputPotentialMembers
+ currentAlignmentMode = toggleInfernalAlignmentModeTrue || alignmentModeInfernal modelconstruction
+ nextModelConstruction = ModelConstruction newIterationNumber (inputFasta modelconstruction) taxEntries upperTaxLimit inputTaxonomicContext (evalueThreshold modelconstruction) currentAlignmentMode inputSelectedQueries potMembers
+
+buildTaxRecords :: [(Sequence,Int,String,Char)] -> Int -> [TaxonomyRecord]
+buildTaxRecords alignmentResults currentIterationNumber = taxonomyRecords
+ where taxIdGroups = groupBy sameTaxIdAlignmentResult alignmentResults
+ taxonomyRecords = map (buildTaxRecord currentIterationNumber) taxIdGroups
+
+sameTaxIdAlignmentResult :: (Sequence,Int,String,Char) -> (Sequence,Int,String,Char) -> Bool
+sameTaxIdAlignmentResult (_,taxId1,_,_) (_,taxId2,_,_) = taxId1 == taxId2
+
+buildTaxRecord :: Int -> [(Sequence,Int,String,Char)] -> TaxonomyRecord
+buildTaxRecord currentIterationNumber entries = taxRecord
+ where recordTaxId = (\(_,currentTaxonomyId,_,_) -> currentTaxonomyId) (head entries)
+ seqRecords = map (buildSeqRecord currentIterationNumber) entries
+ taxRecord = TaxonomyRecord recordTaxId seqRecords
+
+buildSeqRecord :: Int -> (Sequence,Int,String,Char) -> SequenceRecord
+buildSeqRecord currentIterationNumber (parsedFasta,_,seqSubject,seqOrigin) = SequenceRecord parsedFasta currentIterationNumber seqSubject seqOrigin
+
+-- | Partitions sequences by containing a cmsearch hit and extracts the hit region as new sequence
+evaluePartitionTrimCMsearchHits :: Double -> [(CMsearch,(Sequence, Int, String, Char))] -> ([(CMsearch,(Sequence, Int, String, Char))],[(CMsearch,(Sequence, Int, String, Char))],[(CMsearch,(Sequence, Int, String, Char))])
+evaluePartitionTrimCMsearchHits eValueThreshold cmSearchCandidatesWithSequences = (trimmedSelectedCandidates,potentialCandidates,rejectedCandidates)
+ where potentialMemberseValueThreshold = eValueThreshold * 1000
+ (prefilteredCandidates,rejectedCandidates) = partition (\(cmSearchResult,_) -> any (\hitScore' -> potentialMemberseValueThreshold >= hitEvalue hitScore') (cmsearchHits cmSearchResult)) cmSearchCandidatesWithSequences
+ (selectedCandidates,potentialCandidates) = partition (\(cmSearchResult,_) -> any (\hitScore' -> eValueThreshold >= hitEvalue hitScore') (cmsearchHits cmSearchResult)) prefilteredCandidates
+ trimmedSelectedCandidates = map (\(cmSearchResult,inputSequence) -> (cmSearchResult,trimCMsearchHit cmSearchResult inputSequence)) selectedCandidates
+
+
+trimCMsearchHit :: CMsearch -> (Sequence, Int, String, Char) -> (Sequence, Int, String, Char)
+trimCMsearchHit cmSearchResult (inputSequence,b,c,d) = (subSequence,b,c,d)
+ where hitScoreEntry = head (cmsearchHits cmSearchResult)
+ sequenceString = L.unpack (unSD (seqdata inputSequence))
+ sequenceSubstring = cmSearchsubString (hitStart hitScoreEntry) (hitEnd hitScoreEntry) sequenceString
+ --extend original seqheader
+ newSequenceHeader = L.pack (L.unpack (unSL (seqheader inputSequence)) ++ "cmS_" ++ show (hitStart hitScoreEntry) ++ "_" ++ show (hitEnd hitScoreEntry) ++ "_" ++ show (hitStrand hitScoreEntry))
+ subSequence = Seq (SeqLabel newSequenceHeader) (SeqData (L.pack sequenceSubstring)) Nothing
+
+-- | Extract a substring with coordinates from cmsearch, first nucleotide has index 1
+cmSearchsubString :: Int -> Int -> String -> String
+cmSearchsubString startSubString endSubString inputString
+ | startSubString < endSubString = take (endSubString - (startSubString -1))(drop (startSubString - 1) inputString)
+ | startSubString > endSubString = take (reverseEnd - (reverseStart - 1))(drop (reverseStart - 1 ) (reverse inputString))
+ | otherwise = take (endSubString - (startSubString -1))(drop (startSubString - 1) inputString)
+ where stringLength = length inputString
+ reverseStart = stringLength - (startSubString + 1)
+ reverseEnd = stringLength - (endSubString - 1)
+
+extractQueries :: Int -> ModelConstruction -> [Sequence]
+extractQueries foundSequenceNumber modelconstruction
+ | foundSequenceNumber < 3 = [fastaSeqData]
+ | otherwise = querySequences'
+ where fastaSeqData = inputFasta modelconstruction
+ querySeqIds = selectedQueries modelconstruction
+ alignedSequences = fastaSeqData:map nucleotideSequence (concatMap sequenceRecords (taxRecords modelconstruction))
+ querySequences' = concatMap (\querySeqId -> filter (\alignedSeq -> L.unpack (unSL (seqid alignedSeq)) == querySeqId) alignedSequences) querySeqIds
+
+extractQueryCandidates :: [(Sequence,Int,String,Char)] -> V.Vector (Int,Sequence)
+extractQueryCandidates querycandidates = indexedSeqences
+ where sequences = map (\(candidateSequence,_,_,_) -> candidateSequence) querycandidates
+ indexedSeqences = V.map (\(number,candidateSequence) -> (number + 1,candidateSequence))(V.indexed (V.fromList sequences))
+
+buildTaxFilterQuery :: Maybe Int -> Maybe Int -> String
+buildTaxFilterQuery upperTaxLimit lowerTaxLimit
+ | isNothing upperTaxLimit = ""
+ | isNothing lowerTaxLimit = "&ENTREZ_QUERY=" ++ encodedTaxIDQuery (fromJust upperTaxLimit)
+ | otherwise = "&ENTREZ_QUERY=" ++ "%28txid" ++ show (fromJust upperTaxLimit) ++ "%5BORGN%5D%29" ++ "NOT" ++ "%28txid" ++ show (fromJust lowerTaxLimit) ++ "%5BORGN%5D&EQ_OP%29"
+
+buildHitNumberQuery :: String -> String
+buildHitNumberQuery hitNumber
+ | hitNumber == "" = ""
+ | otherwise = "&ALIGNMENTS=" ++ hitNumber
+
+encodedTaxIDQuery :: Int -> String
+encodedTaxIDQuery taxID = "txid" ++ show taxID ++ "%20%5BORGN%5D&EQ_OP"
+
+-- | Adds cm prefix to pseudo random number
+randomid :: Int16 -> String
+randomid number = "cm" ++ show number
+
+createSessionID :: Maybe String -> IO String
+createSessionID sessionIdentificator = do
+ if isJust sessionIdentificator
+ then return (fromJust sessionIdentificator)
+ else do
+ randomNumber <- randomIO :: IO Int16
+ let sessionId = randomid (abs (randomNumber))
+ return sessionId
+
+-- | Run external locarna command and read the output into the corresponding datatype
+systemlocarna :: String -> (String,String,String,String) -> IO ExitCode
+systemlocarna options (inputFilePath1, inputFilePath2, clustalformatoutputFilePath, outputFilePath) = system ("locarna " ++ options ++ " --clustal=" ++ clustalformatoutputFilePath ++ " " ++ inputFilePath1 ++ " " ++ inputFilePath2 ++ " > " ++ outputFilePath)
+
+-- | Run external mlocarna command and read the output into the corresponding datatype, there is also a folder created at the location of the input fasta file
+systemMlocarna :: String -> (String,String) -> IO ExitCode
+systemMlocarna options (inputFilePath, outputFilePath) = system ("mlocarna " ++ options ++ " " ++ inputFilePath ++ " > " ++ outputFilePath)
+
+-- | Run external mlocarna command and read the output into the corresponding datatype, there is also a folder created at the location of the input fasta file, the job is terminated after the timeout provided in seconds
+systemMlocarnaWithTimeout :: String -> String -> (String,String) -> IO ExitCode
+systemMlocarnaWithTimeout timeout options (inputFilePath, outputFilePath) = system ("timeout " ++ timeout ++"s "++ "mlocarna " ++ options ++ " " ++ inputFilePath ++ " > " ++ outputFilePath)
+
+-- | Run external clustalo command and return the Exitcode
+systemClustalw2 :: String -> (String,String,String) -> IO ExitCode
+systemClustalw2 options (inputFilePath, outputFilePath, summaryFilePath) = system ("clustalw2 " ++ options ++ "-INFILE=" ++ inputFilePath ++ " -OUTFILE=" ++ outputFilePath ++ ">" ++ summaryFilePath)
+
+-- | Run external clustalo command and return the Exitcode
+systemClustalo :: String -> (String,String) -> IO ExitCode
+systemClustalo options (inputFilePath, outputFilePath) = system ("clustalo " ++ options ++ "--infile=" ++ inputFilePath ++ " >" ++ outputFilePath)
+
+-- | Run external CMbuild command and read the output into the corresponding datatype
+systemCMbuild :: String -> String -> String -> IO ExitCode
+systemCMbuild alignmentFilepath modelFilepath outputFilePath = system ("cmbuild " ++ modelFilepath ++ " " ++ alignmentFilepath ++ " > " ++ outputFilePath)
+
+-- | Run CMCompare and read the output into the corresponding datatype
+systemCMcompare :: String -> String -> String -> IO ExitCode
+systemCMcompare model1path model2path outputFilePath = system ("CMCompare -q " ++ model1path ++ " " ++ model2path ++ " >" ++ outputFilePath)
+
+-- | Run CMsearch
+systemCMsearch :: Int -> String -> String -> String -> String -> IO ExitCode
+systemCMsearch cpus options covarianceModelPath sequenceFilePath outputPath = system ("cmsearch --notrunc --cpu " ++ show cpus ++ " " ++ options ++ " -g " ++ covarianceModelPath ++ " " ++ sequenceFilePath ++ "> " ++ outputPath)
+
+-- | Run CMstat
+systemCMstat :: String -> String -> IO ExitCode
+systemCMstat covarianceModelPath outputPath = system ("cmstat " ++ covarianceModelPath ++ " > " ++ outputPath)
+
+-- | Run CMcalibrate and return exitcode
+systemCMcalibrate :: String -> Int -> String -> String -> IO ExitCode
+systemCMcalibrate mode cpus covarianceModelPath outputPath
+ | mode == "fast" = system ("cmcalibrate --beta 1E-4 --cpu " ++ show cpus ++ " " ++ covarianceModelPath ++ "> " ++ outputPath)
+ | otherwise = system ("cmcalibrate --cpu " ++ show cpus ++ " " ++ covarianceModelPath ++ "> " ++ outputPath)
+
+
+-- | Run CMcalibrate and return exitcode
+systemCMalign :: String -> String -> String -> String -> IO ExitCode
+systemCMalign options filePathCovarianceModel filePathSequence filePathAlignment = system ("cmalign " ++ options ++ " " ++ filePathCovarianceModel ++ " " ++ filePathSequence ++ "> " ++ filePathAlignment)
+
+compareCM :: String -> String -> String -> IO Double
+compareCM rfamCMPath resultCMpath outputDirectory = do
+ let myOptions = defaultDecodeOptions {
+ decDelimiter = fromIntegral (ord ' ')
+ }
+ let rfamCMFileName = FP.takeBaseName rfamCMPath
+ let resultCMFileName = FP.takeBaseName resultCMpath
+ let cmcompareResultPath = outputDirectory ++ rfamCMFileName ++ resultCMFileName ++ ".cmcompare"
+ _ <- systemCMcompare rfamCMPath resultCMpath cmcompareResultPath
+ inputCMcompare <- readFile cmcompareResultPath
+ let singlespaceCMcompare = unwords(words inputCMcompare)
+ let decodedCmCompareOutput = head (V.toList (fromRight (decodeWith myOptions NoHeader (L.pack singlespaceCMcompare) :: Either String (V.Vector [String]))))
+ --two.cm three.cm 27.996 19.500 CCCAAAGGGCCCAAAGGG (((...)))(((...))) (((...)))(((...))) [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] [11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]
+ let bitscore1 = read (decodedCmCompareOutput !! 2) :: Double
+ let bitscore2 = read (decodedCmCompareOutput !! 3) :: Double
+ let minmax = minimum [bitscore1,bitscore2]
+ return minmax
+
+readInt :: String -> Int
+readInt = read
+
+readDouble :: String -> Double
+readDouble = read
+
+-- | parse from input filePath
+parseCMSearch :: String -> Either ParseError CMsearch
+parseCMSearch = parse genParserCMsearch "parseCMsearch"
+
+-- | parse from input filePath
+readCMSearch :: String -> IO (Either ParseError CMsearch)
+readCMSearch filePath = do
+ parsedFile <- parseFromFile genParserCMsearch filePath
+ CE.evaluate parsedFile
+
+genParserCMsearch :: GenParser Char st CMsearch
+genParserCMsearch = do
+ string "# cmsearch :: search CM(s) against a sequence database"
+ newline
+ string "# INFERNAL "
+ many1 (noneOf "\n")
+ newline
+ string "# Copyright (C) 201"
+ many1 (noneOf "\n")
+ newline
+ string "# Freely distributed under the GNU General Public License (GPLv3)."
+ newline
+ string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"
+ newline
+ string "# query CM file:"
+ many1 space
+ queryCMfile' <- many1 (noneOf "\n")
+ newline
+ string "# target sequence database:"
+ many1 space
+ targetSequenceDatabase' <- many1 (noneOf "\n")
+ newline
+ optional (try (genParserCMsearchHeaderField "# CM configuration"))
+ optional (try (genParserCMsearchHeaderField "# database size is set to"))
+ optional (try (genParserCMsearchHeaderField "# truncated sequence detection"))
+ string "# number of worker threads:"
+ many1 space
+ numberOfWorkerThreads' <- many1 (noneOf "\n")
+ newline
+ string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"
+ newline
+ optional newline
+ string "Query:"
+ many1 (noneOf "\n")
+ newline
+ optional (try (genParserCMsearchHeaderField "Accession"))
+ optional (try (genParserCMsearchHeaderField "Description"))
+ string "Hit scores:"
+ newline
+ choice [try (string " rank"), try (string " rank") , try (string " rank"), try (string " rank"),try (string " rank"),try (string " rank")]
+ many1 space
+ string "E-value"
+ many1 space
+ string "score"
+ many1 space
+ string "bias"
+ many1 space
+ string "sequence"
+ many1 space
+ string "start"
+ many1 space
+ string "end"
+ many1 space
+ string "mdl"
+ many1 space
+ string "trunc"
+ many1 space
+ string "gc"
+ many1 space
+ string "description"
+ newline
+ string " -"
+ many1 (try (oneOf " -"))
+ newline
+ optional (try (string " ------ inclusion threshold ------"))
+ many newline
+ hitScores' <- many (try genParserCMsearchHit) --`endBy` (try (string "Hit alignments:"))
+ optional (try genParserCMsearchEmptyHit)
+ -- this is followed by hit alignments and internal cmsearch statistics which are not parsed
+ many anyChar
+ eof
+ return $ CMsearch queryCMfile' targetSequenceDatabase' numberOfWorkerThreads' hitScores'
+
+genParserCMsearchHeaderField :: String -> GenParser Char st String
+genParserCMsearchHeaderField fieldname = do
+ string (fieldname ++ ":")
+ many1 space
+ many1 (noneOf "\n")
+ newline
+ return []
+
+genParserCMsearchEmptyHit :: GenParser Char st [CMsearchHit]
+genParserCMsearchEmptyHit = do
+ string " [No hits detected that satisfy reporting thresholds]"
+ newline
+ optional (try newline)
+ return []
+
+genParserCMsearchHit :: GenParser Char st CMsearchHit
+genParserCMsearchHit = do
+ many1 space
+ string "("
+ hitRank' <- many1 digit
+ string ")"
+ many1 space
+ hitSignificant' <- choice [char '!', char '?']
+ many1 space
+ hitEValue' <- many1 (oneOf "0123456789.e-")
+ many1 space
+ hitScore' <- many1 (oneOf "0123456789.e-")
+ many1 space
+ hitBias' <- many1 (oneOf "0123456789.e-")
+ many1 space
+ hitSequenceHeader' <- many1 (noneOf " ")
+ many1 space
+ hitStart' <- many1 digit
+ many1 space
+ hitEnd' <- many1 digit
+ many1 space
+ hitStrand' <- choice [char '+', char '-', char '.']
+ many1 space
+ hitModel' <- many1 letter
+ many1 space
+ hitTruncation' <- many1 (choice [alphaNum, char '\''])
+ many1 space
+ hitGCcontent' <- many1 (oneOf "0123456789.e-")
+ many1 space
+ hitDescription' <- many1 (noneOf "\n")
+ newline
+ optional (try (string " ------ inclusion threshold ------"))
+ optional (try newline)
+ return $ CMsearchHit (readInt hitRank') hitSignificant' (readDouble hitEValue') (readDouble hitScore') (readDouble hitBias') (L.pack hitSequenceHeader') (readInt hitStart') (readInt hitEnd') hitStrand' (L.pack hitModel') (L.pack hitTruncation') (readDouble hitGCcontent') (L.pack hitDescription')
+
+-- | parse from input filePath
+parseCMstat :: String -> Either ParseError CMstat
+parseCMstat = parse genParserCMstat "parseCMstat"
+
+-- | parse from input filePath
+readCMstat :: String -> IO (Either ParseError CMstat)
+readCMstat filePath = do
+ parsedFile <- parseFromFile genParserCMstat filePath
+ CE.evaluate parsedFile
+
+genParserCMstat :: GenParser Char st CMstat
+genParserCMstat = do
+ string "# cmstat :: display summary statistics for CMs"
+ newline
+ string "# INFERNAL "
+ many1 (noneOf "\n")
+ newline
+ string "# Copyright (C) 201"
+ many1 (noneOf "\n")
+ newline
+ string "# Freely distributed under the GNU General Public License (GPLv3)."
+ newline
+ string "# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -"
+ newline
+ char '#'
+ many1 (char ' ')
+ string "rel entropy"
+ newline
+ char '#'
+ many1 (char ' ')
+ many1 (char '-')
+ newline
+ char '#'
+ many1 space
+ string "idx"
+ many1 space
+ string "name"
+ many1 space
+ string "accession"
+ many1 space
+ string "nseq"
+ many1 space
+ string "eff_nseq"
+ many1 space
+ string "clen"
+ many1 space
+ string "W"
+ many1 space
+ string "bps"
+ many1 space
+ string "bifs"
+ many1 space
+ string "model"
+ many1 space
+ string "cm"
+ many1 space
+ string "hmm"
+ newline
+ string "#"
+ many1 (try (oneOf " -"))
+ newline
+ many1 space
+ _statIndex <- many1 digit
+ many1 space
+ _statName <- many1 letter
+ many1 space
+ _statAccession <- many1 (noneOf " ")
+ many1 space
+ _statSequenceNumber <- many1 digit
+ many1 space
+ _statEffectiveSequences <- many1 (oneOf "0123456789.e-")
+ many1 space
+ _statConsensusLength <- many digit
+ many1 space
+ _statW <- many1 digit
+ many1 space
+ _statBasepaires <- many1 digit
+ many1 space
+ _statBifurcations <- many1 digit
+ many1 space
+ _statModel <- many1 letter
+ many1 space
+ _relativeEntropyCM <- many1 (oneOf "0123456789.e-")
+ many1 space
+ _relativeEntropyHMM <- many1 (oneOf "0123456789.e-")
+ newline
+ char '#'
+ newline
+ eof
+ return $ CMstat (readInt _statIndex) _statName _statAccession (readInt _statSequenceNumber) (readDouble _statEffectiveSequences) (readInt _statConsensusLength) (readInt _statW) (readInt _statBasepaires) (readInt _statBifurcations) _statModel (readDouble _relativeEntropyCM) (readDouble _relativeEntropyHMM)
+
+extractCandidateSequences :: [(Sequence,Int,String,Char)] -> V.Vector (Int,Sequence)
+extractCandidateSequences candidates' = indexedSeqences
+ where sequences = map (\(inputSequence,_,_,_) -> inputSequence) candidates'
+ indexedSeqences = V.map (\(number,inputSequence) -> (number + 1,inputSequence))(V.indexed (V.fromList sequences))
+
+extractAlignedSequences :: Int -> ModelConstruction -> V.Vector (Int,Sequence)
+extractAlignedSequences iterationnumber modelconstruction
+ | iterationnumber == 0 = V.map (\(number,seq') -> (number + 1,seq')) (V.indexed (V.fromList [inputSequence]))
+ | otherwise = indexedSeqRecords
+ where inputSequence = inputFasta modelconstruction
+ seqRecordsperTaxrecord = map sequenceRecords (taxRecords modelconstruction)
+ seqRecords = concat seqRecordsperTaxrecord
+ --alignedSeqRecords = filter (\seqRec -> (aligned seqRec) > 0) seqRecords
+ indexedSeqRecords = V.map (\(number,seq') -> (number + 1,seq')) (V.indexed (V.fromList (inputSequence : map nucleotideSequence seqRecords)))
+
+filterByParentTaxId :: [(BlastHit,Int)] -> Bool -> [(BlastHit,Int)]
+filterByParentTaxId blastHitsWithParentTaxId singleHitPerParentTaxId
+ | singleHitPerParentTaxId = singleBlastHitperParentTaxId
+ | otherwise = blastHitsWithParentTaxId
+ where blastHitsWithParentTaxIdSortedByParentTaxId = sortBy compareTaxId blastHitsWithParentTaxId
+ blastHitsWithParentTaxIdGroupedByParentTaxId = groupBy sameTaxId blastHitsWithParentTaxIdSortedByParentTaxId
+ singleBlastHitperParentTaxId = map (maximumBy compareHitEValue) blastHitsWithParentTaxIdGroupedByParentTaxId
+
+filterByHitLength :: [BlastHit] -> Int -> Bool -> [BlastHit]
+filterByHitLength blastHits queryLength filterOn
+ | filterOn = filteredBlastHits
+ | otherwise = blastHits
+ where filteredBlastHits = filter (hitLengthCheck queryLength) blastHits
+
+-- | Hits should have a compareable length to query
+hitLengthCheck :: Int -> BlastHit -> Bool
+hitLengthCheck queryLength blastHit = lengthStatus
+ where blastMatches = matches blastHit
+ minHfrom = minimum (map h_from blastMatches)
+ minHfromHSP = fromJust (find (\hsp -> minHfrom == h_from hsp) blastMatches)
+ maxHto = maximum (map h_to blastMatches)
+ maxHtoHSP = fromJust (find (\hsp -> maxHto == h_to hsp) blastMatches)
+ minHonQuery = q_from minHfromHSP
+ maxHonQuery = q_to maxHtoHSP
+ startCoordinate = minHfrom - minHonQuery
+ endCoordinate = maxHto + (queryLength - maxHonQuery)
+ fullSeqLength = endCoordinate - startCoordinate
+ lengthStatus = fullSeqLength < (queryLength * 3)
+
+filterByCoverage :: [BlastHit] -> Int -> Bool -> [BlastHit]
+filterByCoverage blastHits queryLength filterOn
+ | filterOn = filteredBlastHits
+ | otherwise = blastHits
+ where filteredBlastHits = filter (coverageCheck queryLength) blastHits
+
+-- | Hits should have a compareable length to query
+coverageCheck :: Int -> BlastHit -> Bool
+coverageCheck queryLength blastHit = coverageStatus
+ where blastMatches = matches blastHit
+ maxIdentity = fromIntegral (maximum (map (snd . Bio.BlastXML.identity) blastMatches))
+ coverageStatus = (maxIdentity/(fromIntegral queryLength))* (100 :: Double) >= (80 :: Double)
+
+
+-- | Wrapper for retrieveFullSequence that rerequests incomplete return sequees
+retrieveFullSequences :: StaticOptions -> [(String,Int,Int,String,String,Int,String)] -> IO [(Sequence,Int,String)]
+retrieveFullSequences staticOptions requestedSequences = do
+ fullSequences <- mapM (retrieveFullSequence (tempDirPath staticOptions)) requestedSequences
+ if any (isNothing . firstOfTriple) fullSequences
+ then do
+ let fullSequencesWithRequestedSequences = zip fullSequences requestedSequences
+ --let (failedRetrievals, successfulRetrievals) = partition (\x -> L.null (unSD (seqdata (firstOfTriple (fst x))))) fullSequencesWithRequestedSequences
+ let (failedRetrievals, successfulRetrievals) = partition (isNothing . firstOfTriple . fst) fullSequencesWithRequestedSequences
+ --we try to reretrieve failed entries once
+ missingSequences <- mapM (retrieveFullSequence (tempDirPath staticOptions) .snd) failedRetrievals
+ let (stillMissingSequences,reRetrievedSequences) = partition (isNothing . firstOfTriple) missingSequences
+ logWarning ("Sequence retrieval failed: \n" ++ concatMap show stillMissingSequences ++ "\n") (tempDirPath staticOptions)
+ let unwrappedRetrievals = map (\(x,y,z) -> (fromJust x,y,z)) (map fst successfulRetrievals ++ reRetrievedSequences)
+ CE.evaluate unwrappedRetrievals
+ else CE.evaluate (map (\(x,y,z) -> (fromJust x,y,z)) fullSequences)
+
+retrieveFullSequence :: String -> (String,Int,Int,String,String,Int,String) -> IO (Maybe Sequence,Int,String)
+retrieveFullSequence temporaryDirectoryPath (nucleotideId,seqStart,seqStop,strand,_,taxid,subject') = do
+ let program' = Just "efetch"
+ let database' = Just "nucleotide"
+ let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at"
+ let queryString = "id=" ++ nucleotideId ++ "&seq_start=" ++ show seqStart ++ "&seq_stop=" ++ show seqStop ++ "&rettype=fasta" ++ "&strand=" ++ strand ++ registrationInfo
+ let entrezQuery = EntrezHTTPQuery program' database' queryString
+ result <- CE.catch (entrezHTTP entrezQuery)
+ (\e -> do let err = show (e :: CE.IOException)
+ logWarning ("Warning: Full sequence retrieval failed:" ++ " " ++ err) temporaryDirectoryPath
+ return [])
+ if null result
+ then return (Nothing,taxid,subject')
+ else do
+ if null ((mkSeqs . L.lines) (L.pack result))
+ then return (Nothing,taxid,subject')
+ else do
+ let parsedFasta = head ((mkSeqs . L.lines) (L.pack result))
+ if L.null (unSD (seqdata parsedFasta))
+ then return (Nothing,taxid,subject')
+ else CE.evaluate (Just parsedFasta,taxid,subject')
+
+getRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,String,Int,String)
+getRequestedSequenceElement queryLength (blastHit,taxid)
+ | blastHitIsReverseComplement (blastHit,taxid) = getReverseRequestedSequenceElement queryLength (blastHit,taxid)
+ | otherwise = getForwardRequestedSequenceElement queryLength (blastHit,taxid)
+
+blastHitIsReverseComplement :: (BlastHit,Int) -> Bool
+blastHitIsReverseComplement (blastHit,_) = isReverse
+ where blastMatch = head (matches blastHit)
+ firstHSPfrom = h_from blastMatch
+ firstHSPto = h_to blastMatch
+ isReverse = firstHSPfrom > firstHSPto
+
+getForwardRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,String,Int,String)
+getForwardRequestedSequenceElement queryLength (blastHit,taxid) = (geneIdentifier',startcoordinate,endcoordinate,strand,accession',taxid,subjectBlast)
+ where accession' = L.unpack (extractAccession blastHit)
+ subjectBlast = L.unpack (unSL (subject blastHit))
+ geneIdentifier' = extractGeneId blastHit--
+ blastMatch = head (matches blastHit)
+ blastHitOriginSequenceLength = slength blastHit
+ minHfrom = h_from blastMatch
+ maxHto = h_to blastMatch
+ minHonQuery = q_from blastMatch
+ maxHonQuery = q_to blastMatch
+ --unsafe coordinates may exceed length of available sequence
+ unsafestartcoordinate = minHfrom - minHonQuery
+ unsafeendcoordinate = maxHto + (queryLength - maxHonQuery)
+ startcoordinate = lowerBoundryCoordinateSetter 0 unsafestartcoordinate
+ endcoordinate = upperBoundryCoordinateSetter blastHitOriginSequenceLength unsafeendcoordinate
+ strand = "1"
+ ----
+ --blastMatches = matches blastHit
+ --blastHitOriginSequenceLength = slength blastHit
+ --minHfrom = minimum (map h_from blastMatches)
+ --minHfromHSP = fromJust (find (\hsp -> minHfrom == h_from hsp) blastMatches)
+ --maxHto = maximum (map h_to blastMatches)
+ --maxHtoHSP = fromJust (find (\hsp -> maxHto == h_to hsp) blastMatches)
+ --minHonQuery = q_from minHfromHSP
+ --maxHonQuery = q_to maxHtoHSP
+ --unsafe coordinates may exceed length of available sequence
+ --unsafestartcoordinate = minHfrom - minHonQuery
+ --unsafeendcoordinate = maxHto + (queryLength - maxHonQuery)
+ --startcoordinate = lowerBoundryCoordinateSetter 0 unsafestartcoordinate
+ --endcoordinate = upperBoundryCoordinateSetter blastHitOriginSequenceLength unsafeendcoordinate
+ --strand = "1"
+
+lowerBoundryCoordinateSetter :: Int -> Int -> Int
+lowerBoundryCoordinateSetter lowerBoundry currentValue
+ | currentValue < lowerBoundry = lowerBoundry
+ | otherwise = currentValue
+
+upperBoundryCoordinateSetter :: Int -> Int -> Int
+upperBoundryCoordinateSetter upperBoundry currentValue
+ | currentValue > upperBoundry = upperBoundry
+ | otherwise = currentValue
+
+getReverseRequestedSequenceElement :: Int -> (BlastHit,Int) -> (String,Int,Int,String,String,Int,String)
+getReverseRequestedSequenceElement queryLength (blastHit,taxid) = (geneIdentifier',startcoordinate,endcoordinate,strand,accession',taxid,subjectBlast)
+ where accession' = L.unpack (extractAccession blastHit)
+ subjectBlast = L.unpack (unSL (subject blastHit))
+ geneIdentifier' = extractGeneId blastHit
+ blastMatch = head (matches blastHit)
+ blastHitOriginSequenceLength = slength blastHit
+ maxHfrom = h_from blastMatch
+ minHto = h_to blastMatch
+ minHonQuery = q_from blastMatch
+ maxHonQuery = q_to blastMatch
+ --unsafe coordinates may exceed length of avialable sequence
+ unsafestartcoordinate = maxHfrom + minHonQuery
+ unsafeendcoordinate = minHto - (queryLength - maxHonQuery)
+ startcoordinate = lowerBoundryCoordinateSetter 0 unsafeendcoordinate
+ endcoordinate = upperBoundryCoordinateSetter blastHitOriginSequenceLength unsafestartcoordinate
+ strand = "2"
+ --
+ --blastMatches = matches blastHit
+ --blastHitOriginSequenceLength = slength blastHit
+ --maxHfrom = maximum (map h_from blastMatches)
+ --maxHfromHSP = fromJust (find (\hsp -> maxHfrom == h_from hsp) blastMatches)
+ --minHto = minimum (map h_to blastMatches)
+ --minHtoHSP = fromJust (find (\hsp -> minHto == h_to hsp) blastMatches)
+ --minHonQuery = q_from maxHfromHSP
+ --maxHonQuery = q_to minHtoHSP
+ --unsafe coordinates may exeed length of avialable sequence
+ --unsafestartcoordinate = maxHfrom + minHonQuery
+ --unsafeendcoordinate = minHto - (queryLength - maxHonQuery)
+ --startcoordinate = lowerBoundryCoordinateSetter 0 unsafeendcoordinate
+ --endcoordinate = upperBoundryCoordinateSetter blastHitOriginSequenceLength unsafestartcoordinate
+ --strand = "2"
+
+--computeAlignmentSCIs :: [String] -> [String] -> IO ()
+--computeAlignmentSCIs alignmentFilepaths rnazOutputFilepaths = do
+-- let zippedFilepaths = zip alignmentFilepaths rnazOutputFilepaths
+-- mapM_ systemRNAz zippedFilepaths
+
+alignSequences :: String -> String -> [String] -> [String] -> [String] -> [String] -> IO ()
+alignSequences program' options fastaFilepaths fastaFilepaths2 alignmentFilepaths summaryFilepaths = do
+ let zipped4Filepaths = zip4 fastaFilepaths fastaFilepaths2 alignmentFilepaths summaryFilepaths
+ let zipped3Filepaths = zip3 fastaFilepaths alignmentFilepaths summaryFilepaths
+ let zippedFilepaths = zip fastaFilepaths alignmentFilepaths
+ let timeout = "3600"
+ case program' of
+ "locarna" -> mapM_ (systemlocarna options) zipped4Filepaths
+ "mlocarna" -> mapM_ (systemMlocarna options) zippedFilepaths
+ "mlocarnatimeout" -> mapM_ (systemMlocarnaWithTimeout timeout options) zippedFilepaths
+ "clustalo" -> mapM_ (systemClustalo options) zippedFilepaths
+ _ -> mapM_ (systemClustalw2 options ) zipped3Filepaths
+
+constructFastaFilePaths :: String -> (Int, Sequence) -> String
+constructFastaFilePaths currentDirectory (fastaIdentifier, _) = currentDirectory ++ show fastaIdentifier ++".fa"
+
+constructCMsearchFilePaths :: String -> (Int, Sequence) -> String
+constructCMsearchFilePaths currentDirectory (fastaIdentifier, _) = currentDirectory ++ show fastaIdentifier ++".cmsearch"
+
+-- Smaller e-Values are greater, the maximum function is applied
+compareHitEValue :: (BlastHit,Int) -> (BlastHit,Int) -> Ordering
+compareHitEValue (hit1,_) (hit2,_)
+ | hitEValue hit1 > hitEValue hit2 = LT
+ | hitEValue hit1 < hitEValue hit2 = GT
+ -- in case of equal evalues the first hit is selected
+ | hitEValue hit1 == hitEValue hit2 = GT
+-- comparing (hitEValue . Down . fst)
+compareHitEValue (_,_) (_,_) = EQ
+
+compareTaxId :: (BlastHit,Int) -> (BlastHit,Int) -> Ordering
+compareTaxId (_,taxId1) (_,taxId2)
+ | taxId1 > taxId2 = LT
+ | taxId1 < taxId2 = GT
+ -- in case of equal evalues the first hit is selected
+ | taxId1 == taxId2 = EQ
+compareTaxId (_,_) (_,_) = EQ
+
+sameTaxId :: (BlastHit,Int) -> (BlastHit,Int) -> Bool
+sameTaxId (_,taxId1) (_,taxId2) = taxId1 == taxId2
+
+-- | NCBI uses the e-Value of the best HSP as the Hits e-Value
+hitEValue :: BlastHit -> Double
+hitEValue hit = minimum (map e_val (matches hit))
+
+convertFastaFoldStockholm :: Sequence -> String -> String
+convertFastaFoldStockholm fastasequence foldedStructure = stockholmOutput
+ where alnHeader = "# STOCKHOLM 1.0\n\n"
+ --(L.unpack (unSL (seqheader inputFasta')))) ++ "\n" ++ (map toUpper (L.unpack (unSD (seqdata inputFasta')))) ++ "\n"
+ seqIdentifier = L.unpack (unSL (seqheader fastasequence))
+ seqSequence = L.unpack (unSD (seqdata fastasequence))
+ identifierLength = length seqIdentifier
+ spacerLength' = maximum [14,identifierLength + 2]
+ spacer = replicate (spacerLength' - identifierLength) ' '
+ entrystring = seqIdentifier ++ spacer ++ seqSequence ++ "\n"
+ structureString = "#=GC SS_cons" ++ replicate (spacerLength' - 12) ' ' ++ foldedStructure ++ "\n"
+ bottom = "//"
+ stockholmOutput = alnHeader ++ entrystring ++ structureString ++ bottom
+
+convertClustaltoStockholm :: StructuralClustalAlignment -> String
+convertClustaltoStockholm parsedMlocarnaAlignment = stockholmOutput
+ where alnHeader = "# STOCKHOLM 1.0\n\n"
+ clustalAlignment = structuralAlignmentEntries parsedMlocarnaAlignment
+ uniqueIds = nub (map entrySequenceIdentifier clustalAlignment)
+ mergedEntries = map (mergeEntry clustalAlignment) uniqueIds
+ maxIdentifierLenght = maximum (map (length . entrySequenceIdentifier) clustalAlignment)
+ spacerLength' = maxIdentifierLenght + 2
+ stockholmEntries = concatMap (buildStockholmAlignmentEntries spacerLength') mergedEntries
+ structureString = "#=GC SS_cons" ++ replicate (spacerLength' - 12) ' ' ++ secondaryStructureTrack parsedMlocarnaAlignment ++ "\n"
+ bottom = "//"
+ stockholmOutput = alnHeader ++ stockholmEntries ++ structureString ++ bottom
+
+mergeEntry :: [ClustalAlignmentEntry] -> String -> ClustalAlignmentEntry
+mergeEntry clustalAlignment uniqueId = mergedEntry
+ where idEntries = filter (\entry -> entrySequenceIdentifier entry==uniqueId) clustalAlignment
+ mergedSeq = foldr ((++) . entryAlignedSequence) "" idEntries
+ mergedEntry = ClustalAlignmentEntry uniqueId mergedSeq
+
+buildStockholmAlignmentEntries :: Int -> ClustalAlignmentEntry -> String
+buildStockholmAlignmentEntries inputSpacerLength entry = entrystring
+ where idLength = length (filter (/= '\n') (entrySequenceIdentifier entry))
+ spacer = replicate (inputSpacerLength - idLength) ' '
+ entrystring = entrySequenceIdentifier entry ++ spacer ++ entryAlignedSequence entry ++ "\n"
+
+retrieveTaxonomicContextEntrez :: Int -> IO (Maybe Taxon)
+retrieveTaxonomicContextEntrez inputTaxId = do
+ let program' = Just "efetch"
+ let database' = Just "taxonomy"
+ let taxIdString = show inputTaxId
+ let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at"
+ let queryString = "id=" ++ taxIdString ++ registrationInfo
+ let entrezQuery = EntrezHTTPQuery program' database' queryString
+ result <- entrezHTTP entrezQuery
+ if null result
+ then do
+ error "Could not retrieve taxonomic context from NCBI Entrez, cannot proceed."
+ return Nothing
+ else do
+ let taxon = head (readEntrezTaxonSet result)
+ --print taxon
+ if null (lineageEx taxon)
+ then error "Retrieved taxonomic context taxon from NCBI Entrez with empty lineage, cannot proceed."
+ else CE.evaluate (Just taxon)
+
+retrieveParentTaxIdEntrez :: [(BlastHit,Int)] -> IO [(BlastHit,Int)]
+retrieveParentTaxIdEntrez blastHitsWithHitTaxids = do
+ if not (null blastHitsWithHitTaxids)
+ then do
+ let program' = Just "efetch"
+ let database' = Just "taxonomy"
+ let extractedBlastHits = map fst blastHitsWithHitTaxids
+ let taxIds = map snd blastHitsWithHitTaxids
+ let taxIdStrings = map show taxIds
+ let taxIdQuery = intercalate "," taxIdStrings
+ let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at"
+ let queryString = "id=" ++ taxIdQuery ++ registrationInfo
+ let entrezQuery = EntrezHTTPQuery program' database' queryString
+ result <- entrezHTTP entrezQuery
+ let parentTaxIds = readEntrezParentIds result
+ if null parentTaxIds
+ then return []
+ else CE.evaluate (zip extractedBlastHits parentTaxIds)
+ else return []
+
+-- | Wrapper functions that ensures that only 20 queries are sent per request
+retrieveParentTaxIdsEntrez :: [(BlastHit,Int)] -> IO [(BlastHit,Int)]
+retrieveParentTaxIdsEntrez taxIdwithBlastHits = do
+ let splits = portionListElements taxIdwithBlastHits 20
+ taxIdsOutput <- mapM retrieveParentTaxIdEntrez splits
+ return (concat taxIdsOutput)
+
+-- | Wrapper functions that ensures that only 20 queries are sent per request
+retrieveBlastHitsTaxIdEntrez :: [BlastHit] -> IO [([BlastHit],String)]
+retrieveBlastHitsTaxIdEntrez blastHits = do
+ let splits = portionListElements blastHits 20
+ mapM retrieveBlastHitTaxIdEntrez splits
+
+
+retrieveBlastHitTaxIdEntrez :: [BlastHit] -> IO ([BlastHit],String)
+retrieveBlastHitTaxIdEntrez blastHits = do
+ if not (null blastHits)
+ then do
+ let geneIds = map extractGeneId blastHits
+ let idList = intercalate "," geneIds
+ let registrationInfo = buildRegistration "RNAlien" "florian.eggenhofer@univie.ac.at"
+ let query' = "id=" ++ idList ++ registrationInfo
+ let entrezQuery = EntrezHTTPQuery (Just "esummary") (Just "nucleotide") query'
+ threadDelay 10000000
+ result <- entrezHTTP entrezQuery
+ CE.evaluate (blastHits,result)
+ else return (blastHits,"")
+
+extractTaxIdFromEntrySummaries :: String -> [Int]
+extractTaxIdFromEntrySummaries input
+ | null input = []
+ | null parsedResultList = []
+ | otherwise = hitTaxIds
+ where parsedResultList = readEntrezSummaries input
+ parsedResult = head parsedResultList
+ blastHitSummaries = documentSummaries parsedResult
+ hitTaxIdStrings = map extractTaxIdfromDocumentSummary blastHitSummaries
+ hitTaxIds = map readInt hitTaxIdStrings
+
+extractAccession :: BlastHit -> L.ByteString
+extractAccession currentBlastHit = accession'
+ where splitedFields = DS.splitOn "|" (L.unpack (hitId currentBlastHit))
+ accession' = L.pack (splitedFields !! 3)
+
+extractGeneId :: BlastHit -> String
+extractGeneId currentBlastHit = nucleotideId
+ where truncatedId = drop 3 (L.unpack (hitId currentBlastHit))
+ pipeSymbolIndex = fromJust (elemIndex '|' truncatedId)
+ nucleotideId = take pipeSymbolIndex truncatedId
+
+extractTaxIdfromDocumentSummary :: EntrezDocSum -> String
+extractTaxIdfromDocumentSummary documentSummary = itemContent (fromJust (find (\item -> "TaxId" == itemName item) (summaryItems documentSummary)))
+
+getBestHit :: BlastResult -> BlastHit
+getBestHit blastResult
+ | null (concatMap hits (results blastResult)) = error "getBestHit - head: empty list"
+ | otherwise = head (hits (head (results blastResult)))
+
+-- Blast returns low evalues with zero instead of the exact number
+getHitWithFractionEvalue :: BlastResult -> Maybe BlastMatch
+getHitWithFractionEvalue blastResult
+ | null (concatMap hits (results blastResult)) = Nothing
+ | otherwise = find (\match -> e_val match /= (0 ::Double)) (concatMap matches (concatMap hits (results blastResult)))
+
+showlines :: Show a => [a] -> String
+showlines = concatMap (\x -> show x ++ "\n")
+
+logMessage :: String -> String -> IO ()
+logMessage logoutput temporaryDirectoryPath = appendFile (temporaryDirectoryPath ++ "Log") logoutput
+
+logWarning :: String -> String -> IO ()
+logWarning logoutput temporaryDirectoryPath = appendFile (temporaryDirectoryPath ++ "log/warnings") logoutput
+
+logVerboseMessage :: Bool -> String -> String -> IO ()
+logVerboseMessage verboseTrue logoutput temporaryDirectoryPath
+ | verboseTrue = appendFile (temporaryDirectoryPath ++ "Log") (show logoutput)
+ | otherwise = return ()
+
+logEither :: (Show a) => Either a b -> String -> IO ()
+logEither (Left logoutput) temporaryDirectoryPath = appendFile (temporaryDirectoryPath ++ "Log") (show logoutput)
+logEither _ _ = return ()
+
+checkTools :: [String] -> String -> IO (Either String String)
+checkTools tools temporaryDirectoryPath = do
+ -- check if all tools are available via PATH or Left
+ checks <- mapM checkTool tools
+ if not (null (lefts checks))
+ then return (Left (concat (lefts checks)))
+ else do
+ logMessage ("Tools : " ++ intercalate "," tools ++ "\n") temporaryDirectoryPath
+ return (Right "Tools ok")
+
+logToolVersions :: String -> IO ()
+logToolVersions temporaryDirectoryPath = do
+ let clustaloversionpath = temporaryDirectoryPath ++ "log/clustalo.version"
+ let mlocarnaversionpath = temporaryDirectoryPath ++ "log/mlocarna.version"
+ let rnafoldversionpath = temporaryDirectoryPath ++ "log/RNAfold.version"
+ let infernalversionpath = temporaryDirectoryPath ++ "log/Infernal.version"
+ _ <- system ("clustalo --version >" ++ clustaloversionpath)
+ _ <- system ("mlocarna --version >" ++ mlocarnaversionpath)
+ _ <- system ("RNAfold --version >" ++ rnafoldversionpath)
+ _ <- system ("cmcalibrate -h >" ++ infernalversionpath)
+ -- _ <- system ("RNAz" ++ rnazversionpath)
+ -- _ <- system ("CMCompare >" ++ infernalversionpath)
+ clustaloversion <- readFile clustaloversionpath
+ mlocarnaversion <- readFile mlocarnaversionpath
+ rnafoldversion <- readFile rnafoldversionpath
+ infernalversionOutput <- readFile infernalversionpath
+ let infernalversion = lines infernalversionOutput !! 1
+ let messageString = "Clustalo version: " ++ clustaloversion ++ "mlocarna version: " ++ mlocarnaversion ++ "RNAfold version: " ++ rnafoldversion ++ "infernalversion: " ++ infernalversion ++ "\n"
+ logMessage messageString temporaryDirectoryPath
+
+checkTool :: String -> IO (Either String String)
+checkTool tool = do
+ toolcheck <- findExecutable tool
+ if isJust toolcheck
+ then return (Right (fromJust toolcheck))
+ else return (Left ("RNAlien could not find "++ tool ++ " in your $PATH and has to abort.\n"))
+
+constructTaxonomyRecordsCSVTable :: ModelConstruction -> String
+constructTaxonomyRecordsCSVTable modelconstruction = csvtable
+ where tableheader = "Taxonomy Id;Added in Iteration Step;Entry Header\n"
+ tablebody = concatMap constructTaxonomyRecordCSVEntries (taxRecords modelconstruction)
+ csvtable = tableheader ++ tablebody
+
+constructTaxonomyRecordCSVEntries :: TaxonomyRecord -> String
+constructTaxonomyRecordCSVEntries taxRecord = concatMap (\seqrec -> show (recordTaxonomyId taxRecord) ++ ";" ++ show (aligned seqrec) ++ ";" ++ filter (/= ';') (L.unpack (unSL (seqheader (nucleotideSequence seqrec)))) ++ "\n") (sequenceRecords taxRecord)
+
+setVerbose :: Verbosity -> Bool
+setVerbose verbosityLevel
+ | verbosityLevel == Loud = True
+ | otherwise = False
+
+evaluateConstructionResult :: StaticOptions -> Int -> IO String
+evaluateConstructionResult staticOptions entryNumber = do
+ let evaluationDirectoryFilepath = tempDirPath staticOptions ++ "evaluation/"
+ createDirectoryIfMissing False evaluationDirectoryFilepath
+ let fastaFilepath = tempDirPath staticOptions ++ "result.fa"
+ let clustalFilepath = evaluationDirectoryFilepath ++ "result.clustal"
+ let reformatedClustalPath = evaluationDirectoryFilepath ++ "result.clustal.reformated"
+ let cmFilepath = tempDirPath staticOptions ++ "result.cm"
+ systemCMalign ("--outformat=Clustal --cpu " ++ show (cpuThreads staticOptions)) cmFilepath fastaFilepath clustalFilepath
+ let resultModelStatistics = tempDirPath staticOptions ++ "result.cmstat"
+ systemCMstat cmFilepath resultModelStatistics
+ inputcmStat <- readCMstat resultModelStatistics
+ let cmstatString = cmstatEvalOutput inputcmStat
+ if entryNumber > 1
+ then do
+ let resultRNAz = tempDirPath staticOptions ++ "result.rnaz"
+ rnazClustalpath <- preprocessClustalForRNAzExternal clustalFilepath reformatedClustalPath
+ if isRight rnazClustalpath
+ then do
+ systemRNAz "-l" (fromRight rnazClustalpath) resultRNAz
+ inputRNAz <- readRNAz resultRNAz
+ let rnaZString = rnaZEvalOutput inputRNAz
+ return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAz statistics for result alignment: " ++ rnaZString)
+ else do
+ logWarning ("Running RNAz for result evalution encountered a problem:" ++ fromLeft rnazClustalpath) (tempDirPath staticOptions)
+ return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAz statistics for result alignment: Running RNAz for result evalution encountered a problem\n" ++ fromLeft rnazClustalpath)
+ else do
+ logWarning "Message: RNAlien could not find additional covariant sequences\n Could not run RNAz statistics. Could not run RNAz statistics with a single sequence.\n" (tempDirPath staticOptions)
+ return ("\nEvaluation of RNAlien result :\nCMstat statistics for result.cm\n" ++ cmstatString ++ "\nRNAlien could not find additional covariant sequences. Could not run RNAz statistics with a single sequence.\n")
+
+
+cmstatEvalOutput :: Either ParseError CMstat -> String
+cmstatEvalOutput inputcmstat
+ | isRight inputcmstat = cmstatString
+ | otherwise = show (fromLeft inputcmstat)
+ where cmStat = fromRight inputcmstat
+ cmstatString = " Sequence Number: " ++ show (statSequenceNumber cmStat)++ "\n" ++ " Effective Sequences: " ++ show (statEffectiveSequences cmStat)++ "\n" ++ " Consensus length: " ++ show (statConsensusLength cmStat) ++ "\n" ++ " Expected maximum hit-length: " ++ show (statW cmStat) ++ "\n" ++ " Basepairs: " ++ show (statBasepairs cmStat)++ "\n" ++ " Bifurcations: " ++ show (statBifurcations cmStat) ++ "\n" ++ " Modeltype: " ++ show (statModel cmStat) ++ "\n" ++ " Relative Entropy CM: " ++ show (relativeEntropyCM cmStat) ++ "\n" ++ " Relative Entropy HMM: " ++ show (relativeEntropyHMM cmStat) ++ "\n"
+
+rnaZEvalOutput :: Either ParseError RNAz -> String
+rnaZEvalOutput inputRNAz
+ | isRight inputRNAz = rnazString
+ | otherwise = show (fromLeft inputRNAz)
+ where rnaZ = fromRight inputRNAz
+ rnazString = " Mean pairwise identity: " ++ show (meanPairwiseIdentity rnaZ) ++ "\n Shannon entropy: " ++ show (shannonEntropy rnaZ) ++ "\n GC content: " ++ show (gcContent rnaZ) ++ "\n Mean single sequence minimum free energy: " ++ show (meanSingleSequenceMinimumFreeEnergy rnaZ) ++ "\n Consensus minimum free energy: " ++ show (consensusMinimumFreeEnergy rnaZ) ++ "\n Energy contribution: " ++ show (energyContribution rnaZ) ++ "\n Covariance contribution: " ++ show (covarianceContribution rnaZ) ++ "\n Combinations pair: " ++ show (combinationsPair rnaZ) ++ "\n Mean z-score: " ++ show (meanZScore rnaZ) ++ "\n Structure conservation index: " ++ show (structureConservationIndex rnaZ) ++ "\n Background model: " ++ backgroundModel rnaZ ++ "\n Decision model: " ++ decisionModel rnaZ ++ "\n SVM decision value: " ++ show (svmDecisionValue rnaZ) ++ "\n SVM class propability: " ++ show (svmRNAClassProbability rnaZ) ++ "\n Prediction: " ++ prediction rnaZ
+
+-- | Call for external preprocessClustalForRNAz
+preprocessClustalForRNAzExternal :: String -> String -> IO (Either String String)
+preprocessClustalForRNAzExternal clustalFilepath reformatedClustalPath = do
+ clustalString <- readFile clustalFilepath
+ --change clustal format for rnazSelectSeqs.pl
+ let reformatedClustalString = map reformatAln clustalString
+ writeFile reformatedClustalPath reformatedClustalString
+ --select representative entries from result.Clustal with select_sequences
+ let selectedClustalpath = clustalFilepath ++ ".selected"
+ system ("rnazSelectSeqs.pl " ++ reformatedClustalPath ++ " >" ++ selectedClustalpath)
+ return (Right selectedClustalpath)
+
+-- | RNAz can process 500 sequences at max. Using rnazSelectSeqs to isolate representative sample. rnazSelectSeqs only accepts - gap characters, alignment is reformatted accordingly.
+preprocessClustalForRNAz :: String -> String -> IO (Either String String)
+preprocessClustalForRNAz clustalFilepath reformatedClustalPath = do
+ clustalString <- readFile clustalFilepath
+ if length (lines clustalString) > 500
+ then do
+ --change clustal format for rnazSelectSeqs.pl
+ let reformatedClustalString = map reformatAln clustalString
+ writeFile reformatedClustalPath reformatedClustalString
+ --select representative entries from result.Clustal with select_sequences
+ let selectedClustalpath = clustalFilepath ++ ".selected"
+ parsedClustalInput <- readClustalAlignment clustalFilepath
+ if isRight parsedClustalInput
+ then do
+ let filteredClustalInput = rnaZSelectSeqs (fromRight parsedClustalInput) 500 99
+ writeFile selectedClustalpath (show filteredClustalInput)
+ return (Right selectedClustalpath)
+ else return (Left (show (fromLeft parsedClustalInput)))
+ else return (Right clustalFilepath)
+
+-- Iteratively removes sequences with decreasing similarity until target number of alignment entries is reached.
+rnaZSelectSeqs :: ClustalAlignment -> Int -> Double -> ClustalAlignment
+rnaZSelectSeqs currentClustalAlignment targetEntries identityCutoff
+ | targetEntries < numberOfEntries = rnaZSelectSeqs filteredAlignment targetEntries (identityCutoff - 1)
+ | otherwise = currentClustalAlignment
+ where numberOfEntries = length (alignmentEntries currentClustalAlignment)
+ filteredEntries = filterIdenticalAlignmentEntry (alignmentEntries currentClustalAlignment) identityCutoff
+ filteredAlignment = ClustalAlignment filteredEntries (conservationTrack currentClustalAlignment)
+
+reformatAln :: Char -> Char
+reformatAln c
+ | c == '.' = '-'
+ | c == '~' = '-'
+ | c == '_' = '-'
+ | c == 'u' = 'U'
+ | c == 't' = 'T'
+ | c == 'g' = 'G'
+ | c == 'c' = 'C'
+ | c == 'a' = 'A'
+ | otherwise = c
+
+-- | Check if alien can connect to NCBI
+checkNCBIConnection :: IO (Either String String)
+checkNCBIConnection = do
+ response <- simpleHTTP (getRequest "http://www.ncbi.nlm.nih.gov")
+ if isRight response
+ then do
+ let rightResponse = fromRight response
+ if rspCode rightResponse == (2,0,0)
+ then return (Right ("Network connection with NCBI server is ok: " ++ show (rspCode rightResponse)))
+ else return (Left ("Could not connect to NCBI server \"http://www.ncbi.nlm.nih.gov\". Response Code: " ++ show (rspCode rightResponse)))
+ else return (Left ("Could not connect to NCBI server: \"http://www.ncbi.nlm.nih.gov\": " ++ show (fromLeft response)))
+
+-- | Blast evalue is set stricter in inital alignment mode
+setBlastExpectThreshold :: ModelConstruction -> Double
+setBlastExpectThreshold modelConstruction
+ | alignmentModeInfernal modelConstruction = 1 :: Double
+ | otherwise = 0.1 :: Double
+
+reformatFasta :: Sequence -> Sequence
+reformatFasta input = Seq (seqheader input) updatedSequence Nothing
+ where updatedSequence = SeqData (L.pack (map reformatFastaSequence (L.unpack (unSD (seqdata input)))))
+
+reformatFastaSequence :: Char -> Char
+reformatFastaSequence c
+ | c == '.' = '-'
+ | c == '~' = '-'
+ | c == '_' = '-'
+ | c == 'u' = 'T'
+ | c == 't' = 'T'
+ | c == 'g' = 'G'
+ | c == 'c' = 'C'
+ | c == 'a' = 'A'
+ | c == 'U' = 'T'
+ | otherwise = c
+
+setRestrictedTaxonomyLimits :: String -> (Maybe Int,Maybe Int)
+setRestrictedTaxonomyLimits trestriction
+ | trestriction == "bacteria" = (Just (2 :: Int), Nothing)
+ | trestriction == "archea" = (Just (2157 :: Int), Nothing)
+ | trestriction == "eukaryia" = (Just (2759 :: Int), Nothing)
+ | otherwise = (Nothing, Nothing)
+
+checkTaxonomyRestriction :: Maybe String -> Maybe String
+checkTaxonomyRestriction taxonomyRestriction
+ | isJust taxonomyRestriction = checkTaxonomyRestrictionString (fromJust taxonomyRestriction)
+ | otherwise = Nothing
+
+checkTaxonomyRestrictionString :: String -> Maybe String
+checkTaxonomyRestrictionString restrictionString
+ | restrictionString == "archea" = Just "archea"
+ | restrictionString == "bacteria" = Just "bacteria"
+ | restrictionString == "eukaryia" = Just "eukaryia"
+ | otherwise = Nothing
diff --git a/src/Bio/RNAlienStatistics.hs b/src/Bio/RNAlienStatistics.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/RNAlienStatistics.hs
@@ -0,0 +1,256 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+
+-- | Statistics for RNAlien Results
+-- dist/build/RNAlienStatistics/RNAlienStatistics -i /scratch/egg/temp/cm13676/1/model.cm -r /home/mescalin/egg/current/Data/AlienTest/cms/BsrG.cm -g /scratch/egg/temp/AlienSearch/genomes/ -o /scratch/egg/temp/AlienStatistics
+module Main where
+
+import System.Console.CmdArgs
+import Data.Either.Unwrap
+import System.Process
+import qualified Data.ByteString.Lazy.Char8 as L
+import Bio.RNAlienLibrary
+import System.Directory
+import Bio.Core.Sequence
+import Bio.Sequence.Fasta
+import Data.List
+import qualified System.FilePath as FP
+import qualified Data.List.Split as DS
+import Text.Printf
+import Bio.RNAzParser
+
+data Options = Options
+ { alienCovarianceModelPath :: String,
+ alienrnazPath :: String,
+ aliencmstatPath :: String,
+ rfamCovarianceModelPath :: String,
+ rfamFastaFilePath :: String,
+ alienFastaFilePath :: String,
+ rfamModelName :: String,
+ rfamModelId :: String,
+ rfamThreshold :: Double,
+ alienThreshold :: Double,
+ outputDirectoryPath :: String,
+ benchmarkIndex :: Int,
+ threads :: Int
+ } deriving (Show,Data,Typeable)
+
+options :: Options
+options = Options
+ { alienCovarianceModelPath = def &= name "i" &= help "Path to alienCovarianceModelPath",
+ alienrnazPath = def &= name "z" &= help "Path to alienRNAzResult",
+ aliencmstatPath = def &= name "m" &= help "Path to aliencmstatResult",
+ rfamCovarianceModelPath = def &= name "r" &= help "Path to rfamCovarianceModelPath",
+ rfamFastaFilePath = def &= name "g" &= help "Path to rfamFastaFile",
+ rfamModelName = def &= name "n" &= help "Rfam model name",
+ rfamModelId = def &= name "d" &= help "Rfam model id",
+ alienFastaFilePath = def &= name "a" &= help "Path to alienFastaFile",
+ outputDirectoryPath = def &= name "o" &= help "Path to output directory",
+ alienThreshold = 20 &= name "t" &= help "Bitscore threshold for RNAlien model hits on Rfam fasta, default 20",
+ rfamThreshold = 20 &= name "x" &= help "Bitscore threshold for Rfam model hits on Alien fasta, default 20",
+ benchmarkIndex = 1 &= name "b" &= help "Index used to identify sRNA tagged RNA families",
+ threads = 1 &= name "c" &= help "Number of available cpu slots/cores, default 1"
+ } &= summary "RNAlienStatistics devel version" &= help "Florian Eggenhofer - >2013" &= verbosity
+
+--cmSearchFasta threads rfamCovarianceModelPath outputDirectoryPath "Rfam" False genomesDirectoryPath
+cmSearchFasta :: Int -> Double -> Int -> String -> String -> String -> String -> IO [CMsearchHit]
+cmSearchFasta benchmarkIndex thresholdScore cpuThreads covarianceModelPath outputDirectory modelType fastapath = do
+ createDirectoryIfMissing False (outputDirectory ++ "/" ++ modelType)
+ _ <- systemCMsearch cpuThreads "" covarianceModelPath fastapath (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch")
+ result <- readCMSearch (outputDirectory ++ "/" ++ modelType ++ "/" ++ (show benchmarkIndex) ++ ".cmsearch")
+ if (isLeft result)
+ then do
+ print (fromLeft result)
+ return []
+ else do
+ let rightResults = fromRight result
+ let (significantHits,_) = partitionCMsearchHitsByScore thresholdScore rightResults
+ let organismUniquesignificantHits = nubBy cmSearchSameOrganism significantHits
+ return organismUniquesignificantHits
+
+partitionCMsearchHitsByScore :: Double -> CMsearch -> ([CMsearchHit],[CMsearchHit])
+partitionCMsearchHitsByScore thresholdScore cmSearchResult = (selected,rejected)
+ where (selected,rejected) = partition (\hit -> hitScore hit >= thresholdScore) (cmsearchHits cmSearchResult)
+
+trimCMsearchFastaFile :: String -> String -> String -> CMsearch -> String -> IO ()
+trimCMsearchFastaFile genomesDirectory outputFolder modelType cmsearch fastafile = do
+ let fastaInputPath = genomesDirectory ++ "/" ++ fastafile
+ let fastaOutputPath = outputFolder ++ "/" ++ modelType ++ "/" ++ fastafile
+ fastaSequences <- readFasta fastaInputPath
+ let trimmedSequence = trimCMsearchSequence cmsearch (head fastaSequences)
+ writeFasta fastaOutputPath [trimmedSequence]
+
+trimCMsearchSequence :: CMsearch -> Sequence -> Sequence
+trimCMsearchSequence cmSearchResult inputSequence = subSequence
+ where hitScoreEntry = head (cmsearchHits cmSearchResult)
+ sequenceString = L.unpack (unSD (seqdata inputSequence))
+ sequenceSubstring = cmSearchsubString (hitStart hitScoreEntry) (hitEnd hitScoreEntry) sequenceString
+ newSequenceHeader = L.pack ((L.unpack (unSL (seqheader inputSequence))) ++ "cmS_" ++ (show (hitStart hitScoreEntry)) ++ "_" ++ (show (hitEnd hitScoreEntry)) ++ "_" ++ (show (hitStrand hitScoreEntry)))
+ subSequence = Seq (SeqLabel newSequenceHeader) (SeqData (L.pack sequenceSubstring)) Nothing
+
+cmSearchSameOrganism :: CMsearchHit -> CMsearchHit -> Bool
+cmSearchSameOrganism hitscore1 hitscore2
+ | hitOrganism1 == hitOrganism2 = True
+ | otherwise = False
+ where unpackedSeqHeader1 = (L.unpack (hitSequenceHeader hitscore1))
+ unpackedSeqHeader2 = (L.unpack (hitSequenceHeader hitscore2))
+ separationcharacter1 = selectSeparationChar unpackedSeqHeader1
+ separationcharacter2 = selectSeparationChar unpackedSeqHeader2
+ hitOrganism1 = (DS.splitOn separationcharacter1 unpackedSeqHeader1) !! 0
+ hitOrganism2 = (DS.splitOn separationcharacter2 unpackedSeqHeader2) !! 0
+
+selectSeparationChar :: String -> String
+selectSeparationChar inputString
+ | any (\a -> a == ':') inputString = ":"
+ | otherwise = "/"
+
+main :: IO ()
+main = do
+ Options{..} <- cmdArgs options
+ rfamModelExists <- doesFileExist rfamCovarianceModelPath
+ verbose <- getVerbosity
+ rnazString <- rnazOutput verbose alienrnazPath
+ cmStatString <- cmStatOutput verbose aliencmstatPath
+ if rfamModelExists
+ then do
+ --compute linkscore
+ linkscore <- compareCM rfamCovarianceModelPath alienCovarianceModelPath outputDirectoryPath
+ rfamMaxLinkScore <- compareCM rfamCovarianceModelPath rfamCovarianceModelPath outputDirectoryPath
+ alienMaxLinkscore <- compareCM alienCovarianceModelPath alienCovarianceModelPath outputDirectoryPath
+ _ <- system ("cat " ++ rfamFastaFilePath ++ " | grep '>' | wc -l >" ++ outputDirectoryPath ++ FP.takeFileName rfamFastaFilePath ++ ".entries")
+ _ <- system ("cat " ++ alienFastaFilePath ++ " | grep '>' | wc -l >" ++ outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries")
+ rfamFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName rfamFastaFilePath ++ ".entries")
+ alienFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries")
+ let rfamFastaEntriesNumber = read rfamFastaEntries :: Int
+ let alienFastaEntriesNumber = read alienFastaEntries :: Int
+ rfamonAlienResults <- cmSearchFasta benchmarkIndex rfamThreshold threads rfamCovarianceModelPath outputDirectoryPath "rfamOnAlien" alienFastaFilePath
+ alienonRfamResults <- cmSearchFasta benchmarkIndex alienThreshold threads alienCovarianceModelPath outputDirectoryPath "alienOnRfam" rfamFastaFilePath
+ let rfamonAlienResultsNumber = length rfamonAlienResults
+ let alienonRfamResultsNumber = length alienonRfamResults
+ let rfamonAlienRecovery = (fromIntegral rfamonAlienResultsNumber :: Double) / (fromIntegral alienFastaEntriesNumber :: Double)
+ let alienonRfamRecovery = (fromIntegral alienonRfamResultsNumber :: Double) / (fromIntegral rfamFastaEntriesNumber :: Double)
+ if (verbose == Loud)
+ then do
+ putStrLn ("BenchmarkIndex: " ++ show benchmarkIndex)
+ putStrLn ("RfamModelName: " ++ rfamModelName)
+ putStrLn ("RfamModelId: " ++ rfamModelId)
+ putStrLn ("Linkscore: " ++ show linkscore)
+ putStrLn ("rfamMaxLinkScore: " ++ show rfamMaxLinkScore)
+ putStrLn ("alienMaxLinkscore: " ++ show alienMaxLinkscore)
+ putStrLn ("rfamGatheringThreshold: " ++ show rfamThreshold)
+ putStrLn ("alienGatheringThreshold: " ++ show alienThreshold)
+ putStrLn ("rfamFastaEntriesNumber: " ++ show rfamFastaEntriesNumber)
+ putStrLn ("alienFastaEntriesNumber: " ++ show alienFastaEntriesNumber)
+ putStrLn ("rfamonAlienResultsNumber: " ++ show rfamonAlienResultsNumber)
+ putStrLn ("alienonRfamResultsNumber: " ++ show alienonRfamResultsNumber)
+ putStrLn ("RfamonAlienRecovery: " ++ show rfamonAlienRecovery)
+ putStrLn ("AlienonRfamRecovery: " ++ show alienonRfamRecovery)
+ print rnazString
+ print cmStatString
+ else do
+ putStrLn (show benchmarkIndex ++ "\t" ++ rfamModelName ++ "\t" ++ rfamModelId ++ "\t" ++ show linkscore ++ "\t" ++ show rfamMaxLinkScore ++ "\t" ++ show alienMaxLinkscore ++ "\t" ++ show rfamThreshold ++ "\t" ++ show alienThreshold ++ "\t" ++ show rfamFastaEntriesNumber ++ "\t" ++ show alienFastaEntriesNumber ++ "\t" ++ show rfamonAlienResultsNumber ++ "\t" ++ show alienonRfamResultsNumber ++ "\t" ++ printf "%.2f" rfamonAlienRecovery ++ "\t" ++ printf "%.2f" alienonRfamRecovery ++ "\t" ++ rnazString ++ "\t" ++ cmStatString)
+ else do
+ --compute linkscore
+ --linkscore <- compareCM rfamCovarianceModelPath alienCovarianceModelPath outputDirectoryPath
+ --rfamMaxLinkScore <- compareCM rfamCovarianceModelPath rfamCovarianceModelPath outputDirectoryPath
+ alienMaxLinkscore <- compareCM alienCovarianceModelPath alienCovarianceModelPath outputDirectoryPath
+ --_ <- system ("cat " ++ rfamFastaFilePath ++ " | grep '>' | wc -l >" ++ outputDirectoryPath ++ FP.takeFileName rfamFastaFilePath ++ ".entries")
+ _ <- system ("cat " ++ alienFastaFilePath ++ " | grep '>' | wc -l >" ++ outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries")
+ --rfamFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName rfamFastaFilePath ++ ".entries")
+ alienFastaEntries <- readFile (outputDirectoryPath ++ FP.takeFileName alienFastaFilePath ++ ".entries")
+ --let rfamFastaEntriesNumber = read rfamFastaEntries :: Int
+ let alienFastaEntriesNumber = read alienFastaEntries :: Int
+ --rfamonAlienResults <- cmSearchFasta benchmarkIndex rfamThreshold threads rfamCovarianceModelPath outputDirectoryPath "rfamOnAlien" alienFastaFilePath
+ --alienonRfamResults <- cmSearchFasta benchmarkIndex alienThreshold threads alienCovarianceModelPath outputDirectoryPath "alienOnRfam" rfamFastaFilePath
+ --let rfamonAlienResultsNumber = length rfamonAlienResults
+ --let alienonRfamResultsNumber = length alienonRfamResults
+ --let rfamonAlienRecovery = (fromIntegral rfamonAlienResultsNumber :: Double) / (fromIntegral alienFastaEntriesNumber :: Double)
+ --let alienonRfamRecovery = (fromIntegral alienonRfamResultsNumber :: Double) / (fromIntegral rfamFastaEntriesNumber :: Double)
+ if (verbose == Loud)
+ then do
+ putStrLn ("BenchmarkIndex:")
+ putStrLn ("RfamModelName: -")
+ putStrLn ("RfamModelId: -")
+ putStrLn ("Linkscore: -")
+ putStrLn ("rfamMaxLinkScore: -")
+ putStrLn ("alienMaxLinkscore: " ++ show alienMaxLinkscore)
+ putStrLn ("rfamGatheringThreshold: -")
+ putStrLn ("alienGatheringThreshold: -")
+ putStrLn ("rfamFastaEntriesNumber: -")
+ putStrLn ("alienFastaEntriesNumber: " ++ show alienFastaEntriesNumber)
+ putStrLn ("rfamonAlienResultsNumber: -")
+ putStrLn ("alienonRfamResultsNumber: -")
+ putStrLn ("RfamonAlienRecovery: -")
+ putStrLn ("AlienonRfamRecovery: -")
+ print rnazString
+ print cmStatString
+ else do
+ putStrLn (show benchmarkIndex ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ show alienMaxLinkscore ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ show alienFastaEntriesNumber ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ rnazString ++ "\t" ++ cmStatString)
+
+rnazOutput :: Verbosity -> String -> IO String
+rnazOutput verbose rnazPath = do
+ rnazPresent <- doesFileExist rnazPath
+ if rnazPresent
+ then do
+ inputRNAz <- readRNAz rnazPath
+ if isRight inputRNAz
+ then do
+ let rnaZ = fromRight inputRNAz
+ if (verbose == Loud)
+ then do
+ let output = "Mean pairwise identity: " ++ show (meanPairwiseIdentity rnaZ) ++ "\n Shannon entropy: " ++ show (shannonEntropy rnaZ) ++ "\n GC content: " ++ show (gcContent rnaZ) ++ "\n Mean single sequence minimum free energy: " ++ show (meanSingleSequenceMinimumFreeEnergy rnaZ) ++ "\n Consensus minimum free energy: " ++ show (consensusMinimumFreeEnergy rnaZ) ++ "\n Energy contribution: " ++ show (energyContribution rnaZ) ++ "\n Covariance contribution: " ++ show (covarianceContribution rnaZ) ++ "\n Combinations pair: " ++ show (combinationsPair rnaZ) ++ "\n Mean z-score: " ++ show (meanZScore rnaZ) ++ "\n Structure conservation index: " ++ show (structureConservationIndex rnaZ) ++ "\n Background model: " ++ backgroundModel rnaZ ++ "\n Decision model: " ++ decisionModel rnaZ ++ "\n SVM decision value: " ++ show (svmDecisionValue rnaZ) ++ "\n SVM class propability: " ++ show (svmRNAClassProbability rnaZ) ++ "\n Prediction: " ++ (prediction rnaZ)
+ return output
+ else do
+ let output = show (meanPairwiseIdentity rnaZ) ++ "\t" ++ show (shannonEntropy rnaZ) ++ "\t" ++ show (gcContent rnaZ) ++ "\t" ++ show (meanSingleSequenceMinimumFreeEnergy rnaZ) ++ "\t" ++ show (consensusMinimumFreeEnergy rnaZ) ++ "\t" ++ show (energyContribution rnaZ) ++ "\t" ++ show (covarianceContribution rnaZ) ++ "\t" ++ show (combinationsPair rnaZ) ++ "\t" ++ show (meanZScore rnaZ) ++ "\t" ++ show (structureConservationIndex rnaZ) ++ "\t" ++ show (svmDecisionValue rnaZ) ++ "\t" ++ show (svmRNAClassProbability rnaZ) ++ "\t" ++ (prediction rnaZ)
+ return output
+ else do
+ if (verbose == Loud)
+ then do
+ let output = "Mean pairwise identity: " ++ " - \n Shannon entropy: " ++ " - \n GC content: " ++ " - \n Mean single sequence minimum free energy: " ++ " - \n Consensus minimum free energy: " ++ " - \n Energy contribution: " ++ " - \n Covariance contribution: " ++ " - \n Combinations pair: " ++ " - \n Mean z-score: " ++ " - \n Structure conservation index: " ++ " - \n Background model: " ++ " - \n Decision model: " ++ " - \n SVM decision value: " ++ " - \n SVM class propability: " ++ " - \n Prediction: " ++ " - \n"
+ return output
+ else do
+ let output = "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-"
+ return output
+ else do
+ if (verbose == Loud)
+ then do
+ let output = "Mean pairwise identity: " ++ " - \n Shannon entropy: " ++ " - \n GC content: " ++ " - \n Mean single sequence minimum free energy: " ++ " - \n Consensus minimum free energy: " ++ " - \n Energy contribution: " ++ " - \n Covariance contribution: " ++ " - \n Combinations pair: " ++ " - \n Mean z-score: " ++ " - \n Structure conservation index: " ++ " - \n Background model: " ++ " - \n Decision model: " ++ " - \n SVM decision value: " ++ " - \n SVM class propability: " ++ " - \n Prediction: " ++ " - \n"
+ return output
+ else do
+ let output = "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-" ++ "\t" ++ "-"
+ return output
+
+cmStatOutput :: Verbosity -> String -> IO String
+cmStatOutput verbose cmstatPath = do
+ cmstatPresent <- doesFileExist cmstatPath
+ if cmstatPresent
+ then do
+ inputCMstat <- readCMstat cmstatPath
+ if isRight inputCMstat
+ then do
+ let cmStat = fromRight inputCMstat
+ if (verbose == Loud)
+ then do
+ let output = "statSequenceNumber: " ++ (show (statSequenceNumber cmStat)) ++ "\nstatEffectiveSequences: " ++ (show (statEffectiveSequences cmStat)) ++ "\nstatConsensusLength: " ++ (show (statConsensusLength cmStat)) ++ "\nstatW: " ++ show (statW cmStat) ++ "\nstatBasepairs: " ++ show (statBasepairs cmStat) ++ "\nstatBifurcations: " ++ (show (statBifurcations cmStat)) ++ "\nstatModel: " ++ (statModel cmStat) ++ "\nrelativeEntropyCM: " ++ show (relativeEntropyCM cmStat) ++ "\nrelativeEntropyHMM: " ++ show (relativeEntropyHMM cmStat)
+ return output
+ else do
+ let output = (show (statSequenceNumber cmStat)) ++ "\t" ++ (show (statEffectiveSequences cmStat)) ++ "\t" ++ (show (statConsensusLength cmStat)) ++ "\t" ++ show (statW cmStat) ++ "\t" ++ show (statBasepairs cmStat) ++ "\t" ++ (show (statBifurcations cmStat)) ++ "\t" ++ (statModel cmStat) ++ "\t" ++ show (relativeEntropyCM cmStat) ++ "\t" ++ show (relativeEntropyHMM cmStat)
+ return output
+ else do
+ if (verbose == Loud)
+ then do
+ let output = "statSequenceNumber: -" ++ "\nstatEffectiveSequences: -" ++ "\nstatConsensusLength: -" ++ "\nstatW: -" ++ "\nstatBasepairs: -" ++ "\nstatBifurcations: -" ++ "\nstatModel: -" ++ "\nrelativeEntropyCM: -" ++ "\nrelativeEntropyHMM: -"
+ return output
+ else do
+ let output = "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-"
+ return output
+ else do
+ if (verbose == Loud)
+ then do
+ let output = "statSequenceNumber: -" ++ "\nstatEffectiveSequences: -" ++ "\nstatConsensusLength: -" ++ "\nstatW: -" ++ "\nstatBasepairs: -" ++ "\nstatBifurcations: -" ++ "\nstatModel: -" ++ "\nrelativeEntropyCM: -" ++ "\nrelativeEntropyHMM: -"
+ return output
+ else do
+ let output = "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-\t" ++ "-"
+ return output
+
+
diff --git a/src/Bio/SelectSequences.hs b/src/Bio/SelectSequences.hs
new file mode 100644
--- /dev/null
+++ b/src/Bio/SelectSequences.hs
@@ -0,0 +1,37 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+
+-- | Select Sequences
+-- Testcommand: SelectSequences -i /path/to/test.clustal
+module Main where
+
+import System.Console.CmdArgs
+import Bio.RNAlienLibrary
+import Data.Either.Unwrap
+
+data Options = Options
+ { inputClustalPath :: String,
+ toogleExternalSelectSequences :: Bool
+ } deriving (Show,Data,Typeable)
+
+options :: Options
+options = Options
+ { inputClustalPath = def &= name "i" &= help "Path to input clustal file",
+ toogleExternalSelectSequences = False &= name "e" &= help "Use only replacement of alignment characters and external 'selectSequence.pl'. Default: False"
+ } &= summary "SelectSequences" &= help "Florian Eggenhofer 2015" &= verbosity
+
+main :: IO ()
+main = do
+ Options{..} <- cmdArgs options
+ let reformatedClustalPath = inputClustalPath ++ ".reformated"
+ if toogleExternalSelectSequences
+ then do
+ resultStatus <- preprocessClustalForRNAzExternal inputClustalPath reformatedClustalPath
+ if (isRight resultStatus)
+ then (return ())
+ else (print ("A problem occured selecting sequences: " ++ fromLeft resultStatus))
+ else do
+ resultStatus <- preprocessClustalForRNAz inputClustalPath reformatedClustalPath
+ if (isRight resultStatus)
+ then (return ())
+ else (print ("A problem occured selecting sequences: " ++ fromLeft resultStatus))