OSCR

Clinical Improvement and Taxonomic-Functional Gut Microbiome Remodeling After Six Months of Multi-Strain Synbiotic Supplementation in Mexican Children with Autism Spectrum Disorder.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § 2. Materials and Methods › 2.6. Bioinformatic and Statistical Analysis ↔ humann/humann.py, lines 3–36 · score 0.66 · Unified Metabolic, HUMAnN, pipeline, Network, diamond, Metagenomic
  2. [2] § 2. Materials and Methods › 2.6. Bioinformatic and Statistical Analysis ↔ humann/search/prescreen.py, lines 245–386 · score 0.65 · MetaPhlAn, taxonomic profiles, mpa, gene families, chocophlansgb, v4

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,151 lines · 50 KB · MIT · 1 match

  1. #!/usr/bin/env python
  3. """
  4. HUMAnN : HMP Unified Metabolic Analysis Network 3
  6. HUMAnN is a pipeline for efficiently and accurately determining
  7. the coverage and abundance of microbial pathways in a community
  8. from metagenomic data. Sequencing a metagenome typically produces millions
  9. of short DNA/RNA reads.
  11. Dependencies: MetaPhlAn2, ChocoPhlAn, Bowtie2, and ( Diamond or Rapsearch2 or Usearch )
  13. To Run: humann -i <input.fastq> -o <output_dir>
  15. Copyright (c) 2014 Harvard School of Public Health
  17. Permission is hereby granted, free of charge, to any person obtaining a copy
  18. of this software and associated documentation files (the "Software"), to deal
  19. in the Software without restriction, including without limitation the rights
  20. to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  21. copies of the Software, and to permit persons to whom the Software is
  22. furnished to do so, subject to the following conditions:
  24. The above copyright notice and this permission notice shall be included in
  25. all copies or substantial portions of the Software.
  27. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  28. IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  29. FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  30. AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  31. LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  32. OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  33. THE SOFTWARE.
  34. """
  36. import sys
  38. # Try to load one of the humann modules to check the installation
  39. try:
  40. from . import check
  41. except ImportError:
  42. sys.exit("CRITICAL ERROR: Unable to find the HUMAnN python package." +
  43. " Please check your install.")
  45. # Check the python version
  46. check.python_version()
  48. import argparse
  49. import subprocess
  50. import os
  51. import time
  52. import tempfile
  53. import re
  54. import logging
  56. from . import config
  57. from . import store
  58. from . import utilities
  59. from .search import prescreen
  60. from .search import nucleotide
  61. from .search import translated
  62. from .quantify import families
  63. from .quantify import modules
  65. # name global logging instance
  66. logger=logging.getLogger(__name__)
  69. VERSION="4.0.0.alpha.2"
  71. MAX_SIZE_DEMO_INPUT_FILE=35
  73. def parse_arguments(args):
  74. """
  75. Parse the arguments from the user
  76. """
  77. parser = argparse.ArgumentParser(
  78. description= "HUMAnN : HMP Unified Metabolic Analysis Network 3\n",
  79. formatter_class=argparse.RawTextHelpFormatter,
  80. prog="humann")
  82. common_settings=parser.add_argument_group("[0] Common settings")
  84. common_settings.add_argument(
  85. "-i", "--input",
  86. help="input file of type {" +",".join(config.input_format_choices)+ "} \n[REQUIRED]",
  87. metavar="<input.fastq>",
  88. required=True)
  89. common_settings.add_argument(
  90. "-o", "--output",
  91. help="directory to write output files\n[REQUIRED]",
  92. metavar="<output>",
  93. required=True)
  94. common_settings.add_argument(
  95. "--threads",
  96. help="number of threads/processes\n[DEFAULT: " + str(config.threads) + "]",
  97. metavar="<" + str(config.threads) + ">",
  98. type=int,
  99. default=config.threads)
  100. common_settings.add_argument(
  101. "--version",
  102. action="version",
  103. version="%(prog)s v"+VERSION)
  105. workflow_refinement=parser.add_argument_group("[1] Workflow refinement")
  107. workflow_refinement.add_argument(
  108. "-r","--resume",
  109. help="bypass commands if the output files exist\n",
  110. action="store_true",
  111. default=config.resume)
  112. workflow_refinement.add_argument(
  113. "--bypass-nucleotide-index",
  114. help="bypass the nucleotide index step and run on the indexed ChocoPhlAn database\n",
  115. action="store_true",
  116. default=config.bypass_nucleotide_index)
  117. workflow_refinement.add_argument(
  118. "--bypass-nucleotide-search",
  119. help="bypass the nucleotide search steps\n",
  120. action="store_true",
  121. default=config.bypass_nucleotide_search)
  122. workflow_refinement.add_argument(
  123. "--bypass-prescreen",
  124. help="bypass the prescreen step and run on the full ChocoPhlAn database\n",
  125. action="store_true",
  126. default=config.bypass_prescreen)
  127. workflow_refinement.add_argument(
  128. "--bypass-translated-search",
  129. help="bypass the translated search step\n",
  130. action="store_true",
  131. default=config.bypass_translated_search)
  132. workflow_refinement.add_argument(
  133. "--taxonomic-profile",
  134. help="a taxonomic profile (the output file created by metaphlan)\n[DEFAULT: file will be created]",
  135. metavar="<taxonomic_profile.tsv>")
  136. workflow_refinement.add_argument(
  137. "--memory-use",
  138. help="the amount of memory to use\n[DEFAULT: " +
  139. config.memory_use + "]",
  140. default=config.memory_use,
  141. choices=config.memory_use_options)
  142. workflow_refinement.add_argument(
  143. "--input-format",
  144. help="the format of the input file\n[DEFAULT: format identified by software]",
  145. choices=config.input_format_choices)
  146. workflow_refinement.add_argument(
  147. "-v","--verbose",
  148. help="additional output is printed\n",
  149. action="store_true",
  150. default=config.verbose)
  152. tier1_prescreen=parser.add_argument_group("[2] Configure tier 1: prescreen")
  154. tier1_prescreen.add_argument(
  155. "--metaphlan",
  156. help="directory containing the MetaPhlAn software\n[DEFAULT: $PATH]",
  157. metavar="<metaphlan>")
  158. tier1_prescreen.add_argument(
  159. "--metaphlan-options",
  160. help="options to be provided to the MetaPhlAn software\n[DEFAULT: \"" + str(" ".join(config.metaphlan_opts)) + "\"]",
  161. metavar="<metaphlan_options>",
  162. default=" ".join(config.metaphlan_opts))
  163. tier1_prescreen.add_argument(
  164. "--prescreen-threshold",
  165. help="minimum estimated genome coverage for inclusion in pangenome search\n[DEFAULT: "
  166. + str(config.prescreen_threshold) + "]",
  167. metavar="<" + str(config.prescreen_threshold) + ">",
  168. type=float,
  169. default=config.prescreen_threshold)
  170. tier1_prescreen.add_argument(
  171. "--average-read-length",
  172. help="average read length for input file\n[DEFAULT: "
  173. + str(config.average_read_length) + "]",
  174. metavar="<" + str(config.average_read_length) + ">",
  175. type=float,
  176. default=config.average_read_length)
  178. tier2_nucleotide_search=parser.add_argument_group("[3] Configure tier 2: nucleotide search")
  180. tier2_nucleotide_search.add_argument(
  181. "--bowtie2",
  182. help="directory containing the bowtie2 executable\n[DEFAULT: $PATH]",
  183. metavar="<bowtie2>")
  184. tier2_nucleotide_search.add_argument(
  185. "--bowtie-options",
  186. help="options to be provided to the bowtie software\n[DEFAULT: \"" + str(" ".join(config.bowtie2_align_opts)) + "\"]",
  187. metavar="<bowtie_options>",
  188. default=" ".join(config.bowtie2_align_opts))
  189. tier2_nucleotide_search.add_argument(
  190. "--nucleotide-database",
  191. help="directory containing the nucleotide database\n[DEFAULT: "
  192. + config.nucleotide_database + "]",
  193. metavar="<nucleotide_database>")
  194. tier2_nucleotide_search.add_argument(
  195. "--nucleotide-identity-threshold",
  196. help="identity threshold for nucleotide alignments\n[DEFAULT: " + str(config.nucleotide_identity_threshold) + "]",
  197. metavar="<" + str(config.nucleotide_identity_threshold) + ">",
  198. default=config.nucleotide_identity_threshold,
  199. type=float)
  200. tier2_nucleotide_search.add_argument(
  201. "--nucleotide-query-coverage-threshold",
  202. help="query coverage threshold for nucleotide alignments\n[DEFAULT: "
  203. + str(config.nucleotide_query_coverage_threshold) + "]",
  204. metavar="<" + str(config.nucleotide_query_coverage_threshold) + ">",
  205. type=float,
  206. default=config.nucleotide_query_coverage_threshold)
  207. tier2_nucleotide_search.add_argument(
  208. "--nucleotide-subject-coverage-threshold",
  209. help="subject coverage threshold for nucleotide alignments\n[DEFAULT: "
  210. + str(config.nucleotide_subject_coverage_threshold) + "]",
  211. metavar="<" + str(config.nucleotide_subject_coverage_threshold) + ">",
  212. type=float,
  213. default=config.nucleotide_subject_coverage_threshold)
  216. tier3_translated_search=parser.add_argument_group("[3] Configure tier 2: translated search")
  218. tier3_translated_search.add_argument(
  219. "--diamond",
  220. help="directory containing the diamond executable\n[DEFAULT: $PATH]",
  221. metavar="<diamond>")
  222. tier3_translated_search.add_argument(
  223. "--diamond-options",
  224. help="options to be provided to the diamond software\n[DEFAULT: \"" + str(" ".join(config.diamond_opts_uniref50)) + "\"]",
  225. metavar="<diamond_options>",
  226. default=config.diamond_opts_uniref50)
  227. tier3_translated_search.add_argument(
  228. "--evalue",
  229. help="the evalue threshold to use with the translated search\n[DEFAULT: " + str(config.evalue_threshold) + "]",
  230. metavar="<" + str(config.evalue_threshold) + ">",
  231. type=float,
  232. default=config.evalue_threshold)
  233. tier3_translated_search.add_argument(
  234. "--protein-database",
  235. help="directory containing the protein database\n[DEFAULT: "
  236. + config.protein_database + "]",
  237. metavar="<protein_database>")
  238. tier3_translated_search.add_argument(
  239. "--translated-alignment",
  240. help="software to use for translated alignment\n[DEFAULT: " +
  241. config.translated_alignment_selected + "]",
  242. default=config.translated_alignment_selected,
  243. choices=config.translated_alignment_choices)
  244. tier3_translated_search.add_argument(
  245. "--translated-identity-threshold",
  246. help="identity threshold for translated alignments\n[DEFAULT: " + str(config.identity_threshold_uniref50_mode) + "]",
  247. metavar="<" + str(config.identity_threshold_uniref50_mode) + ">",
  248. default=config.identity_threshold_uniref50_mode,
  249. type=float,
  250. dest="identity_threshold")
  251. tier3_translated_search.add_argument(
  252. "--translated-query-coverage-threshold",
  253. help="query coverage threshold for translated alignments\n[DEFAULT: "
  254. + str(config.translated_query_coverage_threshold) + "]",
  255. metavar="<" + str(config.translated_query_coverage_threshold) + ">",
  256. type=float,
  257. default=config.translated_query_coverage_threshold)
  258. tier3_translated_search.add_argument(
  259. "--translated-subject-coverage-threshold",
  260. help="subject coverage threshold for translated alignments\n[DEFAULT: "
  261. + str(config.translated_subject_coverage_threshold) + "]",
  262. metavar="<" + str(config.translated_subject_coverage_threshold) + ">",
  263. type=float,
  264. default=config.translated_subject_coverage_threshold)
  267. gene_and_pathway=parser.add_argument_group("[5] Gene and pathway quantification")
  269. gene_and_pathway.add_argument(
  270. "--count-normalization",
  271. help="normalization mode for results from both nucleotide and translated search\n[DEFAULT: " +
  272. config.count_normalization + "]",
  273. default=config.count_normalization,
  274. choices=config.count_normalization_choices)
  275. gene_and_pathway.add_argument(
  276. "--utility-database",
  277. help="directory containing the utility database\n[DEFAULT: "
  278. + config.utility_mapping_database + "]",
  279. metavar="<utility_mapping_database>")
  280. gene_and_pathway.add_argument(
  281. "--gap-fill",
  282. help="turn on/off the gap fill computation\n[DEFAULT: " +
  283. config.gap_fill_toggle + "]",
  284. default=config.gap_fill_toggle,
  285. choices=config.toggle_choices)
  286. gene_and_pathway.add_argument(
  287. "--minpath",
  288. help="turn on/off the minpath computation\n[DEFAULT: " +
  289. config.minpath_toggle + "]",
  290. default=config.minpath_toggle,
  291. choices=config.toggle_choices)
  292. gene_and_pathway.add_argument(
  293. "--pathways",
  294. help="the database to use for pathway computations\n[DEFAULT: " +
  295. config.pathways_database + "]",
  296. default=config.pathways_database,
  297. choices=config.pathways_database_choices)
  298. gene_and_pathway.add_argument(
  299. "--pathways-database",
  300. help="mapping file (or files, at most two in a comma-delimited list) to use for pathway computations\n[DEFAULT: " +
  301. config.pathways_database + " database ]",
  302. metavar=("<pathways_database.tsv>"),
  303. nargs=1)
  304. gene_and_pathway.add_argument(
  305. "--xipe",
  306. help="turn on/off the xipe computation\n[DEFAULT: " +
  307. config.xipe_toggle + "]",
  308. default=config.xipe_toggle,
  309. choices=config.toggle_choices)
  310. gene_and_pathway.add_argument(
  311. "--annotation-gene-index",
  312. help="the index of the gene in the sequence annotation\n[DEFAULT: "
  313. + ",".join(str(i) for i in config.chocophlan_gene_indexes) + "]",
  314. metavar="<"+",".join(str(i) for i in config.chocophlan_gene_indexes)+">",
  315. default=",".join(str(i) for i in config.chocophlan_gene_indexes))
  316. gene_and_pathway.add_argument(
  317. "--id-mapping",
  318. help="id mapping file for alignments\n[DEFAULT: alignment reference used]",
  319. metavar="<id_mapping.tsv>")
  322. more_output_config=parser.add_argument_group("[6] More output configuration")
  324. more_output_config.add_argument(
  325. "--remove-temp-output",
  326. help="remove temp output files\n" +
  327. "[DEFAULT: temp files are not removed]",
  328. action="store_true")
  329. more_output_config.add_argument(
  330. "--log-level",
  331. help="level of messages to display in log\n" +
  332. "[DEFAULT: " + config.log_level + "]",
  333. default=config.log_level,
  334. choices=config.log_level_choices)
  335. more_output_config.add_argument(
  336. "--o-log",
  337. help="log file\n" +
  338. "[DEFAULT: sample_0.log]",
  339. metavar="<sample_0.log>")
  340. more_output_config.add_argument(
  341. "--output-basename",
  342. help="the basename for the output files\n[DEFAULT: " +
  343. "input file basename]",
  344. default=config.file_basename,
  345. metavar="<sample_name>")
  346. more_output_config.add_argument(
  347. "--output-format",
  348. help="the format of the output files\n[DEFAULT: " +
  349. config.output_format + "]",
  350. default=config.output_format,
  351. choices=config.output_format_choices)
  352. more_output_config.add_argument(
  353. "--output-max-decimals",
  354. help="the number of decimals to output\n[DEFAULT: " +
  355. str(config.output_max_decimals) + "]",
  356. metavar="<" + str(config.output_max_decimals) + ">",
  357. type=int,
  358. default=config.output_max_decimals)
  359. more_output_config.add_argument(
  360. "--remove-column-description-output",
  361. help="remove the description in the output column\n" +
  362. "[DEFAULT: output column includes description]",
  363. action="store_true",
  364. default=config.remove_column_description_output)
  365. more_output_config.add_argument(
  366. "--remove-stratified-output",
  367. help="remove stratification from output\n" +
  368. "[DEFAULT: output is stratified]",
  369. action="store_true",
  370. default=config.remove_stratified_output)
  372. return parser.parse_args()
  374. def update_configuration(args):
  375. """
  376. Update the configuration settings based on the arguments
  377. """
  379. # Use the full path to the input file
  380. args.input=os.path.abspath(args.input)
  382. # set the version header
  383. config.version_header=config.version_header+VERSION
  385. # If set, append paths executable locations
  386. if args.metaphlan:
  387. utilities.add_exe_to_path(os.path.abspath(args.metaphlan))
  389. if args.bowtie2:
  390. utilities.add_exe_to_path(os.path.abspath(args.bowtie2))
  392. if args.diamond:
  393. utilities.add_exe_to_path(os.path.abspath(args.diamond))
  395. # Set the metaphlan options, removing any extra spaces
  396. if args.metaphlan_options:
  397. config.metaphlan_opts=list(filter(None,args.metaphlan_options.split(" ")))
  399. # check for custom diamond options
  400. config.diamond_opts=args.diamond_options
  401. if args.diamond_options and args.diamond_options != config.diamond_opts_uniref50:
  402. config.diamond_options_custom = True
  403. config.diamond_opts=list(filter(None,args.diamond_options.split(" ")))
  405. if args.bowtie_options:
  406. config.bowtie2_align_opts=list(filter(None,args.bowtie_options.split(" ")))
  408. # Set the pathways database selection
  409. if args.pathways == "metacyc":
  410. config.pathways_database_part1=config.metacyc_gene_to_reactions
  411. config.pathways_database_part2=config.metacyc_reactions_to_pathways
  412. config.pathways_ec_column=True
  413. elif args.pathways == "unipathway":
  414. config.pathways_database_part1=config.unipathway_database_part1
  415. config.pathways_database_part2=config.unipathway_database_part2
  416. config.pathways_ec_column=True
  418. # Set the locations of the pathways databases
  419. # If provided by the user, this will take precedence over the pathways database selection
  420. if args.pathways_database:
  421. custom_pathways_files=args.pathways_database[0].split(",")
  422. if len(custom_pathways_files)==2:
  423. config.pathways_database_part1=os.path.abspath(custom_pathways_files[0])
  424. config.pathways_database_part2=os.path.abspath(custom_pathways_files[1])
  425. elif len(custom_pathways_files)==1:
  426. config.pathways_database_part1=None
  427. config.pathways_database_part2=os.path.abspath(custom_pathways_files[0])
  428. else:
  429. sys.exit("ERROR: Please provide one or two pathways files.")
  430. config.pathways_ec_column=False
  432. # Set the locations of the other databases
  433. if args.nucleotide_database:
  434. config.nucleotide_database=os.path.abspath(args.nucleotide_database)
  436. if args.protein_database:
  437. config.protein_database=os.path.abspath(args.protein_database)
  439. # if set, update the config run mode to resume
  440. if args.resume:
  441. config.resume=True
  443. # if set, update the config run mode to verbose
  444. if args.verbose:
  445. config.verbose=True
  447. # if set, update the config run mode to bypass prescreen step
  448. if args.bypass_prescreen:
  449. config.bypass_prescreen=True
  451. # if set, update the config run mode to bypass nucleotide index step
  452. if args.bypass_nucleotide_index:
  453. config.bypass_nucleotide_index=True
  454. config.bypass_prescreen=True
  456. # if set, update the config run mode to bypass translated search step
  457. # set the pick_frames toggle based on the bypass
  458. config.pick_frames_toggle="off"
  459. if args.bypass_translated_search:
  460. config.bypass_translated_search=True
  462. # if set, update the config run mode to bypass nucleotide search steps
  463. if args.bypass_nucleotide_search:
  464. config.bypass_prescreen=True
  465. config.bypass_nucleotide_index=True
  466. config.bypass_nucleotide_search=True
  468. # update the average read length
  469. config.average_read_length=args.average_read_length
  471. # Update thresholds
  472. config.prescreen_threshold=args.prescreen_threshold
  473. config.translated_subject_coverage_threshold=args.translated_subject_coverage_threshold
  474. config.nucleotide_subject_coverage_threshold=args.nucleotide_subject_coverage_threshold
  475. config.translated_query_coverage_threshold=args.translated_query_coverage_threshold
  476. config.nucleotide_query_coverage_threshold=args.nucleotide_query_coverage_threshold
  478. # Update the max decimals output
  479. config.output_max_decimals=args.output_max_decimals
  481. # Update memory use
  482. config.memory_use=args.memory_use
  484. # Update threads
  485. config.threads=args.threads
  487. # Update the evalue threshold
  488. config.evalue_threshold=args.evalue
  490. # Update translated alignment software
  491. config.translated_alignment_selected=args.translated_alignment
  493. # Update the computation toggle choices
  494. config.xipe_toggle=args.xipe
  495. config.minpath_toggle=args.minpath
  496. config.gap_fill_toggle=args.gap_fill
  497. config.count_normalization=args.count_normalization
  499. if args.utility_database:
  500. config.utility_mapping_database=os.path.abspath(args.utility_database)
  502. # Check that the input file exists and is readable
  503. if not os.path.isfile(args.input):
  504. sys.exit("CRITICAL ERROR: Can not find input file selected: "+ args.input)
  506. if not os.access(args.input, os.R_OK):
  507. sys.exit("CRITICAL ERROR: Not able to read input file selected: " + args.input)
  509. # Update the output format
  510. config.remove_stratified_output=args.remove_stratified_output
  511. config.remove_column_description_output=args.remove_column_description_output
  513. # Check that the output directory is writeable
  514. output_dir = os.path.abspath(args.output)
  516. if not os.path.isdir(output_dir):
  517. try:
  518. print("Creating output directory: " + output_dir)
  519. os.mkdir(output_dir)
  520. except EnvironmentError:
  521. sys.exit("CRITICAL ERROR: Unable to create output directory.")
  523. if not os.access(output_dir, os.W_OK):
  524. sys.exit("CRITICAL ERROR: The output directory is not " +
  525. "writeable. This software needs to write files to this directory.\n" +
  526. "Please select another directory.")
  528. print("Output files will be written to: " + output_dir)
  530. # Set the basename of the output files if specified as an option
  531. if args.output_basename:
  532. config.file_basename=args.output_basename
  533. else:
  534. # Determine the basename of the input file to use as output file basename
  535. input_file_basename=os.path.basename(args.input)
  536. # Remove gzip extension if present
  537. if re.search('.gz$',input_file_basename):
  538. input_file_basename='.'.join(input_file_basename.split('.')[:-1])
  539. # Remove input file extension if present
  540. if '.' in input_file_basename:
  541. input_file_basename='.'.join(input_file_basename.split('.')[:-1])
  543. config.file_basename=input_file_basename
  545. # Set the output format
  546. config.output_format=args.output_format
  548. # Set final output file names and location
  549. config.pathabundance_file=os.path.join(output_dir,
  550. config.file_basename + config.pathabundance_file + "." +
  551. config.output_format)
  552. config.pathcoverage_file=os.path.join(output_dir,
  553. config.file_basename + config.pathcoverage_file + "." +
  554. config.output_format)
  555. config.genefamilies_file=os.path.join(output_dir,
  556. config.file_basename + config.genefamilies_file + "." +
  557. config.output_format)
  558. config.profile_file=os.path.join(output_dir,
  559. config.file_basename + config.profile_file + "." +
  560. config.output_format)
  561. config.reactions_file=os.path.join(output_dir,
  562. config.file_basename + config.reactions_file + "." +
  563. config.output_format)
  565. # set the location of the temp directory
  566. if not args.remove_temp_output:
  567. config.temp_dir=os.path.join(output_dir,config.file_basename+"_humann_temp")
  568. if not os.path.isdir(config.temp_dir):
  569. try:
  570. os.mkdir(config.temp_dir)
  571. except EnvironmentError:
  572. sys.exit("Unable to create temp directory: " + config.temp_dir)
  573. else:
  574. config.temp_dir=tempfile.mkdtemp(
  575. prefix=config.file_basename+'_humann_temp_',dir=output_dir)
  577. # create the unnamed temp directory
  578. config.unnamed_temp_dir=tempfile.mkdtemp(dir=config.temp_dir)
  580. # set the name of the log file
  581. log_file=os.path.join(output_dir,config.file_basename+"_0.log")
  583. # change file name if set
  584. if args.o_log:
  585. log_file=args.o_log
  587. # configure the logger
  588. logging.basicConfig(filename=log_file,format='%(asctime)s - %(name)s - %(levelname)s: %(message)s',
  589. level=getattr(logging,args.log_level), filemode='w', datefmt='%m/%d/%Y %I:%M:%S %p')
  591. # write the version of the software to the log
  592. logger.info("Running humann v"+VERSION)
  594. # write the location of the output files to the log
  595. logger.info("Output files will be written to: " + output_dir)
  597. # write the location of the temp file directory to the log
  598. message="Writing temp files to directory: " + config.temp_dir
  599. logger.info(message)
  600. if config.verbose:
  601. print("\n"+message+"\n")
  603. return log_file
  605. def parse_chocophlan_gene_indexes(annotation_gene_index):
  606. """ Parse the chocophlan gene index input """
  608. # Update the chocophlan gene indexes
  609. chocophlan_gene_indexes=[]
  610. for index in annotation_gene_index.split(","):
  611. # Look for array range
  612. if ":" in index:
  613. split_index=index.split(":")
  614. start=split_index[0]
  615. end=split_index[1]
  617. try:
  618. start=int(start)
  619. end=int(end)
  620. chocophlan_gene_indexes+=range(start,end)
  621. except ValueError:
  622. pass
  623. else:
  624. # Convert to int
  625. try:
  626. index=int(index)
  627. chocophlan_gene_indexes.append(index)
  628. except ValueError:
  629. pass
  631. return chocophlan_gene_indexes
  633. def check_requirements(args):
  634. """
  635. Check requirements (file format, dependencies, permissions)
  636. """
  638. # Check the pathways database files exist and are readable
  639. if config.pathways_database_part1:
  640. utilities.file_exists_readable(config.pathways_database_part1)
  641. utilities.file_exists_readable(config.pathways_database_part2)
  643. # Determine the input file format if not provided
  644. if not args.input_format:
  645. args.input_format=utilities.determine_file_format(args.input)
  647. if args.input_format == "unknown":
  648. sys.exit("CRITICAL ERROR: Unable to determine the input file format." +
  649. " Please provide the format with the --input_format argument.")
  651. config.input_format = args.input_format
  653. # If the input file is compressed, then decompress
  654. if args.input_format.endswith(".gz"):
  655. new_file=utilities.gunzip_file(args.input)
  657. if new_file:
  658. args.input=new_file
  659. args.input_format=args.input_format.split(".")[0]
  660. else:
  661. sys.exit("CRITICAL ERROR: Unable to use gzipped input file. " +
  662. " Please check the format of the input file.")
  664. # check if the input file has sequence identifiers of the new illumina casava v1.8+ format
  665. # these have spaces causing the paired end reads to have the same identifier after
  666. # delimiting by space (causing bowtie2/diamond to label the reads with the same identifier)
  667. if args.input_format in ["fasta","fastq"]:
  668. if utilities.space_in_identifier(args.input):
  669. message="Removing spaces from identifiers in input file"
  670. print(message+" ...\n")
  671. logger.info(message)
  672. new_file = utilities.remove_spaces_from_file(args.input)
  674. if new_file:
  675. args.input=new_file
  676. else:
  677. sys.exit("CRITICAL ERROR: Unable to remove spaces from identifiers in input file.")
  679. # If the input format is in binary then convert to sam (tab-delimited text)
  680. if args.input_format == "bam":
  682. # Check for the samtools software
  683. if not utilities.find_exe_in_path("samtools"):
  684. sys.exit("CRITICAL ERROR: The samtools executable can not be found. "
  685. "Please check the install or select another input format.")
  687. new_file=utilities.bam_to_sam(args.input)
  689. if new_file:
  690. args.input=new_file
  691. args.input_format="sam"
  692. else:
  693. sys.exit("CRITICAL ERROR: Unable to convert bam input file to sam.")
  695. # If the input format is in biom then convert to tsv
  696. if args.input_format == "biom":
  698. # Check for the biom software
  699. if not utilities.find_exe_in_path("biom"):
  700. sys.exit("CRITICAL ERROR: The biom executable can not be found. "
  701. "Please check the install or select another input format.")
  703. new_file=utilities.biom_to_tsv(args.input)
  705. if new_file:
  706. args.input=new_file
  707. # determine the format of the file
  708. args.input_format=utilities.determine_file_format(args.input)
  709. else:
  710. sys.exit("CRITICAL ERROR: Unable to convert biom input file to tsv.")
  712. # If the biom output format is selected, check for the biom package
  713. if config.output_format=="biom":
  714. if not utilities.find_exe_in_path("biom"):
  715. sys.exit("CRITICAL ERROR: The biom executable can not be found. "
  716. "Please check the install or select another output format.")
  718. try:
  719. import biom
  720. except ImportError:
  721. sys.exit("Could not find the biom software."+
  722. " This software is required since the output file is a biom file.")
  724. if os.path.basename(config.utility_mapping_database) == "utility_DEMO":
  725. # Check the input file is a demo input if running with demo database
  726. try:
  727. input_file_size=os.path.getsize(args.input)/1024**2
  728. except EnvironmentError:
  729. input_file_size=0
  730. if input_file_size > MAX_SIZE_DEMO_INPUT_FILE:
  731. sys.exit("ERROR: You are using the demo utility database with "
  732. + "a non-demo input file. If you have not already done so, please "
  733. + "run humann_databases to download the full utility database. "
  734. + "If you have downloaded the full database, use the option "
  735. + "--utility-database to provide the location. "
  736. + "You can also run humann_config to update the default "
  737. + "database location. For additional information, please "
  738. + "see the HUMAnN User Manual.")
  740. # If the file is fasta/fastq check for requirements
  741. if args.input_format in ["fasta","fastq"]:
  742. # Check that the chocophlan directory exists
  743. if not config.bypass_nucleotide_index:
  744. if not os.path.isdir(config.nucleotide_database):
  745. if args.nucleotide_database:
  746. sys.exit("CRITICAL ERROR: The directory provided for the ChocoPhlAn database at "
  747. + args.nucleotide_database + " does not exist. Please select another directory.")
  748. else:
  749. sys.exit("CRITICAL ERROR: The default ChocoPhlAn database directory of "
  750. + config.nucleotide_database + " does not exist. Please provide the location "
  751. + "of the ChocoPhlAn directory using the --nucleotide-database option.")
  753. # Check that the files in the chocophlan folder are of the right format
  754. if not config.bypass_nucleotide_index:
  755. valid_format_count=0
  756. for file in os.listdir(config.nucleotide_database):
  757. # expect most of the file names to be of the format g__*s__*
  758. if re.search("^SGB",file):
  759. valid_format_count+=1
  761. if not config.metaphlan_v4_db_matching_uniref in file:
  762. sys.exit("\n\nCRITICAL ERROR: The directory provided for ChocoPhlAn contains files ( "+file+" )"+\
  763. " that are not of the expected version. Please install the latest version"+\
  764. " of the database: "+config.metaphlan_v4_db_matching_uniref)
  766. if valid_format_count == 0:
  767. sys.exit("CRITICAL ERROR: The directory provided for ChocoPhlAn does not "
  768. + "contain files of the expected format (ie \'^SGB\').")
  770. # Check if running with the demo database
  771. if not config.bypass_nucleotide_index:
  772. if os.path.basename(config.nucleotide_database) == "chocophlan_DEMO":
  773. # Check the input file is a demo input if running with demo database
  774. try:
  775. input_file_size=os.path.getsize(args.input)/1024**2
  776. except EnvironmentError:
  777. input_file_size=0
  778. if input_file_size > MAX_SIZE_DEMO_INPUT_FILE:
  779. sys.exit("ERROR: You are using the demo ChocoPhlAn database with "
  780. + "a non-demo input file. If you have not already done so, please "
  781. + "run humann_databases to download the full ChocoPhlAn database. "
  782. + "If you have downloaded the full database, use the option "
  783. + "--nucleotide-database to provide the location. "
  784. + "You can also run humann_config to update the default "
  785. + "database location. For additional information, please "
  786. + "see the HUMAnN User Manual.")
  788. # Check that the metaphlan2 executable can be found
  789. if not config.bypass_prescreen and not config.bypass_nucleotide_index:
  790. if not utilities.find_exe_in_path("metaphlan"):
  791. sys.exit("CRITICAL ERROR: The metaphlan executable can not be found. "
  792. "Please check the install.")
  794. # Check the metaphlan2 version
  795. # utilities.check_software_version("metaphlan",config.metaphlan_version)
  797. # Check that the bowtie2 executable can be found
  798. if not config.bypass_nucleotide_search:
  799. if not utilities.find_exe_in_path("bowtie2"):
  800. sys.exit("CRITICAL ERROR: The bowtie2 executable can not be found. "
  801. "Please check the install.")
  803. # Check the bowtie2 version
  804. utilities.check_software_version("bowtie2", config.bowtie2_version, warning=True)
  806. if not config.bypass_translated_search:
  807. # Check that the protein database directory exists
  808. if not os.path.isdir(config.protein_database):
  809. if args.protein_database:
  810. sys.exit("CRITICAL ERROR: The directory provided for the protein database at "
  811. + args.protein_database + " does not exist. Please select another directory.")
  812. else:
  813. sys.exit("CRITICAL ERROR: The default protein database directory of "
  814. + config.protein_database + " does not exist. Please provide the location "
  815. + "of the directory using the --protein-database option.")
  817. # Check that some files in the protein database folder are of the expected extension
  818. expected_database_extension=""
  819. if config.translated_alignment_selected == "usearch":
  820. expected_database_extension=config.usearch_database_extension
  821. elif config.translated_alignment_selected == "rapsearch":
  822. expected_database_extension=config.rapsearch_database_extension
  823. elif config.translated_alignment_selected == "diamond":
  824. expected_database_extension=config.diamond_database_extension
  826. valid_format_count=0
  827. database_files=os.listdir(config.protein_database)
  828. valid_format_database_files=[]
  829. for file in database_files:
  830. if not config.matching_uniref in file:
  831. sys.exit("\n\nCRITICAL ERROR: The directory provided for the translated database contains files ( "+file+" )"+\
  832. " that are not of the expected version. Please install the latest version"+\
  833. " of the database: "+config.matching_uniref)
  835. if file.endswith(expected_database_extension):
  836. # if rapsearch check for the second database file
  837. if config.translated_alignment_selected == "rapsearch":
  838. database_file=re.sub(config.rapsearch_database_extension+"$","",file)
  839. if database_file in database_files:
  840. valid_format_count+=1
  841. valid_format_database_files.append(database_file)
  842. else:
  843. valid_format_count+=1
  844. valid_format_database_files.append(file)
  846. if valid_format_count == 0:
  847. sys.exit("CRITICAL ERROR: The protein database directory provided ( " + config.protein_database
  848. + " ) does not contain any files that have been formatted to run with"
  849. " the translated alignment software selected ( " +
  850. config.translated_alignment_selected + " ). Please format these files so"
  851. + " they are of the expected extension ( " + expected_database_extension +" ).")
  853. # Check if running with the demo database
  854. if not config.bypass_translated_search:
  855. if os.path.basename(config.protein_database) == "uniref_DEMO":
  856. # Check the input file is a demo input if running with demo database
  857. try:
  858. input_file_size=os.path.getsize(args.input)/1024**2
  859. except EnvironmentError:
  860. input_file_size=0
  861. if input_file_size > MAX_SIZE_DEMO_INPUT_FILE:
  862. sys.exit("ERROR: You are using the demo UniRef database with "
  863. + "a non-demo input file. If you have not already done so, please "
  864. + "run humann_databases to download the full UniRef database. "
  865. + "If you have downloaded the full database, use the option "
  866. + "--protein-database to provide the location. "
  867. + "You can also run humann_config to update the default "
  868. + "database location. For additional information, please "
  869. + "see the HUMAnN User Manual.")
  871. # Check that the translated alignment executable can be found
  872. if not utilities.find_exe_in_path(config.translated_alignment_selected):
  873. sys.exit("CRITICAL ERROR: The " + config.translated_alignment_selected +
  874. " executable can not be found. Please check the install.")
  876. # Check for correct usearch version
  877. if config.translated_alignment_selected == "usearch":
  878. utilities.check_software_version("usearch",config.usearch_version)
  880. # Check for correct rapsearch version
  881. if config.translated_alignment_selected == "rapsearch":
  882. utilities.check_software_version("rapsearch", config.rapsearch_version)
  884. # Check for the correct diamond version
  885. if config.translated_alignment_selected == "diamond":
  886. utilities.check_software_version("diamond", config.diamond_version)
  888. # parse the chocophlan gene index
  889. chocophlan_gene_indexes=parse_chocophlan_gene_indexes(args.annotation_gene_index)
  890. config.chocophlan_gene_indexes = chocophlan_gene_indexes
  892. # set the identity thresholds if provided by the user
  893. config.nucleotide_identity_threshold=args.nucleotide_identity_threshold
  894. config.identity_threshold=args.identity_threshold
  896. def timestamp_message(task, start_time):
  897. """
  898. Print and log a message about the task completed and the time
  899. Log messages are tab delimited for quick task/time access with awk
  900. Return the new start time
  901. """
  902. message="TIMESTAMP: Completed \t" + task + " \t:\t " + \
  903. str(int(round(time.time() - start_time))) + "\t seconds"
  904. logger.info(message)
  905. if config.verbose:
  906. print("\n"+message.replace("\t","")+"\n")
  908. return time.time()
  910. def main():
  911. # Parse arguments from command line
  912. args=parse_arguments(sys.argv)
  914. # Update the configuration settings based on the arguments
  915. log_file=update_configuration(args)
  917. # Check for required files, software, databases, and also permissions
  918. check_requirements(args)
  920. # Write the config settings to the log file
  921. config.log_settings()
  923. # Initialize alignments and gene scores
  924. minimize_memory_use=True
  925. if config.memory_use == "maximum":
  926. minimize_memory_use=False
  928. alignments=store.Alignments(minimize_memory_use=minimize_memory_use)
  929. unaligned_reads_store=store.Reads(minimize_memory_use=minimize_memory_use)
  930. gene_scores=store.GeneScores()
  932. # If id mapping is provided then process
  933. if args.id_mapping:
  934. alignments.process_id_mapping(args.id_mapping)
  936. # Load in the reactions database
  937. reactions_database=None
  938. if config.pathways_database_part1:
  939. reactions_database=store.ReactionsDatabase(config.pathways_database_part1)
  941. message="Load pathways database part 1: " + config.pathways_database_part1
  942. logger.info(message)
  944. # Load in the pathways database
  945. pathways_database=store.PathwaysDatabase(config.pathways_database_part2, reactions_database)
  947. if config.pathways_database_part1:
  948. message="Load pathways database part 2: " + config.pathways_database_part2
  949. else:
  950. message="Load pathways database: " + config.pathways_database_part2
  951. logger.info(message)
  953. # Start timer
  954. start_time=time.time()
  956. # Process fasta or fastq input files
  957. output_files=[]
  958. if args.input_format in ["fasta","fastq"]:
  959. # Run prescreen to identify bugs
  960. bug_file = "Empty"
  961. if args.taxonomic_profile:
  962. bug_file = os.path.abspath(args.taxonomic_profile)
  963. else:
  964. if not config.bypass_prescreen:
  965. bug_file = prescreen.alignment(args.input)
  966. output_files.append(bug_file)
  967. start_time=timestamp_message("prescreen",start_time)
  969. # Create the custom database from the bugs list
  970. custom_database = ""
  971. if not config.bypass_nucleotide_index:
  972. custom_database = prescreen.create_custom_database(config.nucleotide_database, bug_file)
  973. start_time=timestamp_message("custom database creation",start_time)
  974. else:
  975. custom_database = "Bypass"
  977. # Run nucleotide search on custom database
  978. if custom_database != "Empty" and not config.bypass_nucleotide_search:
  979. if not config.bypass_nucleotide_index:
  980. nucleotide_index_file = nucleotide.index(custom_database)
  981. start_time=timestamp_message("database index",start_time)
  982. else:
  983. nucleotide_index_file = nucleotide.find_index(config.nucleotide_database)
  985. nucleotide_alignment_file = nucleotide.alignment(args.input,
  986. nucleotide_index_file)
  988. start_time=timestamp_message("nucleotide alignment",start_time)
  990. # Determine which reads are unaligned and reduce aligned reads file
  991. # Remove the alignment_file as we only need the reduced aligned reads file
  992. [ unaligned_reads_file_fasta, reduced_aligned_reads_file ] = nucleotide.unaligned_reads(
  993. nucleotide_alignment_file, alignments, unaligned_reads_store, keep_sam=True)
  995. start_time=timestamp_message("nucleotide alignment post-processing",start_time)
  997. # Print out total alignments per bug
  998. message="Total bugs from nucleotide alignment: " + str(alignments.count_bugs())
  999. logger.info(message)
  1000. print(message)
  1002. message=alignments.counts_by_bug()
  1003. logger.info("\n"+message)
  1004. print(message)
  1006. message="Total gene families from nucleotide alignment: " + str(alignments.count_genes())
  1007. logger.info(message)
  1008. print("\n"+message)
  1010. # Report reads unaligned
  1011. message="Unaligned reads after nucleotide alignment: " + utilities.estimate_unaligned_reads_stored(
  1012. args.input, unaligned_reads_store) + " %"
  1013. logger.info(message)
  1014. print("\n"+message+"\n")
  1015. else:
  1016. logger.debug("Custom database is empty")
  1017. reduced_aligned_reads_file = "Empty"
  1018. unaligned_reads_file_fasta=args.input
  1019. unaligned_reads_store=store.Reads(unaligned_reads_file_fasta, minimize_memory_use=minimize_memory_use)
  1021. # Do not run if set to bypass translated search in config file
  1022. if not config.bypass_translated_search:
  1023. # Run translated search on UniRef database if unaligned reads exit
  1024. if unaligned_reads_store.count_reads()>0:
  1025. translated_alignment_file = translated.alignment(config.protein_database,
  1026. unaligned_reads_file_fasta)
  1028. start_time=timestamp_message("translated alignment",start_time)
  1030. # Determine which reads are unaligned
  1031. translated_unaligned_reads_file_fastq = translated.unaligned_reads(
  1032. unaligned_reads_store, translated_alignment_file, alignments)
  1034. start_time=timestamp_message("translated alignment post-processing",start_time)
  1036. # Print out total alignments per bug
  1037. message="Total bugs after translated alignment: " + str(alignments.count_bugs())
  1038. logger.info(message)
  1039. print(message)
  1041. message=alignments.counts_by_bug()
  1042. logger.info("\n"+message)
  1043. print(message)
  1045. message="Total gene families after translated alignment: " + str(alignments.count_genes())
  1046. logger.info(message)
  1047. print("\n"+message)
  1049. # Report reads unaligned
  1050. message="Unaligned reads after translated alignment: " + utilities.estimate_unaligned_reads_stored(
  1051. args.input, unaligned_reads_store) + " %"
  1052. logger.info(message)
  1053. print("\n"+message+"\n")
  1054. else:
  1055. message="All reads are aligned so translated alignment will not be run"
  1056. logger.info(message)
  1057. print(message)
  1058. else:
  1059. message="Bypass translated search"
  1060. logger.info(message)
  1061. print(message)
  1063. # Process input files of sam format
  1064. elif args.input_format in ["sam"]:
  1066. # Store the sam mapping results
  1067. message="Process the sam mapping results ..."
  1068. logger.info(message)
  1069. print("\n"+message)
  1071. [unaligned_reads_file_fasta, reduced_aligned_reads_file] = nucleotide.unaligned_reads(
  1072. args.input, alignments, unaligned_reads_store, keep_sam=True)
  1074. start_time=timestamp_message("alignment post-processing",start_time)
  1076. # Process input files of tab-delimited blast format
  1077. elif args.input_format in ["blastm8"]:
  1079. # Store the blastm8 mapping results
  1080. message="Process the blastm8 mapping results ..."
  1081. logger.info(message)
  1082. print("\n"+message)
  1084. translated_unaligned_reads_file_fastq = translated.unaligned_reads(
  1085. unaligned_reads_store, args.input, alignments)
  1087. start_time=timestamp_message("alignment post-processing",start_time)
  1089. # Get the number of remaining unaligned reads
  1090. unaligned_reads_count=unaligned_reads_store.count_reads()
  1092. # Clear all of the unaligned reads as they are no longer needed
  1093. unaligned_reads_store.clear()
  1095. # Compute or load in gene families
  1096. if args.input_format in ["fasta","fastq","sam","blastm8"]:
  1097. # Compute the gene families
  1098. message="Computing gene families ..."
  1099. logger.info(message)
  1100. print("\n"+message)
  1102. families_file=families.gene_families(alignments,gene_scores,unaligned_reads_count)
  1103. output_files.append(families_file)
  1105. start_time=timestamp_message("computing gene families",start_time)
  1107. elif args.input_format in ["genetable"]:
  1108. # Load the gene scores
  1109. message="Process the gene table ..."
  1110. logger.info(message)
  1111. print("\n"+message)
  1113. unaligned_reads_count=gene_scores.add_from_file(args.input,id_mapping_file=args.id_mapping)
  1115. start_time=timestamp_message("processing gene table",start_time)
  1117. # Handle input files of unknown formats
  1118. else:
  1119. sys.exit("CRITICAL ERROR: Input file of unknown format.")
  1121. # Clear all of the alignments data as they are no longer needed
  1122. alignments.clear()
  1124. # Identify reactions and then pathways from the alignments
  1125. message="Computing reaction and pathway abundance ..."
  1126. logger.info(message)
  1127. print("\n"+message)
  1128. pathways_and_reactions_store=modules.identify_reactions_and_pathways(
  1129. gene_scores, reactions_database, pathways_database, unaligned_reads_count)
  1131. # Compute pathway abundance and coverage
  1132. abundance_file, coverage_file, reaction_file=modules.compute_pathways_abundance_and_coverage(
  1133. gene_scores, reactions_database, pathways_and_reactions_store, pathways_database, unaligned_reads_count)
  1134. output_files.append(reaction_file)
  1135. output_files.append(abundance_file)
  1136. output_files.append(log_file)
  1137. #output_files.append(coverage_file)
  1139. start_time=timestamp_message("computing pathways",start_time)
  1141. message="\nOutput files created: \n" + "\n".join(output_files) + "\n"
  1142. logger.info(message)
  1143. print(message)
  1145. # Remove the unnamed temp files
  1146. utilities.remove_directory(config.unnamed_temp_dir)
  1148. # Remove named temp directory
  1149. if args.remove_temp_output:
  1150. utilities.remove_directory(config.temp_dir)

humann.py at commit e07b3a3, under MIT · at the source

Overview

Authors: Amapola De Sales-Millan1, Paulina Reyes-Ferreira2, Rina María González-Cervantes3, Mariana Luna-Álvarez4, Sara Guillén-López5, José F. Cobo-Díaz6, Sandra Ramos7, José Félix Aguirre-Garrido3, José Antonio Velázquez-Aragón7
ORCID iDs: Amapola De Sales-Millan, Rina María González-Cervantes, Mariana Luna-Álvarez, Sara Guillén-López, José F. Cobo-Díaz, Sandra Ramos, José Félix Aguirre-Garrido, José Antonio Velázquez-Aragón
  1. Doctorado en Ciencias Biológicas y de la Salud, Universidad Autónoma Metropolitana, Ciudad de Mexico 04960, Mexico
  2. Departamento de Salud Mental, Instituto Nacional de Pediatría, Ciudad de Mexico 04530, Mexico
  3. Departamento de Ciencias Ambientales, Universidad Autónoma Metropolitana-Lerma, Lerma 52006, Estado de Mexico, Mexico
  4. Grupo Neurológico, Neuroquirúrgico y de Columna, Hospital Ángeles Acoxpa, Ciudad de Mexico 14308, Mexico
  5. Laboratorio de Errores Innatos del Metabolismo y Tamiz, Instituto Nacional de Pediatría, Ciudad de Mexico 04530, Mexico
  6. Departamento de Higiene y Tecnología de los Alimentos, Facultad de Veterinaria, Universidad de León, 24071 Leon, Spain
  7. Laboratorio de Oncología Experimental, Instituto Nacional de Pediatría, Ciudad de Mexico 04530, Mexico
Institutions: Universidad Autónoma Metropolitana (Mexico); Instituto Nacional de Pediatria (Mexico); Universidad de León (Spain)
Journal: Nutrients, volume 18, issue 15, article 2441
Dates: received 27 May 2026; accepted 22 July 2026; published online 26 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3390/nu18152441 · PMID 42588064 · PMCID PMC13468271 · OpenAlex W7171466161
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), autism (population), clinical / translational (subfield)
Methods: Statistics, Preprocessing, Connectivity
Keywords: autism spectrum disorder, probiotics, synbiotics, gut microbiota, microbiome function, ASD severity, gastrointestinal symptoms
MeSH: Autism Spectrum Disorder*, Gastrointestinal Microbiome*, Synbiotics*, Child, Child, Preschool, Dietary Supplements, Dysbiosis, Feces, Female, Humans, Longitudinal Studies, Male, Mexico, Probiotics, RNA, Ribosomal, 16S, Treatment Outcome (* major topic)
Topic: Gut microbiota and health (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: National Institute of Pediatrics
Citations: not cited yet (Europe PMC); 90 references in the paper

Abstract

Background/Objectives: Gut dysbiosis in children with autism spectrum disorder (ASD) has been associated with alterations in microbial ecology and metabolic function that may contribute to gastrointestinal dysfunction and the severity of clinical manifestations. Synbiotic and probiotic supplementation has emerged as a promising microbiome-targeted strategy for ASD; however, its effects on gut microbiome composition, functional potential, and clinical outcomes remain incompletely understood. We conducted a longitudinal study of Mexican children diagnosed with ASD to analyze changes in the composition, diversity, and functional potential of the gut microbiome during six months of multi-strain synbiotic supplementation. Methods: Stool samples were collected from 25 children with ASD at baseline and after 3 and 6 months of multi-strain synbiotic supplementation. Gut microbiome composition and diversity were analyzed by 16S rRNA gene sequencing, whereas whole metagenome sequencing (WMS) was performed in a subset of samples to evaluate the functional potential of the fecal microbiome. Gastrointestinal symptoms were assessed using the Rome IV criteria, and ASD severity was evaluated with the Childhood Autism Rating Scale (CARS). Results: Twenty-five children with ASD completed the 6 months of synbiotic supplementation. Overall, ASD severity decreased, reflected by a reduction in total CARS score, and improvements in several CARS domains. Gastrointestinal symptoms also decreased significantly. Longitudinal microbiome profiling revealed significant taxonomic and diversity changes over the supplementation period, while WMS identified changes in microbial metabolic potential, including enrichment of tryptophan biosynthesis pathways and reduced L-rhamnose degradation. Conclusions: This exploratory research provides proof-of-concept evidence supporting multi-strain synbiotic supplementation in children with ASD. Larger controlled studies are needed to confirm these findings and clarify their relevance to microbiota–gut–brain axis interactions. The observed concordance between clinical improvements and microbiome remodeling supports further investigation of microbiome-targeted interventions according to ASD severity and duration of supplementation.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

biobakery/humann

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: e07b3a34d0b94c09a8ac5d28ff95009611178be2, 10 July 2026
Languages: Python (61), Shell (10)
Size: 206 files, 71 scripts
Software Heritage: archived
Found in: the text, “2.6. Bioinformatic and Statistical Analysis”
Holds: README, license file, environment (setup.py), tests
Not found: CITATION.cff, continuous integration, documentation
Tools: NumPy (4 files), h5py (3 files), SciPy (2 files), Matplotlib (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
73 files
At the source: github.com/biobakery/humann#custom-protein-reference-database

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 71 scripts, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability Statement

All 16S rDNA and Whole Metagenomic Sequence files and corresponding mapping files for the samples utilized in this study have been deposited in the NCBI BioSample repository. Interested parties may access these resources using the provided Accession Number PRJNA1425957, which can be found at the following link: https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA1425957 (accessed on 27 May, 2026).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 7 keywords, 16 MeSH terms, 1 funder, 90 references.

Cite

This paper

De Sales-Millan, A., Reyes-Ferreira, P., González-Cervantes, R. M., Luna-Álvarez, M., Guillén-López, S., Cobo-Díaz, J. F., Ramos, S., Aguirre-Garrido, J. F., & Velázquez-Aragón, J. A. (2026). Clinical Improvement and Taxonomic-Functional Gut Microbiome Remodeling After Six Months of Multi-Strain Synbiotic Supplementation in Mexican Children with Autism Spectrum Disorder. Nutrients, 18(15), 2441. https://doi.org/10.3390/nu18152441

BibTeX

@article{desalesmillan2026clinical,
author = {De Sales-Millan, Amapola and Reyes-Ferreira, Paulina and González-Cervantes, Rina María and Luna-Álvarez, Mariana and Guillén-López, Sara and Cobo-Díaz, José F. and Ramos, Sandra and Aguirre-Garrido, José Félix and Velázquez-Aragón, José Antonio},
title = {{Clinical Improvement and Taxonomic-Functional Gut Microbiome Remodeling After Six Months of Multi-Strain Synbiotic Supplementation in Mexican Children with Autism Spectrum Disorder}},
journal = {Nutrients},
year = {2026},
month = jul,
volume = {18},
number = {15},
pages = {2441},
publisher = {Multidisciplinary Digital Publishing Institute (MDPI)},
issn = {2072-6643},
doi = {10.3390/nu18152441},
url = {https://doi.org/10.3390/nu18152441},
pmid = {42588064},
pmcid = {PMC13468271}
}

RIS

TY - JOUR
AU - De Sales-Millan, Amapola
AU - Reyes-Ferreira, Paulina
AU - González-Cervantes, Rina María
AU - Luna-Álvarez, Mariana
AU - Guillén-López, Sara
AU - Cobo-Díaz, José F.
AU - Ramos, Sandra
AU - Aguirre-Garrido, José Félix
AU - Velázquez-Aragón, José Antonio
TI - Clinical Improvement and Taxonomic-Functional Gut Microbiome Remodeling After Six Months of Multi-Strain Synbiotic Supplementation in Mexican Children with Autism Spectrum Disorder
T2 - Nutrients
J2 - Nutrients
PY - 2026
DA - 2026/07/26
VL - 18
IS - 15
SP - 2441
SN - 2072-6643
PB - Multidisciplinary Digital Publishing Institute (MDPI)
DO - 10.3390/nu18152441
UR - https://doi.org/10.3390/nu18152441
LA - en
ER -

CSL-JSON

{
"id": "10.3390/nu18152441",
"type": "article-journal",
"title": "Clinical Improvement and Taxonomic-Functional Gut Microbiome Remodeling After Six Months of Multi-Strain Synbiotic Supplementation in Mexican Children with Autism Spectrum Disorder",
"container-title": "Nutrients",
"author": [
{
"family": "De Sales-Millan",
"given": "Amapola"
},
{
"family": "Reyes-Ferreira",
"given": "Paulina"
},
{
"family": "González-Cervantes",
"given": "Rina María"
},
{
"family": "Luna-Álvarez",
"given": "Mariana"
},
{
"family": "Guillén-López",
"given": "Sara"
},
{
"family": "Cobo-Díaz",
"given": "José F."
},
{
"family": "Ramos",
"given": "Sandra"
},
{
"family": "Aguirre-Garrido",
"given": "José Félix"
},
{
"family": "Velázquez-Aragón",
"given": "José Antonio"
}
],
"container-title-short": "Nutrients",
"volume": "18",
"issue": "15",
"page": "2441",
"DOI": "10.3390/nu18152441",
"PMID": "42588064",
"PMCID": "PMC13468271",
"ISSN": "2072-6643",
"publisher": "Multidisciplinary Digital Publishing Institute (MDPI)",
"URL": "https://doi.org/10.3390/nu18152441",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
26
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1186/s12866-026-05296-x
Gut bacterial characteristics in children with autism spectrum disorder according to symptom severity: a cross-sectional study.
Journal: BMC microbiology
In common: autism, 4 references
[2] doi:10.1186/s12888-026-08178-8 [code]
Machine learning model to identify gut microbiome-derived metabolites as potential biomarkers of autism spectrum disorder: a pilot study.
Journal: BMC psychiatry
In common: Matplotlib, NumPy, autism, clinical / translational, 2 references
[3] doi:10.1016/j.bbih.2026.101275 [code]
A protocol for the Teen Bugs study: An integrative, multi-omics approach to understanding the role of the gut microbiome and mesocorticolimbic system in adolescent mental health following early adverse caregiving.
Journal: Brain, behavior, & immunity - health
In common: 4 references
[4] doi:10.1016/j.gutmic.2026.100006 [code]
From microbes to milestones: Gut bacterial abundances and functional pathways associate with neurodevelopment following preterm birth.
Journal: Gut microbiology
In common: autism, 3 references
[5] doi:10.1002/hbm.70496 [code]
Transdiagnostic Profiles of BOLD Signal Variability in Autism and Schizophrenia Spectrum Disorders: Associations With Cognition and Functioning.
Journal: Human brain mapping
In common: SciPy, Matplotlib, NumPy, autism, 1 reference
[6] doi:10.3389/fnins.2026.1858005 [code]
Targeted stool metabolomics suggests exploratory catecholamine- and tryptophan-linked metabolic features in autism spectrum disorder.
Journal: Frontiers in neuroscience
In common: autism, clinical / translational, 2 references
[7] doi:10.1016/j.celrep.2026.117590 [code]
Impaired behavioral inhibition in Fmr1 KO mice is linked to disrupted visual cortex theta oscillations.
Journal: Cell reports
In common: SciPy, Matplotlib, NumPy, autism, 1 reference
[8] doi:10.1038/s41467-026-76675-1 [code]
Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: h5py, SciPy, Matplotlib, 1 other tool, autism
[9] doi:10.1038/s41467-026-74320-5 [code]
Spatial architecture of autism pathogenesis reveals mosaic structural disarray during early development.
Journal: Nature communications
In common: h5py, SciPy, Matplotlib, 1 other tool, autism
[10] doi:10.1002/advs.202519479 [code]
Diminished Signal-to-Noise Ratio Disrupts Somatosensory Population Encoding and Drives Tactile Hyposensitivity in the Fmr1&lt;sup&gt;-/y&lt;/sup&gt; Autism Model.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: h5py, SciPy, Matplotlib, 1 other tool, autism

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.