From 414f67e32564a8015d94060616abea51015b32ad Mon Sep 17 00:00:00 2001 From: Harvey Huang Date: Thu, 28 Sep 2023 16:23:36 -0500 Subject: [PATCH] Initial commit of all deidentified code --- .gitignore | 4 + LICENSE.txt | 674 +++++++++++++++++++++++++++++ README.txt | 29 ++ applyCARLARealCCEPs.m | 217 ++++++++++ applyCARLARealCCEPsLoop.m | 153 +++++++ fig1_simCCEPComponents.m | 182 ++++++++ fig2_CARLAonSingleTrial.m | 144 ++++++ fig3_EPConstruct.m | 128 ++++++ fig7_summarizeCARLARealCCEPs.m | 162 +++++++ fig8_summarizeInterChCorr.m | 125 ++++++ figS1_compareCARvCARLA.m | 55 +++ figS3_simCCEPComponentsGlobalSig.m | 191 ++++++++ functions/CARLA.m | 225 ++++++++++ functions/antifind.m | 29 ++ functions/crossChCorr.m | 105 +++++ functions/genRandSig.m | 52 +++ functions/genRandSigOrig.m | 51 +++ functions/getNeighborChs.m | 57 +++ functions/rmBadTrialsAnnots.m | 57 +++ main.m | 179 ++++++++ simCCEPLoop.m | 175 ++++++++ summarizeSimCCEPLoop.m | 197 +++++++++ tab1_dispNumElectrodes.m | 55 +++ 23 files changed, 3246 insertions(+) create mode 100644 .gitignore create mode 100644 LICENSE.txt create mode 100644 README.txt create mode 100644 applyCARLARealCCEPs.m create mode 100644 applyCARLARealCCEPsLoop.m create mode 100644 fig1_simCCEPComponents.m create mode 100644 fig2_CARLAonSingleTrial.m create mode 100644 fig3_EPConstruct.m create mode 100644 fig7_summarizeCARLARealCCEPs.m create mode 100644 fig8_summarizeInterChCorr.m create mode 100644 figS1_compareCARvCARLA.m create mode 100644 figS3_simCCEPComponentsGlobalSig.m create mode 100644 functions/CARLA.m create mode 100644 functions/antifind.m create mode 100644 functions/crossChCorr.m create mode 100644 functions/genRandSig.m create mode 100644 functions/genRandSigOrig.m create mode 100644 functions/getNeighborChs.m create mode 100644 functions/rmBadTrialsAnnots.m create mode 100644 main.m create mode 100644 simCCEPLoop.m create mode 100644 summarizeSimCCEPLoop.m create mode 100644 tab1_dispNumElectrodes.m diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..2d407a9 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +/output +/data +/ignore +*.DS_Store \ No newline at end of file diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..e72bfdd --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. \ No newline at end of file diff --git a/README.txt b/README.txt new file mode 100644 index 0000000..a2f4765 --- /dev/null +++ b/README.txt @@ -0,0 +1,29 @@ +If this code is used in a publication, please cite the manuscript: +"CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, GA Worrell, KJ Miller, and D Hermes. +A preprint is available currently at doi: ___ + +Correspondence: +H Huang: huang.harvey@mayo.edu; D Hermes: hermes.dora@mayo.edu + +***** + +DEPENDENCIES + +- mnl_ieegBasics: https://github.com/MultimodalNeuroimagingLab/mnl_ieegBasics +- matmef: https://github.com/MultimodalNeuroimagingLab/matmef +- mnl_seegview: https://github.com/MultimodalNeuroimagingLab/mnl_seegview +- SPM12 (as dependency for mnl_seegview): https://www.fil.ion.ucl.ac.uk/spm/software/spm12/ +- vistasoft: https://github.com/vistalab/vistasoft +- Freesurfer v7: https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall +- GIfTI: https://github.com/gllmflndn/gifti + +***** + +USAGE + +a. The data used in this analysis will be available upon manuscript publication on OpenNeuro, in BIDS format. Please download all the data and copy into a folder named "data", in the root level code directory. In the meantime before data release, all results pertaining to simulated CCEPs can still be reproduced by the code alone. + +b. Pial, cortical, and subcortical segmentations for each subject were obtained using Freesurfer v7. These were used to localize each stimulation site to a particular tissue type. The relevant Freesurfer outputs for each subject are located in the data subdirectory: data/derivatives/freesurfer + +c. All custom analyses were performed in MATLAB R2023a. Step-by-step code blocks and instructions to generate all manuscript figures and results are in "main.m". Please open "main.m" and follow the code sections contained therein. \ No newline at end of file diff --git a/applyCARLARealCCEPs.m b/applyCARLARealCCEPs.m new file mode 100644 index 0000000..e179bfb --- /dev/null +++ b/applyCARLARealCCEPs.m @@ -0,0 +1,217 @@ +%% This script generates outputs for CARLA applied to single stim sites in real data (Fig 5B-F, Fig 6C-F) +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +%% Load all data for selected subject + +% 'sub' and 'run' variables should have been defined in main script + +elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub)); +elecs = ieeg_readtableRmHyphens(elecsPath); + +mefPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run)); +channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run)); +eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run)); +annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run)); + +% Load CCEP data +mef = ccep_PreprocessMef(mefPath, channelsPath, eventsPath); +mef.loadMefAll; +mef.highpass; +mef.loadMefTrials([-1, 1]); + +mef.plotOutputs(1, [], 200); + +tt = mef.tt; +srate = mef.srate; +channels = mef.channels; +events = mef.evts; +data = mef.data; + +sites = groupby(events.electrical_stimulation_site); + +outdir = fullfile('output', 'realCCEPs'); +mkdir(outdir); + +try + annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t'); + annots = annots.status_description; + assert(length(annots) == height(events), 'Error: length of event annotations does not match events table height'); +catch + warning('Unable to load event annotation for events'); +end + +% show the seizure onset zones which are excluded from stimulation +SOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ')); +disp('SOZs:'); +disp(SOZ'); + +clear *Path* + +%% Extract data from site of interest, NaN out individual bad trials + +% 'site' variable should have been defined in main script + +idxesSite = sites{strcmp(site, sites(:, 1)), 2}; % indices corresponding to this stim site +nTrials = length(idxesSite); + +goodChs = strcmp(channels.type, 'SEEG') & strcmp(channels.status, 'good') & ~ismember(upper(channels.name), getNeighborChs(split(site,'-'), 2)); % exclude sites 2 away from current + +dataSite = data(goodChs, :, idxesSite); +chsSite = channels.name(goodChs); + +% remove bad trials and set channels to nan in individual trials, according to annotations +annotsSite = annots(idxesSite); % get annotations corresponding to this stim site +dataSite = rmBadTrialsAnnots(dataSite, chsSite, annotsSite); + +figure('Position', [200, 200, 1200, 1200]); +subplot(1, 2, 1); % first half of channels +ieeg_plotTrials(tt, mean(dataSite(1:floor(length(chsSite)/2), :, :), 3, 'omitnan')', 200, chsSite(1:floor(length(chsSite)/2))); +xlim([-0.5, 1]); +subplot(1, 2, 2); % first half of channels +ieeg_plotTrials(tt, mean(dataSite(ceil(length(chsSite)/2):end, :, :), 3, 'omitnan')', 200, chsSite(ceil(length(chsSite)/2):end)); +xlim([-0.5, 1]); + +%% Apply CARLA and plot variance, zminmean plots + +rng('default'); + +nanChs = any(isnan(dataSite), [2, 3]); % any channels with nan in any trial + +dataSiteNoNan = dataSite(~nanChs, :, :); +chsSiteNoNan = chsSite(~nanChs); + +[Vcar, CAR, stats] = CARLA(tt, dataSiteNoNan, srate, true); +nCAR = length(stats.chsUsed); + +cmSens = [1, 165/255, 0]; % orange color for max + +% ZminMean plot +figure('Position', [200, 200, 600, 300]); hold on +errorbar(mean(stats.zMinMean, 2), std(stats.zMinMean, 0, 2), 'k.-', 'MarkerSize', 10, 'CapSize', 1); % SD of bootstrapped means +plot(nCAR, mean(stats.zMinMean(nCAR, :), 2), '*', 'Color', cmSens); +xlim([0, length(chsSiteNoNan)+1]); ylim([-1.5, 0]); +yline(0, 'Color', 'k'); +saveas(gcf, fullfile(outdir, sprintf('%s_%s_zMinMean', sub, site)), 'svg'); +saveas(gcf, fullfile(outdir, sprintf('%s_%s_zMinMean', sub, site)), 'png'); + +% Variance plot +vars = stats.vars(stats.order); % sort in order +figure('Position', [200, 200, 600, 300]); hold on +plot(vars, 'k.-', 'MarkerSize', 10); +xline(nCAR + 0.5, 'Color', cmSens); +xlim([0, length(chsSiteNoNan)+1]); +saveas(gcf, fullfile(outdir, sprintf('%s_%s_covar', sub, site)), 'png'); +saveas(gcf, fullfile(outdir, sprintf('%s_%s_covar', sub, site)), 'svg'); + +%% Plot channels sorted by CARLA order and CAR + +% Plot (pre-CAR) data with line noise removed, in order of variance +dataSite2Plot = zeros(size(dataSiteNoNan)); +for tr = 1:size(dataSiteNoNan, 3) + dataSite2Plot(:, :, tr) = ieeg_notch(dataSiteNoNan(:, :, tr)', srate, 60)'; +end +dataSite2Plot = dataSite2Plot(stats.order, :, :); + +yspace = 120; + +% make excluded channels gray +cm = zeros(length(chsSiteNoNan), 3); +cm(nCAR+1:end, :) = repmat([0.5, 0.5, 0.5], length(chsSiteNoNan)-nCAR, 1); + +figure('Position', [200, 200, 1200, 1200]); +nHalf = floor(length(chsSiteNoNan)/2); % number at which to break channels into 2 columns +subplot(1, 2, 1); % first half of channels +ieeg_plotTrials(tt, mean(dataSite2Plot(1:nHalf, :, :), 3)', yspace, 1:nHalf, cm(1:nHalf, :)); +xlim([-0.1, 0.5]); +subplot(1, 2, 2); % second half of channels +ieeg_plotTrials(tt, mean(dataSite2Plot(nHalf+1:end, :, :), 3)', yspace, nHalf+1:length(chsSiteNoNan), cm(nHalf+1:end, :)); +xlim([-0.1, 0.5]); +saveas(gcf, fullfile(outdir, sprintf('%s_%s_sortedChs', sub, site)), 'png'); +print(gcf, fullfile(outdir, sprintf('%s_%s_sortedChs', sub, site)),'-depsc2', '-r300', '-painters'); + +% Plot the adjusted common average itself +figure('Position', [200, 200, 600, 200]); +hold on +plot(tt, CAR, 'Color', [0.5, 0.5, 0.5]); +plot(tt, mean(CAR, 2), 'k-', 'LineWidth', 1.5); +xlim([-0.1, 0.5]); ylim([-1200, 1200]); +saveas(gcf, fullfile(outdir, sprintf('%s_%s_CAR', sub, site)), 'png'); +print(gcf, fullfile(outdir, sprintf('%s_%s_CAR', sub, site)), '-depsc2', '-r300', '-painters'); + +%% Plot all channels in one column (for Figure S1) + +% downsample plotting signals to make less taxing on illustrator +dataSite2PlotDs = nan(size(dataSite2Plot, 1), length(tt)/4, size(dataSite2Plot, 3)); +for ii = 1:size(dataSite2Plot, 3) + dataSite2PlotDs(:, :, ii) = downsample(dataSite2Plot(:, :, ii)', 4)'; +end +ttDs = tt(1:4:end); % downsample time vector + +yspace = 80; + +figure('Position', [200, 200, 400, 1200]); +ys = ieeg_plotTrials(ttDs, mean(dataSite2PlotDs, 3)', yspace, 1:length(chsSiteNoNan), cm); +xlim([-0.1, 0.5]); ylim([ys(end) - 7*yspace, ys(1) + 3*yspace]); +saveas(gcf, fullfile(outdir, sprintf('%s_%s_sortedChs_column', sub, site)), 'png'); +print(gcf, fullfile(outdir, sprintf('%s_%s_sortedChs_column', sub, site)), '-depsc2', '-r300', '-painters'); + +%% Plot examples of channels before and after CARLA + +% chs2Plot should be defined in main script + +for ii = 1:length(chs2Plot) + + ch = chs2Plot(ii); + + fprintf('Ch %s\n', chsSiteNoNan{stats.order(ch)}); + + sigRaw = squeeze(dataSiteNoNan(stats.order(ch), :, :)); % only high-pass, no notch + sigNotch = squeeze(dataSite2Plot(ch, :, :)); % only notch, no re-reference (order already sorted) + sigCar = squeeze(Vcar(stats.order(ch), :, :)); + + % plot raw signal + figure('Position', [200, 200, 450, 450]); + subplot(3, 1, 1); hold on; + plot(tt, sigRaw, 'Color', [0.5, 0.5, 0.5]); + plot(tt, mean(sigRaw, 2), 'k-', 'LineWidth', 1.5); + xlim([-0.1, 0.5]); ylim([-1200, 1200]); + + % plot notched signal + subplot(3, 1, 2); hold on; + plot(tt, sigNotch, 'Color', [0.5, 0.5, 0.5]); + plot(tt, mean(sigNotch, 2), 'k-', 'LineWidth', 1.5); + xlim([-0.1, 0.5]); ylim([-600, 600]); + + % plot CARLA signal + subplot(3, 1, 3); hold on; + plot(tt, sigCar, 'Color', [0.5, 0.5, 0.5]); + plot(tt, mean(sigCar, 2), 'k-', 'LineWidth', 1.5); + xlim([-0.1, 0.5]); ylim([-600, 600]); + + saveas(gcf, fullfile(outdir, sprintf('%s_%s_exCh%d', sub, site, ch)), 'png'); + print(gcf, fullfile(outdir, sprintf('%s_%s_exCh%d', sub, site, ch)), '-depsc2', '-r300', '-painters'); + +end + diff --git a/applyCARLARealCCEPsLoop.m b/applyCARLARealCCEPsLoop.m new file mode 100644 index 0000000..0bc346a --- /dev/null +++ b/applyCARLARealCCEPsLoop.m @@ -0,0 +1,153 @@ +%% This loops through stim sites for each subject and saves CARLA results, with fewer other outputs along the way +% Then run fig7_summarizeCARLARealCCEPs to summarize CARLA results for Figure 7 +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +%% Paths, load data for each CCEP run + +outdir = fullfile('output', 'realCCEPsLoop'); + +% subjects and runs corresponding to each subject to load +subs = {'sub-1', 'sub-2', 'sub-3', 'sub-4'}; +runs = {[1, 2, 3, 5], [12, 13, 14], [2, 3, 4], [1, 2]}; + +cmSens = [1, 165/255, 0]; % orange color used for detecting optimum + +for ss = 1:length(subs) + + sub = subs{ss}; + mkdir(outdir, sub); + + % iterate across runs in subject + for rr = 1:length(runs{ss}) + + run = runs{ss}(rr); + + mefPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run)); + channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run)); + eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run)); + annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run)); + elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub)); + + % Load CCEP data from each subject + mef = ccep_PreprocessMef(mefPath, channelsPath, eventsPath); + mef.loadMefAll; + mef.highpass; + mef.loadMefTrials([-1, 1]); + + tt = mef.tt; + srate = mef.srate; + channels = mef.channels; + events = mef.evts; + data = mef.data; + + sites = groupby(events.electrical_stimulation_site); + + % read manual annotations on events + try + annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t'); + annots = annots.status_description; + catch + warning('Unable to load event annotation for events'); + end + + % check that annots match events and remove bad events from annots + if exist('annots', 'var') + eventsTemp = ieeg_readtableRmHyphens(eventsPath, 'electrical_stimulation_site', 1); % to remove bad events from annotations + assert(length(annots) == height(eventsTemp), 'Error: length of event annotations does not match events table height'); + annots(~strcmpi(eventsTemp.status, 'good')) = []; + end + + elecs = ieeg_readtableRmHyphens(elecsPath); + SOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ')); + disp('SOZs:'); + disp(SOZ'); + + %% Iterate across sites in current run + + for kk = 1:size(sites, 1) + + site = sites{kk}; + idxesSite = sites{kk, 2}; + + if any(ismember(split(site, '-'), SOZ)) + fprintf('Skipping stim site %s in SOZ\n', site); + continue + end + fprintf('Current site: %s\n', site); + + % extract good channels at this stim site + goodChs = strcmp(channels.type, 'SEEG') & strcmp(channels.status, 'good') & ~ismember(upper(channels.name), getNeighborChs(split(site,'-'), 2)); % exclude sites 2 away from current + dataSite = data(goodChs, :, idxesSite); + chsSite = channels.name(goodChs); + + try + % remove bad trials and set channels to nan in singular trials, according to annotations + annotsSite = annots(idxesSite); % get annotations corresponding to this stim site + dataSite = rmBadTrialsAnnots(dataSite, chsSite, annotsSite); + end + + % minimum floor on number of trials for any stim site + if size(dataSite, 3) < 8 + fprintf('Skipping stim site %s with %d trials\n', site, size(dataSite, 3)); + continue + end + + %% Remove channels that have any nan trials, apply CARLA + + rng('default'); + + nanChs = any(isnan(dataSite), [2, 3]); % any channels with nan in any trial + dataSiteNoNan = dataSite(~nanChs, :, :); + chsSiteNoNan = chsSite(~nanChs); + + % calculate average interchannel R-squared at each possible adjusted CAR size + rsqs = crossChCorr(tt, dataSiteNoNan, srate); + save(fullfile(outdir, sub, sprintf('%s_%s_interChRsq.mat', sub, site)), 'rsqs'); + + % Perform CAR + [Vcar, CAR, stats] = CARLA(tt, dataSiteNoNan, srate, true); + nCAR = length(stats.chsUsed); + + % ZminMean plot + f = figure('Position', [200, 200, 600, 300]); hold on + errorbar(mean(stats.zMinMean, 2), std(stats.zMinMean, 0, 2), 'k.-', 'MarkerSize', 10, 'CapSize', 1); % SD of bootstrapped means + plot(nCAR, mean(stats.zMinMean(nCAR, :), 2), '*', 'Color', cmSens); + xlim([0, length(chsSiteNoNan)+1]); ylim([-inf, 0]); + yline(0, 'Color', 'k'); + saveas(gcf, fullfile(outdir, sub, sprintf('%s_%s_zMinMean', sub, site)), 'svg'); + saveas(gcf, fullfile(outdir, sub, sprintf('%s_%s_zMinMean', sub, site)), 'png'); + close(f); + + % save the optimum CAR size and number of total channels to text + fid = fopen(fullfile(outdir, sub, sprintf('%s_%s_nCAR.txt', sub, site)), 'w'); + fprintf(fid, '%d\n%d', nCAR, length(chsSiteNoNan)); + fclose(fid); + + + end + + end + +end \ No newline at end of file diff --git a/fig1_simCCEPComponents.m b/fig1_simCCEPComponents.m new file mode 100644 index 0000000..b72b586 --- /dev/null +++ b/fig1_simCCEPComponents.m @@ -0,0 +1,182 @@ +%% Figure 1. Create simulated CCEP data and their components +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% + +srate = 4800; +tt = -0.5:1/srate:1-1/srate; + +nchs = 10; % number of simulated channels +ntrs = 12; % number of trials +chs = arrayfun(@(x) sprintf('ch%d', x), 1:nchs, 'UniformOutput', false)'; + +outdir = fullfile('output', 'simComponents'); +mkdir(outdir) + +%% i) Create data, add stim artifacts + +% stores the data +V0 = zeros(length(tt), nchs); + +Aart = 50 + rand(nchs, 1)*6; % slightly different artifact amplitudes for each channel +artifact = sin(2*pi*600*tt)'; +artifact(tt < 0 | tt > 0.002) = 0; +V0 = V0 + artifact*Aart'; + +figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V0, 100); +xlabel('Time (s)'); ylabel('Channels'); + +%% A) Add unique input signals at a subset of channels + +rng(14); + +chsResp = 1:4; + +A = 100; +% using previous version of signal generator for this illustration to preserve backwards compatibility +sig = genRandSigOrig(tt, length(chsResp), A); + +V1 = V0; +V1(:, chsResp) = V0(:, chsResp) + sig; + +figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V1, 100, [], [], 'LineWidth', 1.5); +xlabel('Time (s)'); ylabel('Channels'); + +%% B) Add a global signal across all chs and trials when present. Here amplitude Aglobal = 0, so no global signal added + +rng(10); + +Aglobal = 0; + +sigGlob = genRandSigOrig(tt, 1, Aglobal); + +V2 = V1 + sigGlob; + +figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V2, 100, [], [], 'LineWidth', 1.5); +xlabel('Time (s)'); ylabel('Channels'); + +%% C) Add noise common to all channels (line noise and some brown noise from reference) + +rng('default'); + +% add trials, so now put V in ch x time x trial format +V3 = repmat(V2', 1, 1, ntrs); + +% random line noise phases for each trial, trials x harmonic (60, 120, 180) +phLN = rand(ntrs, 3)*2*pi; +LN = zeros(length(tt), ntrs); +for ii = 1:ntrs + LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3)); +end + +% brown noise shared across channels (from ref electrode +BN = cumsum(0.4*randn(2*length(tt), ntrs)); +BN = ieeg_highpass(BN, srate, true); +BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :); + +noiseCommon = LN + BN; + +V3 = V3 + shiftdim(noiseCommon, -1); + +figure('Position', [200, 200, 600, 400]); +subplot(1, 2, 1); ieeg_plotTrials(tt, V3(:, :, 1)', 100, [], [], 'LineWidth', 1.5); +title('V at one trial'); +xlabel('Time (s)'); ylabel('Channels'); + +subplot(1, 2, 2); ieeg_plotTrials(tt, mean(V3, 3)', 100, [], [], 'LineWidth', 1.5); +title('V across trials'); +xlabel('Time (s)'); ylabel('Channels'); + + +%% D) add random brown noise across all channels + +rng('default'); + +noiseRand = cumsum(0.4*randn(nchs, 2*length(tt), ntrs), 2); % give double the number of time points so we can highpass it +for ii = 1:nchs + noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true); +end +noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :); + +V4 = V3 + noiseRand; + +figure('Position', [200, 200, 600, 400]); +subplot(1, 2, 1); ieeg_plotTrials(tt, V4(:, :, 1)', 100, [], [], 'LineWidth', 1.5); +title('V at one trial'); +xlabel('Time (s)'); ylabel('Channels'); + +subplot(1, 2, 2); ieeg_plotTrials(tt, mean(V4, 3)', 100, [], [], 'LineWidth', 1.5); +title('V across trials'); +xlabel('Time (s)'); ylabel('Channels'); + +saveas(gcf, fullfile(outdir, 'V4'), 'svg'); +saveas(gcf, fullfile(outdir, 'V4'), 'png'); + + +%% Save construction components for channels of interest + +% Plot select channels for figure (3 and 8) +chs2Plot = [3, 8]; + +for ii = 1:length(chs2Plot) + + ch = chs2Plot(ii); + + ylims = [-120, 120]; + + % plot full simulation of all trials at data + figure('Position', [200, 200, 600, 300]); hold on + plot(tt, squeeze(V4(ch, :, :)), 'Color', [0.5, 0.5, 0.5]); + plot(tt, mean(squeeze(V4(ch, :, :)), 2), 'k', 'LineWidth', 1.5); + hold off + ylim(ylims); + xlabel('Time (s)'); ylabel('Voltage (\muV)'); + saveas(gcf, fullfile(outdir, sprintf('V4_ch%d', ch)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('V4_ch%d', ch)), 'svg'); + + % plot the true signal + figure('Position', [200, 200, 200, 300]); + plot(tt, V1(:, ch), 'k', 'LineWidth', 2); + ylim(ylims); xlim([0, 0.5]); + xlabel('Time (s)'); ylabel('Voltage (\muV)'); + saveas(gcf, fullfile(outdir, sprintf('V1_ch%d', ch)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('V1_ch%d', ch)), 'svg'); + + % plot examples of the common noise + figure('Position', [200, 200, 200, 300]); + ieeg_plotTrials(tt, noiseCommon(:, 1:3), 50, [], [0.5, 0.5, 0.5]); % plot at one trial + xlim([0, 0.5]); + xlabel('Time (s)'); ylabel('Voltage (\muV)'); + saveas(gcf, fullfile(outdir, 'lineNoise_3trs'), 'png'); + saveas(gcf, fullfile(outdir, 'lineNoise_3trs'), 'svg'); + + % plot the random noise + figure('Position', [200, 200, 200, 300]); hold on + plot(tt, squeeze(noiseRand(ch, :, :)), 'Color', [0.5, 0.5, 0.5]); + hold off + xlim([0, 0.5]); ylim(ylims); + xlabel('Time (s)'); ylabel('Voltage (\muV)'); + saveas(gcf, fullfile(outdir, sprintf('randNoise_ch%d', ch)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('randNoise_ch%d', ch)), 'svg'); +end \ No newline at end of file diff --git a/fig2_CARLAonSingleTrial.m b/fig2_CARLAonSingleTrial.m new file mode 100644 index 0000000..5d01df3 --- /dev/null +++ b/fig2_CARLAonSingleTrial.m @@ -0,0 +1,144 @@ +%% Figure 2: Methods illustration of how CARLA works, on single-trial simulated CCEP data +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +%% A) Create data, add stim artifacts and responses to the last 2 channels + +srate = 4800; +tt = -0.5:1/srate:1-1/srate; + +rng(50); + +nchs = 6; % number of simulated channels +ntrs = 1; % number of trials +chs = arrayfun(@(x) sprintf('ch%d', x), 1:nchs, 'UniformOutput', false)'; + +% stores the data +V0 = zeros(length(tt), nchs); + +Aart = 50 + rand(nchs, 1)*5; % slightly different artifact amplitudes for each channel +artifact = sin(2*pi*600*tt)'; +artifact(tt < 0 | tt > 0.002) = 0; +V0 = V0 + artifact*Aart'; + +chsResp = [nchs-1, nchs]; + +sig = genRandSig(tt, length(chsResp), 150); + +V1 = V0; +V1(:, chsResp) = V0(:, chsResp) + sig; + +figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V1, 80, [], [], 'LineWidth', 1.5); +xlabel('Time (s)'); ylabel('Channels'); + + +%% B) Add common noise (line noise and reference brown noise) to each trial + +rng('default'); + +% no line noise in this demonstration, pretend that this has been removed by filtering +% Create line noise: random noise phase for each harmonic (60, 120, 180) +%phNoise = rand(1, 3)*2*pi; +%LN = sin(2*pi*60*tt - phNoise(1))' + 0.3*sin(2*pi*120*tt - phNoise(2))' + 0.1*sin(2*pi*180*tt - phNoise(3))'; + +% Create some brown noise and low-pass as common signal +brownNoiseRef = cumsum(0.4*randn(2*length(tt), 1)); +brownNoiseRef = ieeg_highpass(brownNoiseRef, srate, true); +brownNoiseRef = brownNoiseRef((0.5*length(tt)+1) : 1.5*length(tt)); + +V3 = V1 + brownNoiseRef; + +figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V3, 80, [], [], 'LineWidth', 1.5); +xlabel('Time (s)'); ylabel('Channels'); + + +%% C) add random brown noise to each channel + +rng(4); + +noiseRand = cumsum(0.4*randn(2*length(tt), nchs)); % give double the number of time points so we can highpass it +noiseRand = ieeg_highpass(noiseRand, srate, true); + +noiseRand = noiseRand((0.5*length(tt)+1) : 1.5*length(tt), :); + +V4 = V3 + noiseRand; + +figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V4, 80, [], [], 'LineWidth', 1.5); +xlim([-0.1, 0.5]) +xlabel('Time (s)'); ylabel('Channels'); + +%% Apply CARLA and show outputs + +rng('default'); + +[Vcar, CAR, stats] = CARLA(tt, V4, srate); + +outdir = fullfile('output', 'simulationTrial'); +mkdir(outdir); + +% Plot in increasing order of variance (Figure 2A) +V4sorted = V4(:, stats.order); +figure('Position', [200, 200, 200, 500]); ieeg_plotTrials(tt, V4sorted, 100, [], [], 'LineWidth', 1.5); +xlabel('Time (s)'); ylabel('Channels'); +xlim([-0.1, 0.5]); +xline(0.01, 'Color', 'r'); +xline(0.3, 'Color', 'r'); +saveas(gcf, fullfile(outdir, 'trialExample'), 'svg'); +saveas(gcf, fullfile(outdir, 'trialExample'), 'png'); + +% Plot zmin (Figure 2C) +figure('Position', [200, 200, 400, 300]); hold on +plot(stats.zMinMean, 'k-o', 'LineWidth', 1, 'MarkerFaceColor', 'k'); +xlim([1, 6]); ylim([-1.5, 0]); +saveas(gcf, fullfile(outdir, 'trialExampleRmin'), 'svg'); +saveas(gcf, fullfile(outdir, 'trialExampleRmin'), 'png'); + +% Plot variance in increasing order (Figure 2D) +vars = stats.vars(stats.order); +figure('Position', [200, 200, 400, 250]); hold on +plot(vars, 'k-o', 'LineWidth', 1, 'MarkerFaceColor', 'k'); +xlim([1, 6]); +saveas(gcf, fullfile(outdir, 'trialExampleVar'), 'svg'); +saveas(gcf, fullfile(outdir, 'trialExampleVar'), 'png'); + +%% Plot CAR-rereferenced signal at a few sizes + +% how many channels to use for common average to illustrate re-referenced signals (Figure 2B) +ns = [2, 4, 5]; + +for ii = 1:length(ns) + + n = ns(ii); + + CARcurr = mean(V4sorted(:, 1:n), 2); + V4Carn = V4sorted - mean(V4sorted(:, 1:n), 2); + + figure('Position', [200, 200, 150, 500]); ieeg_plotTrials(tt, V4Carn, 100, [], [], 'LineWidth', 1.5); + xlabel('Time (s)'); ylabel('Channels'); + xlim([0.01, 0.3]); + saveas(gcf, fullfile(outdir, sprintf('trialExampleCAR%dch', n)), 'svg'); + saveas(gcf, fullfile(outdir, sprintf('trialExampleCAR%dch', n)), 'png'); + +end + diff --git a/fig3_EPConstruct.m b/fig3_EPConstruct.m new file mode 100644 index 0000000..030a8a0 --- /dev/null +++ b/fig3_EPConstruct.m @@ -0,0 +1,128 @@ +%% This script plots an example of enveloped sinusoids that add to form a simulated evoked potential (Figure 3) +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +%% Configure signal parameters + +rng(4); + +srate = 4800; +tt = -0.5:1/srate:1-1/srate; + +A = 100; + +tau1 = 0.01 + rand*0.02; % time constant to dampen the faster sinusoid (s1) +tau2 = 0.06 + rand*0.08; +tauOn1 = 0.005; % time constant of subtracted exponential decay for both +tauOn2 = 0.025; % time constant of subtracted exponential decay for both + +f1 = 8 + rand*4; % freq (hz) of first sinusoid (mean period = 100 ms) +f2 = 1 + rand*2; % freq of second sinusoid, (mean period = 500 ms) + +ph1 = rand*2*pi; % phase of first sinusoid +ph2 = rand*2*pi; % phase of second sinusoid + +outdir = fullfile('output', 'EPConstruct'); +mkdir(outdir); + +%% Construct and save full signal (3C) + +sig = A*( (exp(-tt/tau1) - exp(-tt/tauOn1)).*sin(2*pi*f1*tt - ph1) + (exp(-tt/tau2) - exp(-tt/tauOn2)).*sin(2*pi*f2*tt - ph2) ); +sig(tt < 0) = 0; % make causal + +figure('Position', [200, 200, 200, 100]); +plot(tt, sig, 'k-'); +xlabel('Time (s)'); +xlim([-0.1, 0.5]); +ylim([-80, 80]); + +saveas(gcf, fullfile(outdir, 'signal'), 'png'); +saveas(gcf, fullfile(outdir, 'signal'), 'svg'); + +%% Plot faster part of signal (3A) + +% the direct/faster signal component +env1 = (exp(-tt/tau1) - exp(-tt/tauOn1)); +sig1 = A*sin(2*pi*f1*tt - ph1); +env1(tt < 0) = 0; % make causal +sig1(tt < 0) = 0; + +figure('Position', [200, 200, 200, 100]); +plot(tt, env1, 'Color', [0.5, 0.5, 0.5], 'LineWidth', 2); +xline(tau1, 'Color', 'r'); xline(tauOn1, 'Color', 'r'); +xlim([-0.1, 0.5]); +ylim([-1, 1]); +saveas(gcf, fullfile(outdir, 'env1'), 'png'); +saveas(gcf, fullfile(outdir, 'env1'), 'svg'); + +figure('Position', [200, 200, 200, 100]); +plot(tt, sig1, 'k-'); +xlim([-0.1, 0.5]); +ylim([-120, 120]); +saveas(gcf, fullfile(outdir, 'sig1'), 'png'); +saveas(gcf, fullfile(outdir, 'sig1'), 'svg'); + +figure('Position', [200, 200, 200, 100]); +plot(tt, env1.*sig1, 'k-'); +xlim([-0.1, 0.5]); +ylim([-80, 80]); +saveas(gcf, fullfile(outdir, 'env1sig1'), 'png'); +saveas(gcf, fullfile(outdir, 'env1sig1'), 'svg'); + +%% Plot slower part of signal (3B) + +env2 = (exp(-tt/tau2) - exp(-tt/tauOn2)); +sig2 = A*sin(2*pi*f2*tt - ph2); +env2(tt < 0) = 0; % make causal +sig2(tt < 0) = 0; + +figure('Position', [200, 200, 200, 100]); +plot(tt, env2, 'Color', [0.5, 0.5, 0.5], 'LineWidth', 2); +xlim([-0.1, 0.5]); +ylim([-1, 1]); +xline(tau2, 'Color', 'r'); xline(tauOn2, 'Color', 'r'); +saveas(gcf, fullfile(outdir, 'env2'), 'png'); +saveas(gcf, fullfile(outdir, 'env2'), 'svg'); + +figure('Position', [200, 200, 200, 100]); +plot(tt, sig2, 'k-'); +xlim([-0.1, 0.5]); +ylim([-120, 120]); +saveas(gcf, fullfile(outdir, 'sig2'), 'png'); +saveas(gcf, fullfile(outdir, 'sig2'), 'svg'); + +figure('Position', [200, 200, 200, 100]); +plot(tt, env2.*sig2, 'k-'); +xlim([-0.1, 0.5]); +ylim([-80, 80]); +saveas(gcf, fullfile(outdir, 'env2sig2'), 'png'); +saveas(gcf, fullfile(outdir, 'env2sig2'), 'svg'); + +%% Save full signal (again), using the constructed components + +figure('Position', [200, 200, 200, 100]); +plot(tt, env1.*sig1 + env2.*sig2, 'k-'); +xlabel('Time (s)'); +xlim([-0.1, 0.5]); +ylim([-80, 80]); diff --git a/fig7_summarizeCARLARealCCEPs.m b/fig7_summarizeCARLARealCCEPs.m new file mode 100644 index 0000000..07c8095 --- /dev/null +++ b/fig7_summarizeCARLARealCCEPs.m @@ -0,0 +1,162 @@ +%% This loads the outputs from applyCARLARealCCEPsLoop to summarize the optimal CCEP cutoff size for stim sites +% in each tissue type across the 4 subjects +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% + +% subjects and runs corresponding to each subject to load +subs = {'sub-1', 'sub-2', 'sub-3', 'sub-4'}; +runs = {[1, 2, 3, 5], [12, 13, 14], [2, 3, 4], [1, 2]}; + +for ss = 1:length(subs) + + sub = subs{ss}; + fprintf('Sub: %s\n', sub); + + %% Obtain stim sites across all runs, number of trials + + elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub)); + elecs = ieeg_readtableRmHyphens(elecsPath); + SOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ')); + + sites = cell(0, 2); % name, number of trials + for rr = 1:length(runs{ss}) + run = runs{ss}(rr); + fprintf('%s run %d\n', sub, run); + + eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run)); + annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run)); + + events = ieeg_readtableRmHyphens(eventsPath, 'electrical_stimulation_site', 1); + eventsStatusOrig = events.status; + events(~strcmpi(events.status, 'good'), :) = []; + + sitesRun = groupby(events.electrical_stimulation_site); % all sites in current run + + % read manual annotations on events if it exists + try + annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t'); + annots = strip(annots.status_description); + catch + warning('Unable to load event annotation for events'); + clear annots + end + + if exist('annots', 'var') + annots(~strcmpi(eventsStatusOrig, 'good')) = []; + assert(length(annots) == height(events), 'Error: length of event annotations does not match events table height'); + end + + % Remove stim sites in SOZ + isSoz = any(ismember(split(sitesRun(:, 1), '-', 2), SOZ), 2); + sitesRun(isSoz, :) = []; + + % remove "all" annotated trials and calculate final number of trials + nTrials = cellfun(@length, sitesRun(:, 2)); + if exist('annots', 'var') + for ii = 1:size(sitesRun, 1) + annotsCurr = annots(sitesRun{ii, 2}); + nTrsRm = sum(strcmpi(annotsCurr, 'all')); + nTrials(ii) = nTrials(ii) - nTrsRm; + end + end + + sites = [sites; sitesRun(:, 1), num2cell(nTrials)]; + + end + + % ensure no stim site has less than 8 trials + sites = cell2table(sites, 'VariableNames', {'site', 'trials'}); + sites(sites.trials < 8, :) = []; + + %% Identify Destrieux labels for each stim site + + % subject-specific Freesurfer derivatives folder + FSdir = fullfile(dataPath, 'derivatives', 'freesurfer', sub); + xyzs = ieeg_getPairXyzs(sites.site, elecs); + labs = ieeg_getLabelXyzDestrieux(xyzs, FSdir); % within radius of 3 + + labGroup = repmat({'n/a'}, length(labs), 1); + labGroup(contains(labs, {'Thalamus', 'Hippocampus', 'Amygdala'})) = {'subgray'}; + labGroup(contains(labs, 'White_Matter')) = {'WM'}; + labGroup(startsWith(labs, {'lh', 'rh'})) = {'cortgray'}; + + fprintf('Subcortical Gray: %d\nCortical Gray: %d\nWhite Matter: %d\nn/a: %d\n', ... + sum(strcmp(labGroup, 'subgray')), sum(strcmp(labGroup, 'cortgray')), sum(strcmp(labGroup, 'WM')), sum(strcmp(labGroup, 'n/a'))); + + sites.labs = labs; + sites.labGroup = labGroup; + + %% load the nCAR files and plot distributions + + rng('default'); + + outdir = fullfile('output', ''realCCEPsLoop', sub); + + carSz = zeros(height(sites), 2); + for ii = 1:height(sites) + M = readmatrix(fullfile(outdir, sprintf('%s_%s_nCAR.txt', sub, sites.site{ii}))); + carSz(ii, :) = M; + end + + sites.carSz = carSz(:, 1); + sites.numChs = carSz(:, 2); + sites.pct = 100*carSz(:, 1)./carSz(:, 2); + + % plot overall distribution + figure; histogram(sites.pct, 10); + + % plot distributions by label group and hemisphere + pctByGroup = cell(6, 1); + + % Right hemi + pctByGroup{1} = sites.pct(strcmp(labGroup, 'cortgray') & startsWith(lower(labs), 'r')); + pctByGroup{2} = sites.pct(strcmp(labGroup, 'subgray') & startsWith(lower(labs), 'r')); + pctByGroup{3} = sites.pct(strcmp(labGroup, 'WM') & startsWith(lower(labs), 'r')); + + % Left hemi + pctByGroup{4} = sites.pct(strcmp(labGroup, 'cortgray') & startsWith(lower(labs), 'l')); + pctByGroup{5} = sites.pct(strcmp(labGroup, 'subgray') & startsWith(lower(labs), 'l')); + pctByGroup{6} = sites.pct(strcmp(labGroup, 'WM') & startsWith(lower(labs), 'l')); + + + figure('Position', [200, 200, 250, 300]); hold on + bar(1:3, cellfun(@mean, {[pctByGroup{1}; pctByGroup{4}], [pctByGroup{2}; pctByGroup{5}], [pctByGroup{3}; pctByGroup{6}]}), ... + 'FaceColor', [0.8, 0.8, 0.8]); % mean for stim sites in both hemispheres + hsR = jitterplot(pctByGroup(1:3), brighten([252, 141, 98]/255, -0.5)); + set(hsR, 'Marker', '.', 'MarkerSize', 16); + hsL = jitterplot(pctByGroup(4:6), brighten([102, 194, 165]/255, -0.5)); + set(hsL, 'Marker', '.', 'MarkerSize', 16); + ylim([0, 100]); + set(gca, 'XTick', 1:3, 'XTickLabels', {'cortgray', 'subgray', 'WM'}); + + [p, tbl, stats] = anova1(sites.pct, sites.labGroup, 'off'); + fprintf('ANOVA p = %0.02e\n', p); + + saveas(gcf, fullfile(outdir, 'summaryPctCAR'), 'png'); + saveas(gcf, fullfile(outdir, 'summaryPctCAR'), 'svg'); + + writetable(sites, fullfile(outdir, 'summaryTbl.tsv'), 'Delimiter', '\t', 'FileType', 'text'); + +end diff --git a/fig8_summarizeInterChCorr.m b/fig8_summarizeInterChCorr.m new file mode 100644 index 0000000..e09c8b5 --- /dev/null +++ b/fig8_summarizeInterChCorr.m @@ -0,0 +1,125 @@ +%% This script uses output from applyCARLARealCCEPsLoop to calculate cross-channel rsqs for each adjusted CAR method +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% + +outdir = fullfile('output', 'realCCEPsLoop'); + +% subjects +subs = {'sub-1', 'sub-2', 'sub-3', 'sub-4'}; + +%% Load .mat files and CARLA results + +% sub, stim site, z-transformed interch correlation matrix, CARLA result +sites = cell(0, 4); + +for ss = 1:length(subs) + + sub = subs{ss}; + + mats = glob(fullfile(outdir, sub, sprintf('%s_*-*_interChRsq.mat', sub))); + + for jj = 1:length(mats) + + % note which site is being added currently + eles = split(mats{jj}, '_'); + load(mats{jj}); + M = readmatrix(fullfile(outdir, sub, sprintf('%s_%s_nCAR.txt', sub, eles{2}))); + + % 3rd dim of zs == number of recording channels + 1 (corresponding to no reference) + assert(size(rsqs, 3) == M(2)+1, 'Mismatch between number of channels in CARLA output and correlation output'); + + sites = [sites; {sub}, eles(2), {rsqs}, {M}]; + end +end + +assert(size(sites, 1) == 82, 'Error: expecting 82 total stim sites across the 4 subjects'); + +%% Summarize mean correlation for no re-referencing, 50% covariance, CARLA, and full rereferencing + +func = (@(x) mean(x, 'all', 'omitnan')); + +% each row = stim site. Order of columns = no rereferencing, 50% rereferencing, full rereferencing, CARLA +rsqsAll = nan(size(sites, 1), 5); +for ii = 1:size(sites, 1) + rsqs = sites{ii, 3}; + nChs = sites{ii, 4}(2); % number of recording channels + nCarla = sites{ii, 4}(1); % optimal carla size + + rsqsAll(ii, 1) = func(rsqs(:, :, 1)); % mean correlation without referencing + rsqsAll(ii, 2) = func(rsqs(:, :, end)); % mean correlation at full CAR + rsqsAll(ii, 3) = func(rsqs(:, :, floor(nChs*0.25) + 1)); % mean correlation at 25% + rsqsAll(ii, 4) = func(rsqs(:, :, floor(nChs*0.5) + 1)); % mean correlation at 50% + rsqsAll(ii, 5) = func(rsqs(:, :, nCarla + 1)); % mean correlation using CARLA optimum +end + +figure('Position', [200, 200, 200, 300]); +boxplot(rsqsAll, 'Symbol', '.'); +xlim([0, 6]); ylim([0, 1]); +saveas(gcf, fullfile(outdir, 'rsCompsAllsubs'), 'png'); +saveas(gcf, fullfile(outdir, 'rsCompsAllsubs'), 'svg'); + +ps = nan(4, 1); +for ii = 1:4 + ps(ii) = signrank(rsqsAll(:, ii), rsqsAll(:, 5)); +end +psBonf = ps*4; +disp(psBonf); + +%% Plot histogram of pair-wise stim site differences between each reference condition and CARLA + +rsqsDiff = rsqsAll - rsqsAll(:, 5); +figure('Position', [200, 200, 200, 300]); +boxplot(rsqsDiff(:, 1:4), 'Symbol', '.'); +yline(0, 'Color', [0.2, 0.2, 0.2]); +xlim([0, 5]); ylim([-0.4, 0.8]); +saveas(gcf, fullfile(outdir, 'rsDiffAllsubs'), 'png'); +saveas(gcf, fullfile(outdir, 'rsDiffAllsubs'), 'svg'); + +%% Plot matrices for sub1 stimsite 1 + +sub = 'sub-1'; +site = 'RMO8-RMO9'; + +idx = find(strcmp(sites(:, 1), sub) & strcmp(sites(:, 2), site)); +rmat = sites{idx, 3}; +nChs = sites{idx, 4}(2); +nCAR = sites{idx, 4}(1); + +% no rereferencing, standard CAR, naive low variance CAR, CARLA +labels = {'none', 'CAR', 'bottom25', 'bottom50', 'CARLA'}; +rmatSelect = rmat(:, :, [1, nChs+1, floor(nChs*0.25)+1, floor(nChs*0.5)+1, nCAR+1]); + +for ii = 1:5 + + figure; + imagesc(rmatSelect(:, :, ii), [0, 1]); + axis square; + colormap(hot); + colorbar; + + saveas(gcf, fullfile(outdir, sprintf('rsMap_%s_%s-%s', labels{ii}, sub, site)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('rsMap_%s_%s-%s', labels{ii}, sub, site)), 'svg'); + +end diff --git a/figS1_compareCARvCARLA.m b/figS1_compareCARvCARLA.m new file mode 100644 index 0000000..28dfcbd --- /dev/null +++ b/figS1_compareCARvCARLA.m @@ -0,0 +1,55 @@ +%% Compare CARLA vs CAR re-referenced signals on same set of channels. Run this script as an addendum after running applyCARLARealCCEPs +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% + +%% Sort the CARLA output and calculate the standard CAR-referenced signals + +% sort the CARLA output channels by CARLA order +VcarlaSorted = Vcar(stats.order, :, :); + +% calculate the standard common average and subtract from signals +CARstandard = mean(dataSiteNoNan, 1, 'omitnan'); +VcarStandard = dataSiteNoNan - CARstandard; +VcarStandardSorted = VcarStandard(stats.order, :, :); + +%% Zoom into the first 50 channels, compare side by side CAR vs CARLA + +yspace = 50; +numChsToPlot = 50; + +figure('Position', [200, 200, 1000, 1400]); + +% the CARLA channels +subplot(1, 2, 1); +ieeg_plotTrials(tt, mean(VcarlaSorted(1:numChsToPlot, :, :), 3)', yspace, 1:numChsToPlot, [0, 0, 0]); +xlim([-0.1, 0.5]); + +% the standard CAR channels +subplot(1, 2, 2); +ieeg_plotTrials(tt, mean(VcarStandardSorted(1:numChsToPlot, :, :), 3)', yspace, 1:numChsToPlot, [0, 0, 0]); +xlim([-0.1, 0.5]); + +saveas(gcf, fullfile(outdir, sprintf('%s_%s_sortedChsCARvCARLA_1to%d', sub, site, numChsToPlot)), 'png'); +print(gcf, fullfile(outdir, sprintf('%s_%s_sortedChsCARvCARLA_1to%d', sub, site, numChsToPlot)),'-depsc2', '-r300', '-painters'); diff --git a/figS3_simCCEPComponentsGlobalSig.m b/figS3_simCCEPComponentsGlobalSig.m new file mode 100644 index 0000000..69ce024 --- /dev/null +++ b/figS3_simCCEPComponentsGlobalSig.m @@ -0,0 +1,191 @@ +%% Figure S2. Create simulated CCEP data and their components, like Figure 1, but this time containing a global signal +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% + +srate = 4800; +tt = -0.5:1/srate:1-1/srate; + +nchs = 10; % number of simulated channels +ntrs = 12; % number of trials +chs = arrayfun(@(x) sprintf('ch%d', x), 1:nchs, 'UniformOutput', false)'; + +outdir = fullfile('output', 'simComponentsGlobal'); +mkdir(outdir) + +%% i) Create data, add stim artifacts + +% stores the data +V0 = zeros(length(tt), nchs); + +Aart = 50 + rand(nchs, 1)*6; % slightly different artifact amplitudes for each channel +artifact = sin(2*pi*600*tt)'; +artifact(tt < 0 | tt > 0.002) = 0; +V0 = V0 + artifact*Aart'; + +figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V0, 100); +xlabel('Time (s)'); ylabel('Channels'); + +%% A) Add unique input signals at a subset of channels + +rng(14); + +chsResp = 1:4; + +A = 100; +% using previous version of signal generator for this illustration to preserve backwards compatibility +sig = genRandSigOrig(tt, length(chsResp), A); + +V1 = V0; +V1(:, chsResp) = V0(:, chsResp) + sig; + +figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V1, 100, [], [], 'LineWidth', 1.5); +xlabel('Time (s)'); ylabel('Channels'); + +%% B) Add a global noise across all chs and trials, e.g. from reference. + +rng(10); + +Aglobal = 20; % same amplitude of noise on all channels + +sigGlob = genRandSigOrig(tt, 1, Aglobal); + +V2 = V1 + sigGlob; + +figure('Position', [200, 200, 400, 800]); ieeg_plotTrials(tt, V2, 100, [], [], 'LineWidth', 1.5); +xlabel('Time (s)'); ylabel('Channels'); + +%% C) Add noise common to all channels (line noise and some brown noise from reference) + +rng('default'); + +% add trials, so now put V in ch x time x trial format +V3 = repmat(V2', 1, 1, ntrs); + +% random line noise phases for each trial, trials x harmonic (60, 120, 180) +phLN = rand(ntrs, 3)*2*pi; +LN = zeros(length(tt), ntrs); +for ii = 1:ntrs + LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3)); +end + +% brown noise shared across channels (from ref electrode +BN = cumsum(0.4*randn(2*length(tt), ntrs)); +BN = ieeg_highpass(BN, srate, true); +BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :); + +noiseCommon = LN + BN; + +V3 = V3 + shiftdim(noiseCommon, -1); + +figure('Position', [200, 200, 600, 400]); +subplot(1, 2, 1); ieeg_plotTrials(tt, V3(:, :, 1)', 100, [], [], 'LineWidth', 1.5); +title('V at one trial'); +xlabel('Time (s)'); ylabel('Channels'); + +subplot(1, 2, 2); ieeg_plotTrials(tt, mean(V3, 3)', 100, [], [], 'LineWidth', 1.5); +title('V across trials'); +xlabel('Time (s)'); ylabel('Channels'); + + +%% D) add random brown noise across all channels + +rng('default'); + +noiseRand = cumsum(0.4*randn(nchs, 2*length(tt), ntrs), 2); % give double the number of time points so we can highpass it +for ii = 1:nchs + noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true); +end +noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :); + +V4 = V3 + noiseRand; + +figure('Position', [200, 200, 600, 400]); +subplot(1, 2, 1); ieeg_plotTrials(tt, V4(:, :, 1)', 100, [], [], 'LineWidth', 1.5); +title('V at one trial'); +xlabel('Time (s)'); ylabel('Channels'); + +subplot(1, 2, 2); ieeg_plotTrials(tt, mean(V4, 3)', 100, [], [], 'LineWidth', 1.5); +title('V across trials'); +xlabel('Time (s)'); ylabel('Channels'); + +saveas(gcf, fullfile(outdir, 'V4'), 'svg'); +saveas(gcf, fullfile(outdir, 'V4'), 'png'); + + +%% Save construction components for channels of interest + +ylims = [-120, 120]; + +% Plot select channels for figure (3 and 8) +chs2Plot = [3, 8]; + +for ii = 1:length(chs2Plot) + + ch = chs2Plot(ii); + + % plot full simulation of all trials at data + figure('Position', [200, 200, 600, 300]); hold on + plot(tt, squeeze(V4(ch, :, :)), 'Color', [0.5, 0.5, 0.5]); + plot(tt, mean(squeeze(V4(ch, :, :)), 2), 'k', 'LineWidth', 1.5); + hold off + ylim(ylims); + xlabel('Time (s)'); ylabel('Voltage (\muV)'); + saveas(gcf, fullfile(outdir, sprintf('V4_ch%d', ch)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('V4_ch%d', ch)), 'svg'); + + % plot the true signal + figure('Position', [200, 200, 200, 300]); + plot(tt, V1(:, ch), 'k', 'LineWidth', 2); + ylim(ylims); xlim([0, 0.5]); + xlabel('Time (s)'); ylabel('Voltage (\muV)'); + saveas(gcf, fullfile(outdir, sprintf('V1_ch%d', ch)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('V1_ch%d', ch)), 'svg'); + + % plot the global noise + figure('Position', [200, 200, 200, 300]); + plot(tt, sigGlob, 'k', 'LineWidth', 2); + ylim(ylims); xlim([0, 0.5]); + xlabel('Time (s)'); ylabel('Voltage (\muV)'); + saveas(gcf, fullfile(outdir, 'globalSig'), 'png'); + saveas(gcf, fullfile(outdir, 'globalNoise'), 'svg'); + + % plot examples of the common noise + figure('Position', [200, 200, 200, 300]); + ieeg_plotTrials(tt, noiseCommon(:, 1:3), 50, [], [0.5, 0.5, 0.5]); % plot at one trial + xlim([0, 0.5]); + xlabel('Time (s)'); ylabel('Voltage (\muV)'); + saveas(gcf, fullfile(outdir, 'lineNoise_3trs'), 'png'); + saveas(gcf, fullfile(outdir, 'lineNoise_3trs'), 'svg'); + + % plot the random noise + figure('Position', [200, 200, 200, 300]); hold on + plot(tt, squeeze(noiseRand(ch, :, :)), 'Color', [0.5, 0.5, 0.5]); + hold off + xlim([0, 0.5]); ylim(ylims); + xlabel('Time (s)'); ylabel('Voltage (\muV)'); + saveas(gcf, fullfile(outdir, sprintf('randNoise_ch%d', ch)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('randNoise_ch%d', ch)), 'svg'); + +end diff --git a/functions/CARLA.m b/functions/CARLA.m new file mode 100644 index 0000000..6fd36b2 --- /dev/null +++ b/functions/CARLA.m @@ -0,0 +1,225 @@ +function [Vout, CAR, stats] = CARLA(tt, V, srate, sens, nboot) +% +% This function performs Common Average Re-referencing by Least Anticorrelation (CARLA). +% Channels are ranked in order of increasing covariance across trials (or variance if single trial). +% Next, channels are iteratively added into a CAR, and the mean correlation between each unre-referenced channel and all re-referenced channels is +% calculated, for each bootstrapped mean signal. zmin denotes the mean correlation belonging, at each bootstrap, to the most anticorrelated channel. +% The optimum CAR size corresponds to when there is least anticorrelation (when zmin takes on its least negative value) +% +% Vout = CARLA(tt, V); [NOT RECOMMENDED] +% Vout = CARLA(tt, V, srate); +% [Vout, stats] = CARLA(tt, V, srate, nboot); +% tt = 1 x t num. Time points matching V, in seconds +% V = n x t x k num, or n x t num. Signal data to be re-referenced, with n channels by t time points by k trials (2 dimensional if k=1). +% srate = num. Sampling frequency of data. If not given, this will be estimated from the time points (not recommended because of potential imprecision). +% sens = (optional) boolean. Determines whether the more sensitive cutoff will be used, corresponding to before the first significant decrease in zMinMean. +% Ignored if V is n x t (needs trials to calculate significance). Default = false (returns global maximum in zMinMean) +% nboot = (optional) integer. Number of bootstrapped mean signal samples to generate when calculating correlations. Default = 100. Ignored if data is n x t +% +% RETURNS: +% Vout = n x t x k num, or n x t num. Re-referenced signal data. +% CAR = t x k num, or 1 x t num. The CAR at each trial +% stats = struct, containing fields: +% chsUsed = m x 1 num, m < n. The sorted n channel numbers used in the optimum CAR +% vars = n x 1 num. Trial covariance of each channel (variance if K = 1) +% order = n x 1 num. The channels sorted in order of increasing covariance (variance if K = 1) +% zMin = n x n x nboot num, or n x n num if k = 1. The z (Fisher z-transformed Pearson's r) for the channel with the most negative +% average z, at each bootstrap. Dimensions are (z to each channel) x (number of channels in CAR subset) x (bootstrap) +% zMinMean = n x nboot num, or 1 x n num if k = 1. Mean zMin across the first dimension (averaged across all "target" rereferenced channels) +% for each bootstrap if V was 3 dimensional +% +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +% 2023/04/27 by Harvey Huang +% + % Ensure ch x time points structure if k = 1 + if ismatrix(V) + if size(V, 1) == length(tt), V = V'; end + end + + if nargin < 5 || isempty(nboot), nboot = 100; end + if nargin < 4 || isempty(sens), sens = false; end + + % Estimate sampling frequency from time points if necessary (not recommended) + if nargin < 3 || isempty(srate) + srate = (length(tt) - 1) / (max(tt) - min(tt)); + warning('sampling frequency was estimated to be %0.02f', srate); + end + + nChs = size(V, 1); + nTrs = size(V, 3); + assert(size(V, 2) == length(tt), 'Error: second dimension does not match tt. Data should be channels x timepoints x trials'); + + stats = struct(); + + % Notch filter signals to remove line noise. This is not used in the final output + Vclean = V; + for ff = 60:60:180 + dNotch = designfilt('bandstopiir', 'FilterOrder', 4, ... + 'DesignMethod', 'butter', ... + 'HalfPowerFrequency1', ff-2, ... + 'HalfPowerFrequency2', ff+2, ... + 'SampleRate', srate); + for ii = 1:nTrs % filter per trial + Vclean(:, :, ii) = filtfilt(dNotch, Vclean(:, :, ii)')'; + end + end + + % extract data on response interval to perform analysis on + twin = [0.01, 0.3]; % should this be an input parameter? For now keep hardcoded + Vseg = Vclean(:, tt >= twin(1) & tt <= twin(2), :); + + % Rank channels. If single trial rank based on increasing variance. If more than 1 trial, rank by covariance (self-self excluded). Better estimator of consistent structure + if nTrs == 1 + stats.vars = var(Vseg, 0, 2); + else + stats.vars = nan(nChs, 1); + for ii = 1:nChs + covCurr = cov(squeeze(Vseg(ii, :, :))); + stats.vars(ii) = mean(covCurr(logical(triu(ones(size(covCurr)), 1))), 'all'); + end + end + [~, stats.order] = sort(stats.vars); + + % channels x CARsize x bootstrapped sample + if nTrs == 1 + stats.zMin = nan(nChs, nChs); % no bootstrapping if only single trial + else + stats.zMin = nan(nChs, nChs, nboot); + end + + % Pull subset U (order(1:ii)) from V, growing U each time by 1, to construct CAR + for ii = 2:nChs % min CAR size is 2 + + % Vseg after subtracting the putative CAR on subset U + VsegReref = Vseg - mean(Vseg(stats.order(1:ii), :, :), 1); + + if nTrs == 1 % no bootstrapping, test correlation at the single trial level + Useg = Vseg(stats.order(1:ii), :)'; % UNREREFERENCED channels in subset to correlate with the rereferenced signals + UsegReref = VsegReref(stats.order(1:ii), :)'; % rereferenced signals in subset + + r = corr(Useg, UsegReref); % rows = which unrereferenced channel was used. + r(1:(ii+1):end) = nan; % omit self-self correlations + z = atanh(r); % fisher z-transform + + % choose "most responsive" channel with greatest anticorrelation, on average + [~, kkMost] = min(mean(z, 2, 'omitnan')); + stats.zMin(stats.order(1:ii), ii) = z(kkMost, :); + + continue; + + end + + % if nTrs > 1 + for jj = 1:nboot + + % bootstrap to calculate mean signal + inds = datasample(1:nTrs, nTrs); + Useg = mean(Vseg(stats.order(1:ii), :, inds), 3)'; + UsegReref = mean(VsegReref(stats.order(1:ii), :, inds), 3)'; + + r = corr(Useg, UsegReref); % rows = which unrereferenced signal was used. + r(1:(ii+1):end) = nan; % omit self-self correlations + z = atanh(r); % fisher z-transform + + % choose "most responsive" channel with greatest anticorrelation, on average, across channels, for this bootstrapped sample + [~, kkMost] = min(mean(z, 2, 'omitnan')); + stats.zMin(stats.order(1:ii), ii, jj) = z(kkMost, :); + + end + + end + + % calculate mean across the unreferenced target channels + stats.zMinMean = mean(stats.zMin, 1, 'omitnan'); + + if sens && nTrs > 1 % find sensitive optimum + + nMin = ceil(0.1*nChs); % minimum number of channels to start at, for stability. Hardcoded at 10 % + zMMxTrs = mean(stats.zMinMean(1, :, :), 3); % mean ZMinMean across bootstrapped samples for each n + + ii = nMin; + while ii <= nChs + + if ii == nChs % at the end. no more comparisons needed + nOptimum = ii; + break + end + + % loop until next local maximum found + if zMMxTrs(ii + 1) > zMMxTrs(ii) + ii = ii + 1; + continue; + end + + zMMxTrs(1:ii-1) = nan; % set all earlier values to nan, to avoid them being found as nextGreater or as trough + + % find index of the next greater zMinmean (ignores earlier nans), if it exists + nextGreater = find(zMMxTrs > zMMxTrs(ii), 1, 'First'); + if isempty(nextGreater) % no more samples greater than current peak. Optimum is identical to global optimum. Done + nOptimum = ii; + break + end + zMMCurr = squeeze(stats.zMinMean(1, ii, :)); % all bootstrapped samples at current max + [~, idx] = min(zMMxTrs(1:nextGreater-1)); % locate Zminmean trough before the next max (ignores earlier nans) + zMMTrough = squeeze(stats.zMinMean(1, idx, :)); % all bootstrapped samples at trough before next greater value + + % all pairwise differences in the mean correlation estimate between the trough and current peak. + % Note that zMinMean is the estimated correlation of the mean signal. We are testing the difference of zMinMean directly, + % Not the mean of zMinMean, which would be redundant and inappropriate because the bootstrapped samples are not independent. + [X, Y] = meshgrid(zMMTrough, zMMCurr); + diffs = X - Y; diffs = diffs(:); + conf = [-inf, prctile(diffs, 95)]; % left-tailed 95% confidence interval (H_a = difference is less than 0) + + % found significant decrease. Current max is optimal. + if conf(2) < 0 + nOptimum = ii; + break + end + + % Difference was not significant. Start at next greater sample + ii = nextGreater; + + end + + if nOptimum == nMin, warning('Optimum CAR was detected at minimal 10% floor'); end + + else + % Find optimum n at max (least negative) zminMean, averaged across trials. + [~, nOptimum] = max(mean(stats.zMinMean, 3)); + end + + % delete first single dimension if 3D, after finding nOptimum + stats.zMinMean = squeeze(stats.zMinMean); + + % which channels were used to make CAR + stats.chsUsed = sort(stats.order(1:nOptimum)); + + CAR = mean(V(stats.chsUsed, :, :)); + + % Using the original (unfiltered signals) + Vout = V - CAR; + + CAR = squeeze(CAR); % squeeze out the first dimension to make t x k, if 3 dimensional + +end diff --git a/functions/antifind.m b/functions/antifind.m new file mode 100644 index 0000000..109b676 --- /dev/null +++ b/functions/antifind.m @@ -0,0 +1,29 @@ +%% Converts numerical indices to logical indices (opposite of built-in find) +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +function idxBool = antifind(idx, len) + idxBool = zeros(len, 1, 'logical'); + idxBool(idx) = true; +end \ No newline at end of file diff --git a/functions/crossChCorr.m b/functions/crossChCorr.m new file mode 100644 index 0000000..fad3698 --- /dev/null +++ b/functions/crossChCorr.m @@ -0,0 +1,105 @@ +function rsqs = crossChCorr(tt, V, srate) +% +% This function calculates how mean cross-channel R-squared change at each CAR size +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% + % Ensure ch x time points structure if k = 1 + if ismatrix(V) + if size(V, 1) == length(tt), V = V'; end + end + + % Estimate sampling frequency from time points if necessary (not recommended) + if nargin < 3 || isempty(srate) + srate = (length(tt) - 1) / (max(tt) - min(tt)); + warning('sampling frequency was estimated to be %0.02f', srate); + end + + nChs = size(V, 1); + nTrs = size(V, 3); + assert(size(V, 2) == length(tt), 'Error: second dimension does not match tt. Data should be channels x timepoints x trials'); + + stats = struct(); + + % Create 60, 120, 180 Hz notch filters + for ff = 60:60:180 + dNotch(ff/60) = designfilt('bandstopiir', 'FilterOrder', 4, ... + 'DesignMethod', 'butter', ... + 'HalfPowerFrequency1', ff-2, ... + 'HalfPowerFrequency2', ff+2, ... + 'SampleRate', srate); + end + + % notch filter signals to use for covariance calculating + Vclean = V; + for ff = 60:60:180 + for ii = 1:nTrs % filter per trial + Vclean(:, :, ii) = filtfilt(dNotch(ff/60), Vclean(:, :, ii)')'; + end + end + + % extract data on response interval to perform analysis on + twin = [0.01, 0.3]; % should this be an input parameter? For now keep hardcoded + Vseg = Vclean(:, tt >= twin(1) & tt <= twin(2), :); + VorigSeg = V(:, tt >= twin(1) & tt <= twin(2), :); + + % Rank channels. If single trial rank based on increasing variance. If more than 1 trial, rank by covariance (self-self excluded). Better estimator of consistent structure + stats.vars = nan(nChs, 1); + for ii = 1:nChs + covCurr = cov(squeeze(Vseg(ii, :, :))); + stats.vars(ii) = mean(covCurr(logical(triu(ones(size(covCurr)), 1))), 'all'); + end + [~, stats.order] = sort(stats.vars); + + % + rsqs = nan(nChs, nChs, nChs+1); % first column corresponds to no rereferencing + + % Pull subset U (order(1:ii)) from V, growing U each time by 1, to construct CAR + for ii = 0:nChs % min CAR size is 2 + + % rereferenced signal is calculated from the non-notch data + if ii == 0 + VsegReref = VorigSeg; + else + VsegReref = VorigSeg - mean(VorigSeg(stats.order(1:ii), :, :), 1); + end + + % apply notch filter to the rereferenced signal before calculating correlations, to reduce effect of line noise + for ff = 60:60:180 + for jj = 1:nTrs % filter per trial + VsegReref(:, :, jj) = filtfilt(dNotch(ff/60), VsegReref(:, :, jj)')'; + end + end + + VsegRerefMean = mean(VsegReref, 3); % mean across trials + + r = corr(VsegRerefMean'); + r(1:(nChs+1):end) = nan; % omit self-self correlations + + rsqs(:, :, ii+1) = r.^2; % convert to R^2 value + + + end + +end diff --git a/functions/genRandSig.m b/functions/genRandSig.m new file mode 100644 index 0000000..c5126c5 --- /dev/null +++ b/functions/genRandSig.m @@ -0,0 +1,52 @@ +%% File to simulate an evoked potential as part of the code to test CARLA +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +function sig = genRandSig(tt, nChs, amp) + % tt = time points + % amp = average amplitude of s1 + % nChs = number of channels to generate + % Adds s1 and s2. s1 represents faster, direct response. s2 represents slower, indirect response + + % signal parameters + A = 0.8*amp + 0.4*amp*rand(nChs, 1); % amplitudes + + tau1 = 0.01 + rand(nChs, 1)*0.02; % time constant to dampen the faster sinusoid (s1) + tau2 = 0.06 + rand(nChs, 1)*0.08; + tauOn1 = 0.005; % time constant of subtracted exponential decay for both + tauOn2 = 0.025; % time constant of subtracted exponential decay for both + + f1 = 8 + rand(nChs, 1)*4; % freq (hz) of first sinusoid (mean period = 100 ms) + f2 = 1 + rand(nChs, 1)*2; % freq of second sinusoid, (mean period = 500 ms) + + ph1 = rand(nChs, 1)*2*pi; % phase of first sinusoid + ph2 = rand(nChs, 1)*2*pi; % phase of second sinusoid + + sig = zeros(length(tt), nChs); + for ii = 1:nChs + sig(:, ii) = A(ii)*( (exp(-tt/tau1(ii)) - exp(-tt/tauOn1)).*sin(2*pi*f1(ii)*tt - ph1(ii)) + (exp(-tt/tau2(ii)) - exp(-tt/tauOn2)).*sin(2*pi*f2(ii)*tt - ph2(ii)) ); + end + sig(tt < 0, :) = 0; % make causal + +end \ No newline at end of file diff --git a/functions/genRandSigOrig.m b/functions/genRandSigOrig.m new file mode 100644 index 0000000..a8603a6 --- /dev/null +++ b/functions/genRandSigOrig.m @@ -0,0 +1,51 @@ +%% Older, alternative version of generating random signals. Used for schematic purposes, not for quantification +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +function sig = genRandSigOrig(tt, nChs, amp) + % tt = time points + % amp = average amplitude of s1 + % nChs = number of channels to generate + % Adds s1 and s2. s1 represents faster, direct response. s2 represents slower, indirect response + + % signal parameters + A = 0.8*amp + 0.4*amp*rand(nChs, 1); % amplitudes + + tau1 = 0.05 + rand(nChs, 1)*0.06; % time constant of slower exp decay + tau2 = ones(nChs, 1)*0.01; % time constant of second exp decay (shorter, fixed) + tau3 = 0.025; % time constant to dampen the nested, faster sinusoid (s1), corresponds to half the half-period of + + f1 = 8 + rand(nChs, 1)*4; % freq (hz) of first sinusoid (mean period = 100 ms) + f2 = 2 + rand(nChs, 1)*2; % freq of second sinusoid, (mean period = 250 ms) + + ph1 = rand(nChs, 1)*2*pi; % phase of first sinusoid + ph2 = rand(nChs, 1)*2*pi; % phase of second sinusoid + + sig = zeros(length(tt), nChs); + for ii = 1:nChs + sig(:, ii) = A(ii)*( (exp(-tt/tau1(ii)) - exp(-tt/tau2(ii))) .* (exp(-tt/tau3).*sin(2*pi*f1(ii)*tt - ph1(ii)) + sin(2*pi*f2(ii)*tt - ph2(ii)))); + end + sig(tt < 0, :) = 0; % make causal + +end \ No newline at end of file diff --git a/functions/getNeighborChs.m b/functions/getNeighborChs.m new file mode 100644 index 0000000..e576fae --- /dev/null +++ b/functions/getNeighborChs.m @@ -0,0 +1,57 @@ +%% Returns the neighbors within n contacts for input channels +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +function neighbors = getNeighborChs(chs, n, selfincl) + if nargin < 3, selfincl = true; end % whether or not the input channel(s) should be included as neighbors + if nargin < 2 || isempty(n), n = 1; end + + isDigit = @(x) x > 47 & x < 58; % returns true for char array elements that are digits (0 - 9) + + if ischar(chs), chs = {chs}; end % if only 1 channel is requested + + neighbors = {}; + for ii = 1:length(chs) + ch = strip(upper(chs{ii})); + + % split channel name into lead and contact + lead = ch(~isDigit(ch)); + contact = str2double(ch(isDigit(ch))); + + % all neighoring contact integers within n + neighborcontacts = contact + (-n:n); + neighborcontacts(neighborcontacts < 1 | neighborcontacts > 18) = []; + + neighbors = [neighbors, arrayfun(@(x) sprintf('%s%d', lead, x), neighborcontacts, 'UniformOutput', false)]; + end + + % remove redundant entries + neighbors = unique(neighbors, 'stable'); + + % remove input channels if necessary + if ~selfincl + neighbors = setdiff(neighbors, chs); + end + +end \ No newline at end of file diff --git a/functions/rmBadTrialsAnnots.m b/functions/rmBadTrialsAnnots.m new file mode 100644 index 0000000..ad50720 --- /dev/null +++ b/functions/rmBadTrialsAnnots.m @@ -0,0 +1,57 @@ +function dataOut = rmBadTrialsAnnots(dataIn, chs, annots) +%% Removes bad trials according to text in a cell array of chars, annotations +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% + assert(size(dataIn, 1) == length(chs), 'Error: First dimension of input data does not match channel names'); + assert(size(dataIn, 3) == length(annots), 'Error: Third dimension of input data does not match annotations length'); + + dataOut = dataIn; + + annots = upper(strip(annots)); % clean up, make uppercase + + % Remove trials that are bad across all channels + trs2Rm = strcmp(annots, 'ALL'); + fprintf('Removing %d/%d trials\n', sum(trs2Rm), length(annots)); + dataOut(:, :, trs2Rm) = []; + annots(trs2Rm) = []; + + % find channels bad in individual trials and set as nan, to be removed later + for ii = 1:length(annots) + try + annot = annots{ii}; + + % Channels mentioned in description. matches strs starting with L/R, followed by 1 or more A-Z letters, then 1 or more numerical digits + out = regexp(annot, '[LR][A-Z]+\d+', 'match'); + if isempty(out), continue; end + + % identify which channels are bad at this trial + idxbad = ismember(chs, out); + fprintf('Setting %d bad channels at trial number %d to NaN values: %s\n', sum(idxbad), ii, string(join(out, ', '))); + dataOut(idxbad, :, ii) = nan; + catch + end + end + +end \ No newline at end of file diff --git a/main.m b/main.m new file mode 100644 index 0000000..5253c9a --- /dev/null +++ b/main.m @@ -0,0 +1,179 @@ +%% Use each section of this main script to output all results and figures for the manuscript. +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package: main script +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Unless specified, each section of code below can be run directly from this script to generate all outputs. +% The scripts listed may also be opened so that individual sections within can be separately examined and run. +% Note that in output filenames, "zminmean" refers to zeta in the manuscript +% +% Subject CCEP and anatomical data will be available in BIDS format upon manuscript publication. +% In the meantime, all results on simulated data can be reproduced with the appropriate section below. +% +%% Configure paths and dependencies + +% Make sure you are currently working in the root manuscript directory (CARLA) +addpath('functions'); + +% Data should be downloaded from OpenNeuro and copied to a folder here (./data) +dataPath = 'data'; + +% Add paths to other dependencies. Change the dummy paths listed here to where they are located on your system. +addpath(genpath('path/to/mnl_ieegBasics')); +addpath(genpath('path/to/matmef')); +addpath('path/to/mnl_seegview'); +addpath(genpath('path/to/vistasoft')); +addpath('path/to/spm12'); addpath(fullfile('path/to/spm12', 'config')); +addpath(genpath('path/to/freesurfer/7.1.1/matlab')); +addpath(genpath('path/to/gifti')); + +% To save figures as vectorized objects +set(0, 'DefaultFigureRenderer', 'painters'); + +%% i) Generate T1 MRI slices with stimulation site electrodes overlaid (Figures 4A, 6C) + +% This was done using sliceGUI.m in the mnl_seegview repository: +% https://github.com/MultimodalNeuroimagingLab/mnl_seegview + +% Define paths to defaced T1 MRI and matching electrodes table +niiPath = fullfile(dataPath, 'derivatives', 'freesurfer', 'sub-1', 'sub-1_ses-mri01_T1w_acpc_deFaced.nii'); +electrodes = fullfile(dataPath, 'sub-1', 'ses-ieeg01', 'ieeg', 'sub-1_ses-ieeg01_electrodes.tsv'); + +% Run sliceGUI, and manually interact with it +sliceGUI(niiPath, electrodes); + +%% Table 1: Calculate the number of measurement electrodes (channels) used in each subject + +% outputs are directly printed to command window +tab1_dispNumElectrodes + +%% Figure 1: Save simulated CCEP channels and their theoretical components + +% outputs saved to .\output\simComponents +fig1_simCCEPComponents + +%% Figure 2: Methods illustration of how CARLA works, on single-trial simulated CCEP data + +% outputs saved to .\output\simulationTrial +fig2_CARLAonSingleTrial + +%% Figure 3: Save example of enveloped sinusoids used to simulate evoked potentials + +% outputs saved to .\output\EPConstruct +fig3_EPConstruct + +%% Figure 4: Simulate 50 channels with responsiveness varying from 0 to 50, with 30 repetitions each, and apply CARLA + +% First run this to generate outputs for each number of responsive channels +% outputs saved to .\output\simLoop +Aglobal = 0; % no global signal desired (0 amplitude) +simCCEPLoop; + +% Then, keeping Aglobal in the workspace, run this to make summary figures (4C, D), as well as identify median candidates for 4A, B +% .\output\simLoop will be renamed .\output\simLoopNoGlob, where outputs are saved. +summarizeSimCCEPLoop + +%% Figure 5: Apply CARLA to real data (sub 1, stim site 1) + +% Configure the subject, data run, and stim site +sub = 'sub-1'; +run = 2; +site = 'RMO8-RMO9'; + +% which channels to plot in panel F +chs2Plot = [100, 203]; + +% outputs saved to .\output\realCCEPs\. These correspond to figures 5B-F. To generate 5A (stim site overlaid on T1 MRI), see section i), above. +applyCARLARealCCEPs + +%% Figure S1: Average channel responses for sub 1, stim site 1 after applying CAR vs. CARLA + +% This supplemental figure is a direct addendum to Figure 5, showing the mean signal for the first 50 channels after rereferencing by CARLA or CAR +% Run after the previous section, keeping all workspace variables +figS1_compareCARvCARLA + +%% Figure S2: Apply CARLA to real data at another 4 stim sites in subjects 1-4 + +% config parameters for each subject shown in S1 +config = struct('sub', {'sub-1', 'sub-2', 'sub-3', 'sub-4'}, ... + 'run', {1, 12, 3, 2}, ... + 'site', {'ROC11-ROC12', 'RZ4-RZ5', 'LA2-LA3', 'RO3-RO4'}); +chs2Plot = 1; % dummy value, no individual channels are plotted in Figure S1 + +for ii = 1:length(config) + sub = config(ii).sub; + run = config(ii).run; + site = config(ii).site; + + % outputs saved to .\output\realCCEPs\. + % Omission of the middle section of sorted channels (subjects 1-3) in Figure S1 was done manually in illustrator. + applyCARLARealCCEPs +end + +%% Figure S3: Save simulated CCEP channels, this time with a global signal, and their theoretical components + +% outputs saved to .\output\simComponentsGlobal +figS3_simCCEPComponentsGlobalSig + +%% Figure 6A, B Simulate 50 channels with responsiveness from 0 to 50, like Figure 4, but now with a global signal added + +% First run this to generate outputs for each number of responsive channels +% outputs saved to .\output\simLoop +Aglobal = 25; % global signal with mean amplitude 25 +simCCEPLoop; + +% Then, keeping Aglobal in the workspace, run this to make summary figure (6B), and to identify a trial candidate for 6A +% .\output\simLoop will be renamed .\output\simLoopWithGlob, where outputs are saved. +summarizeSimCCEPLoop + +% To generate figure 6C, see section i) on using sliceGUI, above. + +%% Figure 6D-F: Apply CARLA to real data containing global signal (sub 1, stim site 2) + +% Configure the subject, data run, and stim site +sub = 'sub-1'; +run = 3; +site = 'RK5-RK6'; + +% which channels to plot in panel F +chs2Plot = [23, 51, 205]; + +% outputs saved to .\output\realCCEPs\ +applyCARLARealCCEPs + +%% Figure 7: Loop through stimulation sites in each subject and compute optimal CARLA cutoff +% This also calculates the cross-channel R-squared values for all possible adjusted common average sizes, to be summarized in Figure 8 + +% First run this script. Outputs are saved to .\output\realCCEPsLoop\. +applyCARLARealCCEPsLoop + +% Then run this script to generate outputs for Figure 7 +fig7_summarizeCARLARealCCEPs + +%% Figure 8: Calculate cross-channel R-squared values for no CAR vs different types of fixed CARs vs CARLA + +% outputs saved to .\output\realCCEPsLoop\. +fig8_summarizeInterChCorr + diff --git a/simCCEPLoop.m b/simCCEPLoop.m new file mode 100644 index 0000000..adffd52 --- /dev/null +++ b/simCCEPLoop.m @@ -0,0 +1,175 @@ +%% This script creates N channels with a certain number of responsive channels. +% At each loop it tests results and accuracy of CARLA. +% Set Aglobal to 0 for outputs to Figure 4 and set Aglobal to 25 for outputs to Figure 6A,B +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +%% First control the parameters of the simulated channels + +srate = 4800; +tt = -0.5:1/srate:1-1/srate; + +nchs = 50; % number of simulated channels +ntrs = 12; % number of trials + +cmSens = [1, 165/255, 0]; % use orange as color for the sensitive optimum + +%% Loop through number of responsive channels + +rng('default'); + +for nresp = 0:1:50 + + %% Create output dir for current responsiveness + + fprintf('%d of %d channels responsive\n', nresp, nchs); + + outdir = fullfile('output', 'simLoop', sprintf('nchs%d-%d', nchs, nresp)); + mkdir(outdir); + + chsResp = 1:nresp; % first x channels responsive, for ease. order doesn't matter. will sort and color them red + + reps = 30; + + %% + for rr = 1:reps % multiple repetitions at each responsive level + + fprintf('.'); % loading bar + + % i) artifact only + V0 = zeros(length(tt), nchs); + Aart = 50 + rand(nchs, 1)*5; % slightly different artifact amplitudes for each channel + artifact = sin(2*pi*600*tt)'; + artifact(tt < 0 | tt > 0.002) = 0; + V0 = V0 + artifact*Aart'; + + % A) Add the evoked potentials + A = 100; + V1 = V0; + sig = genRandSig(tt, length(chsResp), A); + V1(:, chsResp) = V0(:, chsResp) + sig; + + % B) Option to add a global noise. + % Configure in main script for Aglobal to be 0 (for figure 4) or 25 (for figure 6A-B) + if Aglobal + sigCommon = genRandSig(tt, 1, Aglobal); + else + sigCommon = 0; % no global signal + end + V2 = V1 + sigCommon; + + % C) Add common noise to all channels at each trial + V3 = repmat(V2', 1, 1, ntrs); % ch x time points x trial + phLN = rand(ntrs, 3)*2*pi; % LN phases + LN = zeros(length(tt), ntrs); + for ii = 1:ntrs + LN(:, ii) = 8*sin(2*pi*60*tt - phLN(ii, 1)) + 2*sin(2*pi*120*tt - phLN(ii, 2)) + 1*sin(2*pi*180*tt - phLN(ii, 3)); + end + + % brown noise component of common noise shared across channels + BN = cumsum(0.4*randn(2*length(tt), ntrs)); + BN = ieeg_highpass(BN, srate, true); + BN = BN((0.5*length(tt)+1) : 1.5*length(tt), :); + + noiseCommon = LN + BN; + V3 = V3 + shiftdim(noiseCommon, -1); + + % D) add random brown noise + noiseRand = cumsum(0.4*randn(nchs, 2*length(tt), ntrs), 2); % give double the number of time points so we can highpass it + for ii = 1:nchs + noiseRand(ii, :, :) = ieeg_highpass(squeeze(noiseRand(ii, :, :)), srate, true); + end + noiseRand = noiseRand(:, (0.5*length(tt)+1) : 1.5*length(tt), :); + V4 = V3 + noiseRand; + + %% Apply CARLA and plot outputs + + [Vout, CAR, stats] = CARLA(tt, V4, srate, true); % get the sensitive output + + % number of channels used for the CAR + nCAR = length(stats.chsUsed); + [~, nCARglob] = max(mean(stats.zMinMean, 2)); % number of channels at global maximum + + % Plot average zmin across trials + figure('Position', [200, 200, 400, 300]); hold on + errorbar(mean(stats.zMinMean, 2), std(stats.zMinMean, 0, 2), 'k-o'); + plot(nCARglob, mean(stats.zMinMean(nCARglob, :), 2), 'b*'); % global max as blue + if nCARglob ~= nCAR; plot(nCAR, mean(stats.zMinMean(nCAR, :), 2), '*', 'Color', cmSens); end + yline(0, 'Color', 'k'); + saveas(gcf, fullfile(outdir, sprintf('zmin_%d-%d_rep%d', nchs, nresp, rr)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('zmin_%d-%d_rep%d', nchs, nresp, rr)), 'svg'); + + % Plot the variance of channels, sorted in increasing order + vars = stats.vars(stats.order); + figure('Position', [200, 200, 400, 300]); hold on + plot(vars, 'k-o', 'LineWidth', 1, 'MarkerFaceColor', 'k'); + xline(nCARglob + 0.5, 'Color', 'b'); % global threshold + if nCARglob ~= nCAR; xline(nCAR + 0.5, 'Color', cmSens); end % sensitive threshold + xlim([0, nchs+1]); + saveas(gcf, fullfile(outdir, sprintf('covar_%d-%d_rep%d', nchs, nresp, rr)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('covar_%d-%d_rep%d', nchs, nresp, rr)), 'svg'); + + % Sort and plot channels by increasing covariance, draw line at cutoff + V4MeanSorted = mean(V4(stats.order, :, :), 3); + + % create colormap where responsive channels are red + respBool = antifind(chsResp, nchs); + respBool = respBool(stats.order); % logical array of where responsive channels are + cm = zeros(nchs, 3); + cm(respBool, 1) = 1; % make red + figure('Position', [200, 200, 250, 600]); + yspace = 80; + ys = ieeg_plotTrials(tt, V4MeanSorted', 80, [], cm, 'LineWidth', 1); + yline(ys(nCARglob)-yspace/2, 'Color', 'b', 'LineWidth', 1.5); + if nCARglob ~= nCAR; yline(ys(nCAR)-yspace/2, 'Color', cmSens, 'LineWidth', 1.5); end + xlim([-0.1, 0.5]); set(gca, 'xtick', [0, 0.5]); + xlabel('Time (s)'); ylabel('Channels'); + saveas(gcf, fullfile(outdir, sprintf('chsSorted_%d-%d_rep%d', nchs, nresp, rr)), 'png'); + saveas(gcf, fullfile(outdir, sprintf('chsSorted_%d-%d_rep%d', nchs, nresp, rr)), 'svg'); + + close all; + + % Accuracy values. positive means responsive/excluded from CAR + TP = sum(find(respBool) > nCAR); % responsive channels successfully excluded from CAR (above the cutoff) + TN = sum(find(~respBool) <= nCAR); % NR channels successfully below or at cutoff + FN = sum(find(respBool) <= nCAR); % responsive channels incorrectly included in CAR. *This matters most + FP = sum(find(~respBool) > nCAR); % NR channels incorrectly excluded from CAR + + % same for the global threshold + TPglob = sum(find(respBool) > nCARglob); + TNglob = sum(find(~respBool) <= nCARglob); + FNglob = sum(find(respBool) <= nCARglob); + FPglob = sum(find(~respBool) > nCARglob); + + fid = fopen(fullfile(outdir, sprintf('accuracy_%d-%d_rep%d.txt', nchs, nresp, rr)), 'w'); + fprintf(fid, 'TP\t%d\nTN\t%d\nFN\t%d\nFP\t%d\n', TP, TN, FN, FP); + fprintf(fid, 'TPglob\t%d\nTNglob\t%d\nFNglob\t%d\nFPglob\t%d', TPglob, TNglob, FNglob, FPglob); + fclose(fid); + + end + + fprintf('\n'); + +end + diff --git a/summarizeSimCCEPLoop.m b/summarizeSimCCEPLoop.m new file mode 100644 index 0000000..c3ee41c --- /dev/null +++ b/summarizeSimCCEPLoop.m @@ -0,0 +1,197 @@ +%% This script takes outputs from simCCEPLoop, and evaluates accuracy at each level of responsiveness. +% Generates outputs for Figure 4C,D and 6B +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% + +%% Configure loading directories +% +nchs = 50; % total number of channels + +% rename the output folder depending on whether Aglobal was 0 to avoid overwriting Figure 4 outputs with Figure 6 outputs +if ~Aglobal % no global signal + rootdir = fullfile('output', 'simLoopNoGlob'); +else % yes global signal + rootdir = fullfile('output', 'simLoopWithGlob'); +end +movefile(fullfile('output', 'simLoop'), rootdir); + +%% Read and store accuracy values at each subset of responsiveness + +% each field is reps x nresp (starting at 0). XXglob = at cutoff found using global maximum +accuracy = struct(); +accuracy.TP = nan(30, nchs+1); +accuracy.TN = nan(30, nchs+1); +accuracy.FN = nan(30, nchs+1); +accuracy.FP = nan(30, nchs+1); +accuracy.TPglob = nan(30, nchs+1); +accuracy.TNglob = nan(30, nchs+1); +accuracy.FNglob = nan(30, nchs+1); +accuracy.FPglob = nan(30, nchs+1); + +for nresp = 0:nchs + + fprintf('.'); + + indir = fullfile(rootdir, sprintf('nchs%d-%d', nchs, nresp)); + + for rep = 1:30 + + T_acc = readtable(fullfile(indir, sprintf('accuracy_%d-%d_rep%d.txt', nchs, nresp, rep)), 'FileType', 'text', 'Delimiter', '\t'); + accuracy.TP(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'TP')); + accuracy.TN(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'TN')); + accuracy.FN(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'FN')); + accuracy.FP(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'FP')); + accuracy.TPglob(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'TPglob')); + accuracy.TNglob(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'TNglob')); + accuracy.FNglob(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'FNglob')); + accuracy.FPglob(rep, nresp+1) = T_acc.Var2(strcmp(T_acc.Var1, 'FPglob')); + + end + +end +fprintf('\n'); + +%% Identify examples that correspond to median values at each # responsive channels + +FNmed = median(accuracy.FN); +FPmed = median(accuracy.FP); +FNglobMed = median(accuracy.FNglob); +FPglobMed = median(accuracy.FPglob); + +if contains(rootdir, 'simLoopWithGlob') + idxes = 26; % 50% for main figure with global signal + FNmedSelect = FNmed(idxes); + FPmedSelect = FPmed(idxes); + FNglobmedSelect = FNglobMed(idxes); + FPglobmedSelect = FPglobMed(idxes); + fprintf('Median FN at 25 resp chs = %0.1f\n', FNmedSelect); + fprintf('Median FP at 25 = %0.1f\n', FPmedSelect); + fprintf('Median FN (global) at 25 = %0.1f\n', FNglobmedSelect); + fprintf('Median FP (global) at 25 = %0.1f\n', FPglobmedSelect); + +elseif contains(rootdir, 'simLoopNoGlob') + idxes = [1, 11, 21, 31, 41, 46]; % every 20%, for main figure with no global signal + FNmedSelect = FNmed(idxes); + FPmedSelect = FPmed(idxes); + FNglobmedSelect = FNglobMed(idxes); + FPglobmedSelect = FPglobMed(idxes); + fprintf('Median FN at 0, 10, 20, 30, 40, 45 resp chs = %0.1f, %0.1f, %0.1f, %0.1f, %0.1f, %0.1f\n', FNmedSelect(:)); + fprintf('Median FP at 0, 10, 20, 30, 40, 45 = %0.1f, %0.1f, %0.1f, %0.1f, %0.1f, %0.1f\n', FPmedSelect(:)); + fprintf('Median FN (global) at 0, 10, 20, 30, 40, 45 = %0.1f, %0.1f, %0.1f, %0.1f, %0.1f, %0.1f\n', FNglobmedSelect(:)); + fprintf('Median FP (global) at 0, 10, 20, 30, 40, 45 = %0.1f, %0.1f, %0.1f, %0.1f, %0.1f, %0.1f\n', FPglobmedSelect(:)); + +end + +% which rep to use that best matches the median +repSelect = zeros(size(idxes)); +for ii = 1:length(idxes) + rep = find(accuracy.FN(:, idxes(ii)) == floor(FNmedSelect(ii)) & ... + accuracy.FP(:, idxes(ii)) == floor(FPmedSelect(ii)) & ... + accuracy.FNglob(:, idxes(ii)) == floor(FNglobmedSelect(ii)) & ... + accuracy.FPglob(:, idxes(ii)) == floor(FPglobmedSelect(ii)), 1, 'first'); + if isempty(rep), rep = nan; end + repSelect(ii) = rep; +end +fprintf('Rep to use at responsiveness: %d \n', repSelect(:)); + +%% Absolute values of FP, FN + +cmSens = [1, 165/255, 0]; % orange = sensitive threshold + +% Using global optimum +% up to 90% responsiveness (45/50 responsive channels) +figure('Position', [200, 200, 600, 600]); +subplot(2, 1, 1); % FN = responsive channels included in the CAR +boxplot(accuracy.FNglob(:, 1:46), 'Positions', 0:45, 'Colors', 'b', 'PlotStyle', 'compact'); +xlim([-1, 46]); ylim([-1, 50]); +set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off'); +ylabel('Number of incorrectly assigned channels'); +title('False Negative (Resp. Chs in CAR)'); + +subplot(2, 1, 2); +boxplot(accuracy.FPglob(:, 1:46), 'Positions', 0:45, 'Colors', 'b', 'PlotStyle', 'compact'); +xlim([-1, 46]); ylim([-1, 50]); +set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off'); +xlabel('Number of responsive channels / 50'); +ylabel('Number of incorrectly assigned channels'); +title('False Positive (Non-Resp. Chs in CAR)'); +saveas(gcf, fullfile(rootdir, sprintf('nchs%d_FNFPglob', nchs)), 'png'); +saveas(gcf, fullfile(rootdir, sprintf('nchs%d_FNFPglob', nchs)), 'svg'); + +% For sensitive method +% up to 90% responsiveness (45/50 responsive channels) +figure('Position', [200, 200, 600, 600]); +subplot(2, 1, 1); % FN = responsive channels included in the CAR +boxplot(accuracy.FN(:, 1:46), 'Positions', 0:45, 'Colors', cmSens, 'PlotStyle', 'compact'); +xlim([-1, 46]); ylim([-1, 50]); +set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off'); +ylabel('Number of incorrectly assigned channels'); +title('False Negative (Resp. Chs in CAR)'); + +subplot(2, 1, 2); +boxplot(accuracy.FP(:, 1:46), 'Positions', 0:45, 'Colors', cmSens, 'PlotStyle', 'compact'); +xlim([-1, 46]); ylim([-1, 50]); +set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50, 'box', 'off'); +xlabel('Number of responsive channels / 50'); +ylabel('Number of incorrectly assigned channels'); +title('False Positive (Non-Resp. Chs in CAR)'); +saveas(gcf, fullfile(rootdir, sprintf('nchs%d_FNFP', nchs)), 'png'); +saveas(gcf, fullfile(rootdir, sprintf('nchs%d_FNFP', nchs)), 'svg'); + +%% Sensitivity specificity plot + +% sensitivity and specificity +sensglob = accuracy.TPglob ./ (accuracy.TPglob + accuracy.FNglob); +specglob = accuracy.TNglob ./ (accuracy.TNglob + accuracy.FPglob); +sens = accuracy.TP ./ (accuracy.TP + accuracy.FN); % more important. need to sensitively detect responsive channels +spec = accuracy.TN ./ (accuracy.TN + accuracy.FP); + +% balanced accuracy +accglob = (sensglob + specglob)/2; +acc = (sens + spec)/2; + + +% Using global optimum +cmBlue = [0, 0, 139; 115, 147, 179]/255; +figure('Position', [200, 200, 600, 300]); hold on +errorbar(0:45, mean(sensglob(:, 1:46)), std(sensglob(:, 1:46)), '-o', 'MarkerSize', 5, 'MarkerFaceColor', cmBlue(1, :), 'Color', cmBlue(1, :)); +errorbar(0:45, mean(specglob(:, 1:46)), std(specglob(:, 1:46)), '--s', 'MarkerSize', 6, 'Color', cmBlue(2, :)); +xlim([-1, 46]); ylim([-0.2, 1.2]); +set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50); +xlabel('Number of responsive channels / 50'); +saveas(gcf, fullfile(rootdir, sprintf('nchs%d_sensSpecglob', nchs)), 'png'); +saveas(gcf, fullfile(rootdir, sprintf('nchs%d_sensSpecglob', nchs)), 'svg'); + +% Using sensitive method +cmOrange = [139, 64, 0; 255, 172, 28]/255; +figure('Position', [200, 200, 600, 300]); hold on +errorbar(0:45, mean(sens(:, 1:46)), std(sens(:, 1:46)), '-o', 'MarkerSize', 5, 'MarkerFaceColor', cmOrange(1, :), 'Color', cmOrange(1, :)); +errorbar(0:45, mean(spec(:, 1:46)), std(spec(:, 1:46)), '--s', 'MarkerSize', 6, 'Color', cmOrange(2, :)); +xlim([-1, 46]); ylim([-0.2, 1.2]); +set(gca, 'xtick', 0:5:50, 'xticklabels', 0:5:50); +xlabel('Number of responsive channels / 50'); +saveas(gcf, fullfile(rootdir, sprintf('nchs%d_sensSpec', nchs)), 'png'); +saveas(gcf, fullfile(rootdir, sprintf('nchs%d_sensSpec', nchs)), 'svg'); + diff --git a/tab1_dispNumElectrodes.m b/tab1_dispNumElectrodes.m new file mode 100644 index 0000000..f04c50e --- /dev/null +++ b/tab1_dispNumElectrodes.m @@ -0,0 +1,55 @@ +%% This script calculates and prints the number of channels used in the analysis for each subject ("good" in at least 1 run of data) +% +% 2023/09/27 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% 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 . +% +% path to subject data +dataPath = 'data'; + +% subjects and runs corresponding to each subject to load +subs = {'sub-1', 'sub-2', 'sub-3', 'sub-4'}; +runs = {[1, 2, 3, 5], [12, 13, 14], [2, 3, 4], [1, 2]}; + +% iterate through each subject, then each run, and calculate number of good channels +for ss = 1:4 + + sub = subs{ss}; + + seegChStatus = {}; + + % iterate across runs for each subject + for rr = 1:length(runs{ss}) + + run = runs{ss}(rr); + + channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run)); + channels = ieeg_readtableRmHyphens(channelsPath); + seegChStatus = [seegChStatus, channels.status(strcmp(channels.type, 'SEEG'))]; + + end + + % Good channels are defined as any that are not bad across all runs (and therefore not used in any analysis) + goodChs = any(strcmp(seegChStatus, 'good'), 2); + fprintf('Number of good sEEG channels across all runs in %s = %d\n', sub, sum(goodChs)); + +end