RNAlien (empty) → 1.0.0
raw patch · 8 files changed
+3108/−0 lines, 8 filesdep +BlastHTTPdep +ClustalParserdep +EntrezHTTPsetup-changed
Dependencies added: BlastHTTP, ClustalParser, EntrezHTTP, HTTP, RNAlien, Taxonomy, ViennaRNAParser, base, biocore, biofasta, blastxml, bytestring, cassava, cmdargs, containers, directory, edit-distance, either-unwrap, filepath, hierarchical-clustering, matrix, parsec, process, random, split, time, vector
Files
- LICENSE +675/−0
- RNAlien.cabal +75/−0
- Setup.hs +2/−0
- src/Bio/RNAlien.hs +100/−0
- src/Bio/RNAlienData.hs +144/−0
- src/Bio/RNAlienLibrary.hs +1819/−0
- src/Bio/RNAlienStatistics.hs +256/−0
- src/Bio/SelectSequences.hs +37/−0
+ LICENSE view
@@ -0,0 +1,675 @@+ GNU GENERAL PUBLIC LICENSE+ Version 3, 29 June 2007++ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>+ 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.++ <one line to give the program's name and a brief idea of what it does.>+ Copyright (C) <year> <name of author>++ 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 <http://www.gnu.org/licenses/>.++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:++ <program> Copyright (C) <year> <name of author>+ 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+<http://www.gnu.org/licenses/>.++ 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+<http://www.gnu.org/philosophy/why-not-lgpl.html>.+
+ RNAlien.cabal view
@@ -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):+ .+ * <https://github.com/eggzilla/RNAlien RNAlien>+ .+ * <https://github.com/eggzilla/AlienServer AlienServer>+ .+ TaxonomyTools which can be used to visualise the organisms included in a RNAlien result can be found here (License GPL-3):+ .+ * <https://github.com/eggzilla/TaxonomyTools TaxonomyTools-Github>+ .+ * <https://hackage.haskell.org/package/Taxonomy TaxonomyTools-Hackage>+ .+ For instruction how to use RNAlien please see the <http://rna.tbi.univie.ac.at/rnalien/help Help page>.+ .+ Dependencies:+ .+ * <http://infernal.janelia.org/ Infernal>+ .+ * <http://www.clustal.org/omega/#Download clustal-omega>+ .+ * <http://www.bioinf.uni-freiburg.de/Software/LocARNA/#download Locarna>+ .+ * <https://www.tbi.univie.ac.at/~wash/RNAz/ RNAz>+ .+ * <http://www.tbi.univie.ac.at/RNA/index.html#download ViennaRNA package>+ .+ 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++
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ src/Bio/RNAlien.hs view
@@ -0,0 +1,100 @@+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveDataTypeable #-}++-- | Unsupervized construction of RNA family models+-- For more information on RNA family models consult <http://>+-- 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") ""
+ src/Bio/RNAlienData.hs view
@@ -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"+
+ src/Bio/RNAlienLibrary.hs view
@@ -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
+ src/Bio/RNAlienStatistics.hs view
@@ -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++
+ src/Bio/SelectSequences.hs view
@@ -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))