OSCR

Non-concussive head impacts sustained during American football correlate with changes in gut microbiome diversity and composition.

Code ↔ Paper

12 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 12 matches
  1. [1] § Materials and methods › Fecal sample collection lifestyle questionnaires ↔ TBI_CLR.Rmd, lines 30–80 · score 0.82 · miceRanger, Bristol Stool, orthopedic injury, illness, nicotine, vomiting
  2. [2] § Materials and methods › Fecal sample collection lifestyle questionnaires ↔ TBI_Microbiome_R.Rmd, lines 31–83 · score 0.82 · miceRanger, Bristol Stool, orthopedic injury, illness, nicotine, vomiting
  3. [3] § Results › Bray-Curtis dissimilarity increases three days following substantial head impact exposure ↔ TBI Microbiome.ipynb, lines 1033–1165 · score 0.73 · 72–96 hours, 0–24 hours, Bray Curtis Dissimilarity, 48–72 hours, head impact, Nemenyi
  4. [4] § Results › Bray-Curtis dissimilarity increases three days following substantial head impact exposure ↔ TBI_Microbiome.ipynb, lines 978–1109 · score 0.73 · 72–96 hours, 0–24 hours, Bray Curtis Dissimilarity, 48–72 hours, head impact, Nemenyi
  5. [5] § Results › Bray-Curtis dissimilarity increases three days following substantial head impact exposure ↔ TBI Microbiome.ipynb, lines 1033–1165 · score 0.68 · 72–96 hours, 0–24 hour, Bray Curtis Dissimilarity, 48–72, head impact, Nemenyi
  6. [6] § Results › Bray-Curtis dissimilarity increases three days following substantial head impact exposure ↔ TBI_Microbiome.ipynb, lines 978–1109 · score 0.68 · 72–96 hours, 0–24 hour, Bray Curtis Dissimilarity, 48–72, head impact, Nemenyi
  7. [7] § Materials and methods › Data analysis – Data slicing and repeated-measures analysis ↔ TBI Microbiome.ipynb, lines 920–1029 · score 0.57 · post hoc pairwise, Rank, ANOVA, Nemenyi, Friedman, microbial
  8. [8] § Materials and methods › Data analysis – Data slicing and repeated-measures analysis ↔ TBI_Microbiome.ipynb, lines 865–974 · score 0.57 · post hoc pairwise, Rank, ANOVA, Nemenyi, Friedman, microbial
  9. [9] § Materials and methods › Data processing – Taxonomic data ↔ TBI Microbiome.ipynb, lines 273–291 · score 0.55 · Silva V4, classifier provided, taxonomic, sequence, microbial
  10. [10] § Materials and methods › Data processing – Taxonomic data ↔ TBI_Microbiome.ipynb, lines 280–298 · score 0.55 · Silva V4, classifier provided, taxonomic, sequence, microbial
  11. [11] § Results › Gut microbiome composition changes 48–72 hours post head impact ↔ TBI_CLR.Rmd, lines 251–321 · score 0.55 · CLR, circumstances, Orthopedic, caffeine, Coriobacteriales, sleep
  12. [12] § Results › Gut microbiome composition changes 48–72 hours post head impact ↔ TBI_Microbiome_R.Rmd, lines 312–385 · score 0.51 · circumstances, Orthopedic, caffeine, Coriobacteriales, sleep, Verrucomicrobiales

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

Jupyter notebook · 1,633 lines · 70 KB · no license · 4 matches

  1. # %% [markdown]
  2. # # General Notes
  4. # %% [markdown]
  5. # This Jupyter Notebook mostly uses command line commands instead of actual Python code to utilize QIIME 2. In the first line of code cells, you can often see `%%bash`, which is used to denote that the chunk contains commands rather than Python code.
  6. #
  7. # QIIME 2 has its own file formats: `.qza` and `.qzv`. You can consider `.qza` as files that contain data and `.qzv` files as visualizations of results. To view results stored in `.qzv` files, open [QIIME View](https://view.qiime2.org) (https://view.qiime2.org) and drag your `.qzv` file into the browser window. The file will load and display results.
  8. #
  9. # Make sure you are running this notebook with your QIIME 2 virtual environment. If you have followed QIIME 2's installation guide and installed successfully, you should have a dedicated Anaconda/Miniconda virtual environment for QIIME 2. If you are using this notebook with VSCode, select the environment that corresponds to QIIME 2 on the upper right side.
  10. #
  11. # If you have any questions about QIIME 2, they have detailed [documentation](https://docs.qiime2.org) online. Specifically for this notebook (time series analysis), many of the commands use the `longitudinal` module. You can check the tutorial on [longitudinal analysis](https://docs.qiime2.org/2022.11/tutorials/longitudinal/) specifically if you want more details on the commands.
  13. # %% [markdown]
  14. # # Setup
  16. # %% [markdown]
  17. # ### Import Libraries
  19. # %%
  20. %matplotlib inline
  21. import pandas as pd
  22. import os
  23. import numpy as np
  24. #import xlrd
  25. import matplotlib.pyplot as plt
  26. from sklearn import metrics
  27. from sklearn.linear_model import LinearRegression
  28. #import statsmodels.api as sm
  29. from sklearn.metrics import r2_score
  30. import seaborn as sns
  31. #import statsmodels.formula.api as smf
  32. #import patsy
  33. from typing import Union
  34. from scipy import stats
  35. from itertools import combinations
  36. from scipy.spatial.distance import braycurtis, pdist, squareform
  37. from scipy.stats import chi2, rankdata
  38. #from statsmodels.stats.anova import AnovaRM
  39. #from statsmodels.stats.libqsturng import psturng
  40. #from skbio import DistanceMatrix
  41. #from skbio.stats.ordination import pcoa
  42. #import scikit_posthocs
  43. import math
  44. import datetime
  45. import warnings
  46. warnings.filterwarnings('ignore')
  49. # %% [markdown]
  50. # ### Change working directory to your folder
  52. # %%
  53. #workdir='/Users/Aziz/Documents/ApMa Thesis/2023_03_14_012623KBillcus515F_Raw_Data_UDI/data_files' # Set this to your working directory
  54. workdir = '/Users/aziz/Documents/Ay Lab/2023_03_14_012623KBillcus515F_Raw_Data_UDI/data_files'
  55. %cd $workdir
  57. # %% [markdown]
  58. # # Load Files
  60. # %% [markdown]
  61. # ### Zip fastq Files
  63. # %% [markdown]
  64. # If you have files in .fastq format, run this chunk to zip it to .fastq.gz format so that importing is more efficient.
  65. #
  66. # If you already have files in .fastq.gz format, skip this chunk.
  68. # %%
  69. %%bash
  70. gzip FastQ/*.fastq
  72. # %% [markdown]
  73. # ### Construct Manifest File and convert Metadata
  75. # %% [markdown]
  76. # You need to manually create a manifest file, which specifies the file paths of .fastq.gz files for each sample.
  77. #
  78. # If you have paired end data, create this file in excel with three columns: "sampleid", "forward-absolute-filepath", "reverse-absolute-filepath". If you have single end data, you only need the first two columns.
  79. #
  80. # List the sample names in the "sampleid" column, and put the absolute file paths to the sample's forward and reverse reads in the next two columns. Then save the excel file as manifest.csv (make sure it is .csv), and run this chunk.
  81. #
  82. # For a reference of what a manifest file looks like, check QIIME 2 [documentation](https://docs.qiime2.org/2022.2/tutorials/importing/#fastq-manifest-formats).
  83. #
  84. # Below, there is some code to generate a manifest file using file names.
  86. # %%
  87. list_of_files = os.listdir(workdir)
  88. print(list_of_files)
  89. if '.DS_Store' in list_of_files:
  90. list_of_files.remove('.DS_Store')
  91. if "manifest.tsv" in list_of_files:
  92. os.remove("manifest.tsv")
  93. #print(list_of_files)
  94. separator = "_S1_"
  95. sample_ids = [filename.split(separator, 1)[0] for filename in list_of_files]
  96. sample_ids = [*set(sample_ids)]
  97. #print(sample_ids)
  98. forward_absolute_files = []
  99. reverse_absolute_files = []
  100. for id in sample_ids:
  101. #find the forward and reverse for this id
  102. forward_reverse = [filename for filename in list_of_files if id in filename]
  103. #print(forward_reverse)
  104. if "R1" in forward_reverse[0]:
  105. forward_absolute_files.append(os.path.join(workdir, forward_reverse[0]))
  106. reverse_absolute_files.append(os.path.join(workdir, forward_reverse[1]))
  107. else:
  108. forward_absolute_files.append(os.path.join(workdir, forward_reverse[1]))
  109. reverse_absolute_files.append(os.path.join(workdir, forward_reverse[0]))
  111. manifest_dict = {"sample-id": sample_ids, "forward-absolute-filepath": forward_absolute_files,
  112. "reverse-absolute-filepath": reverse_absolute_files}
  113. manifest = pd.DataFrame(data=manifest_dict)
  114. manifest.to_csv("manifest.tsv", sep = "\t", index = False)
  117. # %%
  118. manifest = pd.read_csv("manifest.csv")
  119. # manifest = manifest.iloc[:,:3] # This line of code keeps the first three columns, making sure there's no empty columns messing up the manifest. May not be necessary
  120. manifest.to_csv("manifest.tsv", sep="\t", index=False)
  122. # %% [markdown]
  123. # ### Import Sequences and View Quality Information
  125. # %% [markdown]
  126. # This chunk imports the sequences into QIIME.
  127. #
  128. # Make sure the `type` and `input-format` parameters are appropriate for your data. If you have single end data, change `type` to `SampleData[SequencesWithQuality]` and `input-format` to `SingleEndFastqManifestPhred33V2`.
  129. #
  130. # Make sure the offset for sequence quality score is 33. You should have this information from the people who conducted the sequencing when you get the sequences back. If the offset is 64 instead of 33, change the `33` to `64` in `input-format`.
  132. # %%
  133. %%bash
  134. qiime tools import --type 'SampleData[PairedEndSequencesWithQuality]' --input-path manifest.tsv --output-path demux.qza --input-format PairedEndFastqManifestPhred33V2
  136. # %% [markdown]
  137. # This next chunk visualizes the quality information. After you run this chunk, you should have a `demux.qzv` file in your folder. Open [QIIME View](https://view.qiime2.org) in your browser and drag `demux.qzv` in. Click on the `Interactive Quality Plot` tab at the upper left to view quality information.
  139. # %%
  140. %%bash
  141. qiime demux summarize \
  142. --i-data demux.qza \
  143. --o-visualization demux.qzv
  145. # %% [markdown]
  146. # # Sequence Trim and Denoise
  148. # %% [markdown]
  149. # There are two choices of denoising methods within QIIME: DADA2 and Deblur. There is no consensus on which is better, but you should have preferences depending on your data type.
  150. #
  151. # DADA2 works natively with paired end data, while Deblur does not. To use with Deblur, paired end sequences must be joined first with `qiime vsearch join-pairs` function (or some other joining action) to become joined sequences. However, parameters for the `vsearch join-pairs` operation are very tricky. They are not easily chosen and may lead to loss of read data.
  152. #
  153. # Typically, these two methods will give similar results for weighted analysis, but may have variations in unweighted analysis and alpha diversity metrics. DADA2 usually identifies noticeably more ASVs than Deblur in their output.
  154. #
  155. # DADA2 takes significantly longer and more memory to run than Deblur, so if you have huge single end read files, you may want to try Deblur first.
  156. #
  157. # Be aware that the outputted file names are different for these two methods in order to distinguish (e.g. dada2-table.qza and deblur-table.qza), so change file names accordingly for later analysis.
  159. # %% [markdown]
  160. # ### DADA2:
  162. # %% [markdown]
  163. # After viewing quality information in `demux.qzv`, you can determine the truncating spot on the left for forward (`p-trim-left-f`) and reverse (`p-trim-left-r`) reads. `p-trunc-len-f` and `p-trunc-len-r`are for truncating spots on the right for forward and reverse reads. `p-n-threads` determines how many cores of you CPU will be used. If you would like to use all cores, put `0` for this parameter.
  165. # %%
  166. %%bash
  167. qiime dada2 denoise-paired \
  168. --i-demultiplexed-seqs demux.qza \
  169. --p-trim-left-f 0 \
  170. --p-trim-left-r 0 \
  171. --p-trunc-len-f 220 \
  172. --p-trunc-len-r 180 \
  173. --p-n-threads 8 \
  174. --o-table dada2-table.qza \
  175. --o-representative-sequences dada2-rep-seqs.qza \
  176. --o-denoising-stats dada2-denoising-stats.qza
  178. # %%
  179. %%bash
  180. qiime feature-table summarize \
  181. --i-table dada2-table.qza \
  182. --o-visualization dada2-table.qzv \
  183. --m-sample-metadata-file metadata.tsv
  185. qiime feature-table tabulate-seqs \
  186. --i-data dada2-rep-seqs.qza \
  187. --o-visualization dada2-rep-seqs.qzv
  189. qiime metadata tabulate \
  190. --m-input-file dada2-denoising-stats.qza \
  191. --o-visualization dada2-denoising-stats.qzv
  193. # %% [markdown]
  194. # ### Deblur:
  196. # %% [markdown]
  197. # If you have paired end reads, join them first by running this chunk.
  198. #
  199. # For `qiime vsearch join-pairs` parameters, refer to [documentation](https://docs.qiime2.org/2022.2/plugins/available/vsearch/join-pairs/).
  201. # %%
  202. %%bash
  203. qiime vsearch join-pairs \
  204. --i-demultiplexed-seqs demux.qza \
  205. --p-truncqual 10 \
  206. --p-threads 0 \
  207. --o-joined-sequences deblur-demux.qza
  209. qiime demux summarize \
  210. --i-data deblur-demux.qza \
  211. --o-visualization deblur-demux.qzv
  213. # %% [markdown]
  214. # For `qiime deblur denoise-16S` parameters, check quality information in `deblur-demux.qzv`. `p-left-trim-len` is the truncating spot on the left, and `p-trim-length` is the truncating spot on the right. `p-jobs-to-start` is the number of CPU cores being used by this command.
  216. # %%
  217. %%bash
  218. qiime deblur denoise-16S \
  219. --i-demultiplexed-seqs joined-demux.qza \
  220. --p-trim-length 240 \
  221. --p-left-trim-len 13 \
  222. --p-jobs-to-start 4 \
  223. --o-table deblur-table.qza \
  224. --o-representative-sequences deblur-rep-seqs.qza \
  225. --o-stats deblur-denoising-stats.qza
  227. # %% [markdown]
  228. # # Process the Feature Table
  230. # %% [markdown]
  231. # Now we have obtained feature table and representative sequences, we can further process and clean them.
  232. #
  233. # One option is to transform the feature table into relative abundance table. Number of sequence reads may vary greatly across samples/individuals. Transforming the table into relative abundance effectively controls the library size (sum of reads for a sample). Additionally, in some cases microbiome data should be considered compositional and therefore should be represented in percentages.
  235. # %%
  236. %%bash
  237. qiime feature-table relative-frequency \
  238. --i-table dada2-filtered-table.qza \
  239. --o-relative-frequency-table dada2-relative-table.qza \
  241. # %% [markdown]
  242. # You can also filter by feature prevalence. In the code below, we filter out features whose sum of frequency across all samples is less than 10, and also features present in fewer than 3 samples. You should choose these parameters after inspecting your `table.qzv` file. If you chose to convert your table to relative abundance already, make sure to pass in the relative table in `i-table`, and change `min-frequency` to something like `0.0001`.
  244. # %%
  245. %%bash
  246. qiime feature-table filter-features \
  247. --i-table dada2-table.qza \
  248. --p-min-frequency 10 \
  249. --p-min-samples 3 \
  250. --o-filtered-table dada2-filtered-table.qza
  252. # %% [markdown]
  253. # # Collapsing by Taxa
  255. # %% [markdown]
  258. # %% [markdown]
  259. # If you need to look at the feature table as a tsv (viewable in Excel), you can run the below two commands. The first one converts the `.qza` file to a biom table. Then we can convert the biom table to a tsv.
  261. # %%
  262. %%bash
  263. qiime tools export \
  264. --input-path dada2-table.qza \
  265. --output-path dada2-table
  267. # %%
  268. %%bash
  269. biom convert \
  270. -i dada2-table/feature-table.biom \
  271. -o dada2-table.tsv --to-tsv
  273. # %% [markdown]
  274. # # Assign Taxonomy
  276. # %% [markdown]
  277. # Next, we perform taxonomic assignment. QIIME provides classifiers trained on two gene databases: SILVA and Greengenes.
  278. #
  279. # SILVA database has a bigger gene tree, is regularly updated, and probably will take more time to run. Greengenes is more compact, was last updated in 2013 (so maybe relatively out-dated), but sometimes it is considered to be more accurate when dealing with specifically human gut microbiome.
  280. #
  281. # Below we download two pre-trained Naive Bayes classifiers provided by QIIME. One is trained on SILVA and the other on Greengenes database. Both of these are trained on sequences from only the V4 region in 16S rRNA sequencing. If your sequences are full-length 16S rRNA reads, then get rid of the `-515-806` part of the download link to download classifiers on the full region.
  283. # %%
  284. %%bash
  285. wget \
  286. -O 'Silva-V4-classifier.qza' \
  287. 'https://data.qiime2.org/2020.6/common/silva-138-99-515-806-nb-classifier.qza'
  289. #wget \
  290. # -O 'Greengenes-V4-classifier.qza' \
  291. # 'https://data.qiime2.org/2022.2/common/gg-13-8-99-515-806-nb-classifier.qza'
  293. # %% [markdown]
  294. # Now we use the classifier to assign taxonomy. The code here uses SILVA database. If you want to use Greengenes, change `i-classifier` to `Greengenes-V4-classifier.qza`.
  296. # %%
  297. %%bash
  298. qiime feature-classifier classify-sklearn \
  299. --i-classifier Silva-V4-classifier.qza \
  300. --i-reads dada2-rep-seqs.qza \
  301. --o-classification taxonomy.qza
  303. # %% [markdown]
  304. # # Collapsing by Taxa
  305. # You can also collapse the feature table to different taxonomic levels, i.e. the columns of your feature table will be genera/families rather than species. `p-level` denotes which taxonomic level you want to collapse to. `7` is the lowest, the species level.
  307. # %% [markdown]
  308. # Remember: King Philip came over for good soup
  310. # %%
  311. %%bash
  312. qiime taxa collapse \
  313. --i-table dada2-table.qza \
  314. --i-taxonomy taxonomy.qza \
  315. --p-level 2 \
  316. --o-collapsed-table phyla-table.qza
  318. # %% [markdown]
  319. # After the taxonomic classifier is ready, you can prepare data table for different taxa levels. It is necessary to store both the "filtered" version which has raw counts and the "relative" version which has those counts normalized.
  320. #
  321. # 1. Genus Level.
  323. # %%
  324. %%bash
  325. qiime feature-table filter-features \
  326. --i-table genus-table.qza \
  327. --p-min-frequency 10 \
  328. --p-min-samples 3 \
  329. --o-filtered-table RawCount_Tables/genus-filtered-table.qza
  331. # %%
  332. %%bash
  333. qiime feature-table relative-frequency \
  334. --i-table RawCount_Tables/genus-filtered-table.qza \
  335. --o-relative-frequency-table Relative_Tables/genus-relative-table.qza \
  337. # %%
  338. %%bash
  339. qiime tools export \
  340. --input-path Relative_Tables/genus-relative-table.qza \
  341. --output-path Relative_Tables/genus-table
  343. # %%
  344. %%bash
  345. biom convert \
  346. -i Relative_Tables/genus-table/feature-table.biom \
  347. -o Relative_Tables/genus-table.tsv --to-tsv
  349. # %%
  350. %%bash
  351. qiime tools export \
  352. --input-path RawCount_Tables/genus-filtered-table.qza \
  353. --output-path RawCount_Tables/genus-table
  355. # %%
  356. %%bash
  357. biom convert \
  358. -i RawCount_Tables/genus-table/feature-table.biom \
  359. -o RawCount_Tables/genus-table.tsv --to-tsv
  361. # %% [markdown]
  362. # 2. Family Level
  364. # %%
  365. %%bash
  366. qiime feature-table filter-features \
  367. --i-table family-table.qza \
  368. --p-min-frequency 10 \
  369. --p-min-samples 3 \
  370. --o-filtered-table RawCount_Tables/family-filtered-table.qza
  372. # %%
  373. %%bash
  374. qiime feature-table relative-frequency \
  375. --i-table RawCount_Tables/family-filtered-table.qza \
  376. --o-relative-frequency-table Relative_Tables/family-relative-table.qza \
  378. # %%
  379. %%bash
  380. qiime tools export \
  381. --input-path Relative_Tables/family-relative-table.qza \
  382. --output-path Relative_Tables/family-table
  384. # %%
  385. %%bash
  386. biom convert \
  387. -i Relative_Tables/family-table/feature-table.biom \
  388. -o Relative_Tables/family-table.tsv --to-tsv
  390. # %%
  391. %%bash
  392. qiime tools export \
  393. --input-path RawCount_Tables/family-filtered-table.qza \
  394. --output-path RawCount_Tables/family-table
  396. # %%
  397. %%bash
  398. biom convert \
  399. -i RawCount_Tables/family-table/feature-table.biom \
  400. -o RawCount_Tables/family-table.tsv --to-tsv
  402. # %% [markdown]
  403. # 3. Order Level
  405. # %%
  406. %%bash
  407. qiime feature-table filter-features \
  408. --i-table order-table.qza \
  409. --p-min-frequency 10 \
  410. --p-min-samples 3 \
  411. --o-filtered-table RawCount_Tables/order-filtered-table.qza
  413. # %%
  414. %%bash
  415. qiime feature-table relative-frequency \
  416. --i-table RawCount_Tables/order-filtered-table.qza \
  417. --o-relative-frequency-table Relative_Tables/order-relative-table.qza \
  419. # %%
  420. %%bash
  421. qiime tools export \
  422. --input-path Relative_Tables/order-relative-table.qza \
  423. --output-path Relative_Tables/order-table
  425. # %%
  426. %%bash
  427. biom convert \
  428. -i Relative_Tables/order-table/feature-table.biom \
  429. -o Relative_Tables/order-table.tsv --to-tsv
  431. # %%
  432. %%bash
  433. qiime tools export \
  434. --input-path RawCount_Tables/order-filtered-table.qza \
  435. --output-path RawCount_Tables/order-table
  437. # %%
  438. %%bash
  439. biom convert \
  440. -i RawCount_Tables/order-table/feature-table.biom \
  441. -o RawCount_Tables/order-table.tsv --to-tsv
  443. # %% [markdown]
  444. # 4. Species Level
  446. # %%
  447. %%bash
  448. qiime feature-table filter-features \
  449. --i-table species-table.qza \
  450. --p-min-frequency 10 \
  451. --p-min-samples 3 \
  452. --o-filtered-table RawCount_Tables/species-filtered-table.qza
  454. # %%
  455. %%bash
  456. qiime feature-table relative-frequency \
  457. --i-table RawCount_Tables/species-filtered-table.qza \
  458. --o-relative-frequency-table Relative_Tables/species-relative-table.qza \
  460. # %%
  461. %%bash
  462. qiime tools export \
  463. --input-path Relative_Tables/species-relative-table.qza \
  464. --output-path Relative_Tables/species-table
  466. # %%
  467. %%bash
  468. biom convert \
  469. -i Relative_Tables/species-table/feature-table.biom \
  470. -o Relative_Tables/species-table.tsv --to-tsv
  472. # %%
  473. %%bash
  474. qiime tools export \
  475. --input-path RawCount_Tables/species-filtered-table.qza \
  476. --output-path RawCount_Tables/species-table
  478. # %%
  479. %%bash
  480. biom convert \
  481. -i RawCount_Tables/species-table/feature-table.biom \
  482. -o RawCount_Tables/species-table.tsv --to-tsv
  484. # %% [markdown]
  485. # 5. Phylum Level
  487. # %%
  488. %%bash
  489. qiime feature-table filter-features \
  490. --i-table phyla-table.qza \
  491. --p-min-frequency 10 \
  492. --p-min-samples 3 \
  493. --o-filtered-table RawCount_Tables/phyla-filtered-table.qza
  495. # %%
  496. %%bash
  497. qiime feature-table relative-frequency \
  498. --i-table RawCount_Tables/phyla-filtered-table.qza \
  499. --o-relative-frequency-table Relative_Tables/phyla-relative-table.qza \
  501. # %%
  502. %%bash
  503. qiime tools export \
  504. --input-path Relative_Tables/phyla-relative-table.qza \
  505. --output-path Relative_Tables/phyla-table
  507. # %%
  508. %%bash
  509. biom convert \
  510. -i Relative_Tables/phyla-table/feature-table.biom \
  511. -o Relative_Tables/phyla-table.tsv --to-tsv
  513. # %%
  514. %%bash
  515. qiime tools export \
  516. --input-path RawCount_Tables/phyla-filtered-table.qza \
  517. --output-path RawCount_Tables/phyla-table
  519. # %%
  520. %%bash
  521. biom convert \
  522. -i RawCount_Tables/phyla-table/feature-table.biom \
  523. -o RawCount_Tables/phyla-table.tsv --to-tsv
  525. # %% [markdown]
  526. # # Create Phylogenetic Tree
  528. # %% [markdown]
  529. # We construct a phylogenetic tree based on the representative sequences information. Three algorithms for constructing trees are available in QIIME 2: FastTree, IQ-Tree, and RAxML. FastTree runs the fastest but returns poorest quality. IQ-Tree and RAxML will give somewhat different results, but how exactly they differ varies between datasets. Here the code uses RAxML. If you wish to use any of the other method, simply change `raxml` in the command to `fasttree` or `iqtree`.
  531. # %%
  532. %%bash
  533. qiime phylogeny align-to-tree-mafft-raxml \
  534. --i-sequences dada2-rep-seqs.qza \
  535. --p-n-threads 4 \
  536. --output-dir phylogeny-trees
  538. # %%
  539. %%bash
  540. qiime feature-table filter-seqs \
  541. --i-data dada2-rep-seqs.qza \
  542. --i-table dada2-filtered-table.qza \
  543. --o-filtered-data filtered-rep-seqs.qza
  545. # %% [markdown]
  546. # # Calculating Diversity Metrics
  547. #
  548. # 1.1 Alpha Diversity (Faith's Phylogenetic Diversity)
  550. # %%
  551. %%bash
  552. qiime diversity alpha-phylogenetic \
  553. --i-table dada2-table.qza \
  554. --i-phylogeny phylogeny-trees/rooted_tree.qza \
  555. --p-metric faith_pd \
  556. --o-alpha-diversity species-faith-table.qza
  558. # %%
  559. %%bash
  560. qiime tools export \
  561. --input-path species-faith-table.qza \
  562. --output-path species-faith-table
  564. # %%
  565. %%bash
  566. biom convert \
  567. -i species-faith-table/feature-table.biom \
  568. -o --to-tsv
  570. # %% [markdown]
  571. # 1.2 Simpson's Alpha Diversity
  573. # %%
  574. %%bash
  575. qiime diversity alpha \
  576. --i-table species-filtered-table.qza \
  577. --p-metric simpson \
  578. --o-alpha-diversity species-simpson-table.qza
  580. # %%
  581. %%bash
  582. qiime tools export \
  583. --input-path species-simpson-table.qza \
  584. --output-path species-simpson-table
  586. # %%
  587. %%bash
  588. biom convert \
  589. -i species-simpson-table/feature-table.biom \
  590. -o --to-tsv
  592. # %% [markdown]
  593. # 2.1 Beta Diversity (Unweighted UniFrac -- Phylogenetic)
  595. # %%
  596. %%bash
  597. qiime diversity beta-phylogenetic \
  598. --i-table dada2-table.qza \
  599. --i-phylogeny phylogeny-trees/rooted_tree.qza \
  600. --p-metric unweighted_unifrac \
  601. --o-distance-matrix unweighted_unifrac_distance_matrix.qza
  603. # %%
  604. %%bash
  605. qiime tools export \
  606. --input-path unweighted_unifrac_distance_matrix.qza \
  607. --output-path unweighted_unifrac_distance_matrix
  609. # %% [markdown]
  610. # 2.11 Processing Beta Diversity Distance Matrix
  612. # %%
  613. unifrac_distance_matrix = pd.read_csv("unweighted_unifrac_distance_matrix/distance-matrix.tsv",sep="\t", index_col=0, header=0)
  614. unifrac_distance_matrix.drop(axis = 0, index=unifrac_distance_matrix.columns[-5:], inplace = True)
  615. unifrac_distance_matrix.drop(axis = 1, columns=unifrac_distance_matrix.columns[-5:], inplace = True)
  617. uni_1 = unifrac_distance_matrix.loc[unifrac_distance_matrix.columns.str.startswith("01."),unifrac_distance_matrix.columns.str.startswith("01.")].iloc[:,0]
  618. uni_4 = unifrac_distance_matrix.loc[unifrac_distance_matrix.columns.str.startswith("04."),unifrac_distance_matrix.columns.str.startswith("04.")].iloc[:,0]
  619. uni_5 = unifrac_distance_matrix.loc[unifrac_distance_matrix.columns.str.startswith("05."),unifrac_distance_matrix.columns.str.startswith("05.")].iloc[:,0]
  620. uni_8 = unifrac_distance_matrix.loc[unifrac_distance_matrix.columns.str.startswith("08."),unifrac_distance_matrix.columns.str.startswith("08.")].iloc[:,0]
  621. uni_9 = unifrac_distance_matrix.loc[unifrac_distance_matrix.columns.str.startswith("09."),unifrac_distance_matrix.columns.str.startswith("09.")].iloc[:,0]
  622. uni_16 = unifrac_distance_matrix.loc[unifrac_distance_matrix.columns.str.startswith("16."),unifrac_distance_matrix.columns.str.startswith("16.")].iloc[:,0]
  624. unifrac_distance = pd.concat([uni_1,uni_4, uni_5, uni_8, uni_9, uni_16], axis = 0).to_frame(name = "UniFrac_Distance")
  625. unifrac_distance.at["09.0808", "UniFrac_Distance"] = (unifrac_distance.at["09.0808.1", "UniFrac_Distance"] + unifrac_distance.at["09.0808.2", "UniFrac_Distance"]) /2
  627. print(unifrac_distance)
  631. # %% [markdown]
  632. # 2.2 Bray Curtis Dissimilarity
  634. # %%
  635. species_table = pd.read_csv("Relative_Tables/species-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  636. species_table = species_table.T
  637. species_table.index.names = ["Sample_ID"]
  638. species_table.drop(species_table.tail(5).index,inplace=True) # drop last n rows
  640. species = species_table.columns
  641. def apply_BC(df):
  642. first_row = df.loc[df.index[0], :].values.flatten().tolist()[0:-2]
  643. for index in df.index.values:
  644. row = df.loc[index, :].values.flatten().tolist()[0:-2]
  645. df.at[index, "BC_Dissimilarity"] = braycurtis(row, first_row)
  646. return df
  647. species_table["Player"] = [i.split(".")[0] for i in species_table.index.tolist()]
  648. species_table["BC_Dissimilarity"] = np.nan
  649. species_table = species_table.groupby("Player").apply(lambda df: apply_BC(df))
  650. species_table.index = species_table.index.droplevel(0)
  651. print(species_table)
  653. species_table.at["09.0808", "BC_Dissimilarity"] = (species_table.at["09.0808.1", "BC_Dissimilarity"] + species_table.at["09.0808.2", "BC_Dissimilarity"]) /2
  654. print(species_table.loc[:,"BC_Dissimilarity"])
  656. # %% [markdown]
  657. # # Combine Diversity Data, Fecal Sample data, and Survey Data
  659. # %%
  660. master_data = pd.read_csv("../Master Spreadsheet - Master Data Sheet.tsv",sep = "\t",dtype={"ID/Study Date ":"str"})
  661. master_data = master_data.set_index(master_data.columns[0])
  663. #Load Faith's Diversity
  664. alpha_diversity_table = pd.read_csv("species-faith-table/alpha-diversity.tsv", sep = "\t",
  665. index_col=0)
  666. alpha_diversity_table.drop(alpha_diversity_table.tail(5).index,inplace=True) # drop last n rows
  667. alpha_diversity_table.at["09.0808", "faith_pd"] = (alpha_diversity_table.at["09.0808.1", "faith_pd"] + alpha_diversity_table.at["09.0808.2", "faith_pd"]) /2
  670. #Load Simpson's Diversity
  671. simpson_table = pd.read_csv("species-simpson-table/alpha-diversity.tsv", sep = "\t",
  672. index_col=0)
  673. simpson_table.drop(simpson_table.tail(5).index,inplace=True) # drop last n rows
  674. simpson_table.at["09.0808", "simpson"] = (simpson_table.at["09.0808.1", "simpson"] + simpson_table.at["09.0808.2", "simpson"]) /2
  677. #Combine all four diveristy metrics
  678. alpha_diversity_table = alpha_diversity_table.join(unifrac_distance)
  679. alpha_diversity_table["BC_Dissimilarity"] = species_table["BC_Dissimilarity"]
  680. alpha_diversity_table = alpha_diversity_table.join(simpson_table)
  682. #Obtain Days
  683. # alpha_diversity_table["Player"]= [i.split(".")[0] for i in alpha_diversity_table.index.tolist()]
  685. # alpha_diversity_table["Date"] = [datetime.datetime.strptime(i.split(".")[1] + "2022", "%m%d%Y").date() for i in alpha_diversity_table.index.tolist()]
  686. # alpha_diversity_table["Date"] = pd.to_datetime(alpha_diversity_table["Date"])
  688. # alpha_diversity_table["Days"] = alpha_diversity_table.groupby("Player")["Date"].apply(lambda x: (x - x.iloc[0]).dt.days)
  689. # print(alpha_diversity_table)
  690. # alpha_diversity_table["Days"] = [i.days for i in alpha_diversity_table["Days"]]
  692. alpha_diversity_table["Player"] = [i.split(".")[0] for i in alpha_diversity_table.index.tolist()]
  693. alpha_diversity_table["Date"] = pd.to_datetime(alpha_diversity_table.index.str.split('.').str[1] + '2022', format='%m%d%Y')
  694. alpha_diversity_table['Days'] = alpha_diversity_table.groupby('Player')['Date'].transform(lambda x: (x - x.min()).dt.days)
  695. print(alpha_diversity_table)
  696. alpha_master_data = master_data.join(alpha_diversity_table)
  697. alpha_master_data.index.name = "Sample_ID"
  698. alpha_master_data.to_csv("../diversity_master.tsv", sep ="\t")
  699. #print(alpha_master_data)
  701. # %% [markdown]
  702. # ## Perform same process for different taxa levels
  704. # %%
  705. def pct_change_manual(df, colnames):
  706. n_rows = df.shape[0]
  707. for i in range(n_rows-1,-1,-1):
  708. for col in colnames:
  710. if i ==0:
  711. df.at[df.index.values[i], col] = np.nan
  712. elif df.at[df.index.values[i], "TWO FECAL SAMPLES FOLLOWING HEAD IMPACT?"] == 1:
  713. ##Two days prior
  714. if np.isnan((df.at[df.index.values[i-2], col])):
  715. df.at[df.index.values[i], col] = np.nan
  716. else:
  717. df.at[df.index.values[i], col] = (df.at[df.index.values[i], col] - df.at[df.index.values[i-2], col])/(df.at[df.index.values[i-2], col] + 0.01)
  718. else:
  719. if np.isnan((df.at[df.index.values[i-1], col])):
  720. df.at[df.index.values[i], col] = np.nan
  721. else:
  722. df.at[df.index.values[i], col] = (df.at[df.index.values[i], col] - df.at[df.index.values[i-1], col])/(df.at[df.index.values[i-1], col] + 0.01)
  723. return df
  725. # %% [markdown]
  726. # # Processing taxa data for R -- Relative Data Tables (0-1 normalized)
  728. # %%
  730. order_table = pd.read_csv("Relative_Tables/order-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  731. order_table = order_table.T
  732. order_table.index.names = ["Sample_ID"]
  733. order_table.drop(order_table.tail(5).index,inplace=True) # drop last n rows
  734. orders = order_table.columns
  735. for i in range(order_table.shape[1]):
  736. column_name = order_table.columns[i]
  737. if "o__" in column_name:
  738. order_table.rename(columns = {column_name:column_name.split("o__", 1)[1]}, inplace= True)
  740. order_table.loc["09.0808"] = (order_table.loc["09.0808.1"] + order_table.loc["09.0808.2"]) /2
  741. order_table.to_csv("../processed_data/order_processed.tsv", sep ="\t")
  745. genus_table = pd.read_csv("Relative_Tables/genus-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  746. genus_table = genus_table.T
  747. genus_table.index.names = ["Sample_ID"]
  748. genus_table.drop(genus_table.tail(5).index,inplace=True) # drop last n rows
  749. for i in range(genus_table.shape[1]):
  750. column_name = genus_table.columns[i]
  751. if "g__" in column_name:
  752. genus_table.rename(columns = {column_name:column_name.split("g__", 1)[1]}, inplace= True)
  754. genus_table.loc["09.0808"] = (genus_table.loc["09.0808.1"] + genus_table.loc["09.0808.2"]) /2
  756. genus_table.to_csv("../processed_data/genus_processed.tsv", sep ="\t")
  760. family_table = pd.read_csv("Relative_Tables/family-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  761. family_table = family_table.T
  762. family_table.index.names = ["Sample_ID"]
  763. family_table.drop(family_table.tail(5).index,inplace=True) # drop last n rows
  764. for i in range(family_table.shape[1]):
  765. column_name = family_table.columns[i]
  766. if "f__" in column_name:
  767. family_table.rename(columns = {column_name:column_name.split("f__", 1)[1]}, inplace= True)
  770. family_table.loc["09.0808"] = (family_table.loc["09.0808.1"] + family_table.loc["09.0808.2"]) /2
  772. family_table.to_csv("../processed_data/family_processed.tsv", sep ="\t")
  774. species_table = pd.read_csv("Relative_Tables/species-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  775. species_table = species_table.T
  776. species_table.index.names = ["Sample_ID"]
  777. species_table.drop(species_table.tail(5).index,inplace=True) # drop last n rows
  778. for i in range(species_table.shape[1]):
  779. column_name = species_table.columns[i]
  780. if "s__" in column_name:
  781. species_table.rename(columns = {column_name:column_name.split("s__", 1)[1]}, inplace= True)
  784. species_table.loc["09.0808"] = (species_table.loc["09.0808.1"] + species_table.loc["09.0808.2"]) /2
  786. species_table.to_csv("../processed_data/species_processed.tsv", sep ="\t")
  789. # %% [markdown]
  790. # # Processing Taxa data for R -- Raw Counts (for CLR transform)
  792. # %%
  793. order_table = pd.read_csv("RawCount_Tables/order-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  794. order_table = order_table.T
  795. order_table.index.names = ["Sample_ID"]
  796. order_table.drop(order_table.tail(5).index,inplace=True) # drop last n rows
  797. orders = order_table.columns
  798. for i in range(order_table.shape[1]):
  799. column_name = order_table.columns[i]
  800. if "o__" in column_name:
  801. order_table.rename(columns = {column_name:column_name.split("o__", 1)[1]}, inplace= True)
  803. order_table.loc["09.0808"] = (order_table.loc["09.0808.1"] + order_table.loc["09.0808.2"]) /2
  804. order_table.to_csv("../processed_data/CLR/order_processed.tsv", sep ="\t")
  808. genus_table = pd.read_csv("RawCount_Tables/genus-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  809. genus_table = genus_table.T
  810. genus_table.index.names = ["Sample_ID"]
  811. genus_table.drop(genus_table.tail(5).index,inplace=True) # drop last n rows
  812. for i in range(genus_table.shape[1]):
  813. column_name = genus_table.columns[i]
  814. if "g__" in column_name:
  815. genus_table.rename(columns = {column_name:column_name.split("g__", 1)[1]}, inplace= True)
  817. genus_table.loc["09.0808"] = (genus_table.loc["09.0808.1"] + genus_table.loc["09.0808.2"]) /2
  819. genus_table.to_csv("../processed_data/CLR/genus_processed.tsv", sep ="\t")
  823. family_table = pd.read_csv("RawCount_Tables/family-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  824. family_table = family_table.T
  825. family_table.index.names = ["Sample_ID"]
  826. family_table.drop(family_table.tail(5).index,inplace=True) # drop last n rows
  827. for i in range(family_table.shape[1]):
  828. column_name = family_table.columns[i]
  829. if "f__" in column_name:
  830. family_table.rename(columns = {column_name:column_name.split("f__", 1)[1]}, inplace= True)
  833. family_table.loc["09.0808"] = (family_table.loc["09.0808.1"] + family_table.loc["09.0808.2"]) /2
  835. family_table.to_csv("../processed_data/CLR/family_processed.tsv", sep ="\t")
  837. species_table = pd.read_csv("RawCount_Tables/species-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  838. species_table = species_table.T
  839. species_table.index.names = ["Sample_ID"]
  840. species_table.drop(species_table.tail(5).index,inplace=True) # drop last n rows
  841. for i in range(species_table.shape[1]):
  842. column_name = species_table.columns[i]
  843. if "s__" in column_name:
  844. species_table.rename(columns = {column_name:column_name.split("s__", 1)[1]}, inplace= True)
  847. species_table.loc["09.0808"] = (species_table.loc["09.0808.1"] + species_table.loc["09.0808.2"]) /2
  849. species_table.to_csv("../processed_data/CLR/species_processed.tsv", sep ="\t")
  852. # %% [markdown]
  853. # # Acute Changes after Head Impact
  855. # %%
  856. order_master_no_9 = pd.read_csv("../processed_data/order_master_without_9.tsv", sep ="\t", dtype={"Sample_ID":"str"})
  857. order_master_no_9 = order_master_no_9.set_index(order_master_no_9.columns[0])
  859. genus_master_no_9 = pd.read_csv("../processed_data/genus_master_without_9.tsv", sep ="\t", dtype={"Sample_ID":"str"})
  860. genus_master_no_9 = genus_master_no_9.set_index(genus_master_no_9.columns[0])
  862. family_master_no_9 = pd.read_csv("../processed_data/family_master_without_9.tsv", sep ="\t", dtype={"Sample_ID":"str"})
  863. family_master_no_9 = family_master_no_9.set_index(family_master_no_9.columns[0])
  865. species_master_no_9 = pd.read_csv("../processed_data/species_master_without_9.tsv", sep ="\t", dtype={"Sample_ID":"str"})
  866. species_master_no_9 = species_master_no_9.set_index(species_master_no_9.columns[0])
  870. # %%
  871. def __convert_to_block_df(
  872. a,
  873. y_col: str = None,
  874. group_col: str = None,
  875. block_col: str = None,
  876. melted: bool = False) -> pd.DataFrame:
  877. # TODO: refactor conversion of block data to DataFrame
  878. if melted and not all([i is not None for i in [block_col, group_col, y_col]]):
  879. raise ValueError(
  880. '`block_col`, `group_col`, `y_col` should be explicitly specified if using melted data')
  882. if isinstance(a, pd.DataFrame) and not melted:
  883. x = a.copy(deep=True)
  884. group_col = 'groups'
  885. block_col = 'blocks'
  886. y_col = 'y'
  887. x.columns.name = group_col
  888. x.index.name = block_col
  889. x = x.reset_index().melt(id_vars=block_col, var_name=group_col, value_name=y_col)
  891. elif isinstance(a, pd.DataFrame) and melted:
  892. x = pd.DataFrame.from_dict({'groups': a[group_col],
  893. 'blocks': a[block_col],
  894. 'y': a[y_col]})
  896. elif not isinstance(a, pd.DataFrame):
  897. x = np.array(a)
  898. x = pd.DataFrame(x, index=np.arange(
  899. x.shape[0]), columns=np.arange(x.shape[1]))
  901. if not melted:
  902. group_col = 'groups'
  903. block_col = 'blocks'
  904. y_col = 'y'
  905. x.columns.name = group_col
  906. x.index.name = block_col
  907. x = x.reset_index().melt(id_vars=block_col, var_name=group_col, value_name=y_col)
  909. else:
  910. x.rename(columns={group_col: 'groups',
  911. block_col: 'blocks', y_col: 'y'}, inplace=True)
  913. group_col = 'groups'
  914. block_col = 'blocks'
  915. y_col = 'y'
  917. return x, y_col, group_col, block_col
  920. def posthoc_nemenyi_friedman_q(
  921. a: Union[list, np.ndarray, pd.DataFrame],
  922. y_col: str = None,
  923. block_col: str = None,
  924. group_col: str = None,
  925. melted: bool = False,
  926. sort: bool = False) -> pd.DataFrame:
  927. '''Calculate pairwise comparisons using Nemenyi post hoc test for
  928. unreplicated blocked data. This test is usually conducted post hoc if
  929. significant results of the Friedman's test are obtained. The statistics
  930. refer to upper quantiles of the studentized range distribution (Tukey) [1]_,
  931. [2]_, [3]_.
  933. Parameters
  934. ----------
  935. a : array_like or pandas DataFrame object
  936. An array, any object exposing the array interface or a pandas
  937. DataFrame.
  939. If `melted` is set to False (default), `a` is a typical matrix of
  940. block design, i.e. rows are blocks, and columns are groups. In this
  941. case you do not need to specify col arguments.
  943. If `a` is an array and `melted` is set to True,
  944. y_col, block_col and group_col must specify the indices of columns
  945. containing elements of correspondary type.
  947. If `a` is a Pandas DataFrame and `melted` is set to True,
  948. y_col, block_col and group_col must specify columns names (strings).
  950. y_col : str or int
  951. Must be specified if `a` is a pandas DataFrame object.
  952. Name of the column that contains y data.
  954. block_col : str or int
  955. Must be specified if `a` is a pandas DataFrame object.
  956. Name of the column that contains blocking factor values.
  958. group_col : str or int
  959. Must be specified if `a` is a pandas DataFrame object.
  960. Name of the column that contains treatment (group) factor values.
  962. melted : bool, optional
  963. Specifies if data are given as melted columns "y", "blocks", and
  964. "groups".
  966. sort : bool, optional
  967. If True, sort data by block and group columns.
  969. Returns
  970. -------
  971. result : pandas.DataFrame
  972. P values.
  974. Notes
  975. -----
  976. A one-way ANOVA with repeated measures that is also referred to as ANOVA
  977. with unreplicated block design can also be conducted via Friedman's
  978. test. The consequent post hoc pairwise multiple comparison test
  979. according to Nemenyi is conducted with this function.
  981. This function does not test for ties.
  983. References
  984. ----------
  985. .. [1] J. Demsar (2006), Statistical comparisons of classifiers over
  986. multiple data sets, Journal of Machine Learning Research, 7, 1-30.
  988. .. [2] P. Nemenyi (1963) Distribution-free Multiple Comparisons. Ph.D.
  989. thesis, Princeton University.
  991. .. [3] L. Sachs (1997), Angewandte Statistik. Berlin: Springer.
  992. Pages: 668-675.
  994. Examples
  995. --------
  996. >>> # Non-melted case, x is a block design matrix, i.e. rows are blocks
  997. >>> # and columns are groups.
  998. >>> x = np.array([[31,27,24],[31,28,31],[45,29,46],[21,18,48],[42,36,46],[32,17,40]])
  999. >>> sp.posthoc_nemenyi_friedman(x)
  1000. '''
  1001. def compare_stats(i, j):
  1002. dif = np.abs(R[groups[i]] - R[groups[j]])
  1003. qval = dif / np.sqrt(k * (k + 1.) * 0.975 / (6. * n))
  1004. return qval
  1006. x, _y_col, _group_col, _block_col = __convert_to_block_df(
  1007. a, y_col, group_col, block_col, melted)
  1008. x = x.sort_values(by=[_group_col, _block_col],
  1009. ascending=True) if sort else x
  1010. x.dropna(inplace=True)
  1012. groups = x[_group_col].unique()
  1013. k = groups.size
  1014. n = x[_block_col].unique().size
  1016. x['mat'] = x.groupby(_block_col)[_y_col].rank()
  1017. R = x.groupby(_group_col)['mat'].mean()
  1018. vs = np.zeros((k, k))
  1019. combs = combinations(range(k), 2)
  1021. tri_upper = np.triu_indices(vs.shape[0], 1)
  1022. tri_lower = np.tril_indices(vs.shape[0], -1)
  1023. vs[:, :] = 0
  1025. for i, j in combs:
  1026. vs[i, j] = compare_stats(i, j)
  1027. qs = vs
  1029. return pd.DataFrame(qs, index=groups, columns=groups)
  1033. # %%
  1034. diversity_master_no_9 = pd.read_csv("../diversity_master_without_9.tsv", sep ="\t", dtype={"Sample_ID":"str"})
  1035. diversity_master_no_9 = diversity_master_no_9.set_index(diversity_master_no_9.columns[0])
  1036. #Data Slicing
  1039. def acute_changes(df, metric, axis_label):
  1040. df_new = df
  1041. df_new = df_new[df_new['BC_Dissimilarity'].notna()]
  1043. #threshold = df_new["Impact load sustained 0-24 hours prior"].quantile(0.75)
  1044. #df_new["Sustained_Impact"] = np.where(df_new["Impact load sustained 0-24 hours prior"]> threshold, 1,0)
  1046. threshold = df_new["Impact.load.sustained.0.24.hours.prior"].quantile(0.75)
  1047. print(threshold)
  1048. df_new["Sustained_Impact"] = np.where(df_new["Impact.load.sustained.0.24.hours.prior"]> threshold, 1,0)
  1052. global data_slices
  1053. data_slices = pd.DataFrame(columns=["Day of Impact", "1 day after", "2 days after", "3 days after", "Player"])
  1055. df_new.groupby("Player").apply(lambda df: slice_data(df, metric))
  1058. data_slices_long = data_slices
  1059. data_slices_long["ID"] = data_slices_long.index.values
  1060. data_slices_long = pd.melt(data_slices_long, value_vars=["Day of Impact", "1 day after", "2 days after", "3 days after"],
  1061. var_name="Time", id_vars= ["ID", "Player"], value_name="Value")
  1062. data_slices_long['Player'] = data_slices_long['Player'].astype(str)
  1064. print(data_slices_long)
  1066. boxplot = sns.boxplot(x='Time',y='Value',data=data_slices_long, color = "white", palette="Greys", width=0.6,showfliers = False)
  1067. boxplot = sns.stripplot(x='Time',y='Value',hue='Player',data=data_slices_long, s= 6, jitter =False)
  1068. boxplot = sns.lineplot(
  1069. data=data_slices_long, x="Time", y="Value", units="ID",sort =False,
  1070. color=".6", estimator=None,alpha=0.7
  1071. )
  1072. boxplot.set_xlabel("Hours since Head Impact",fontsize=15)
  1074. boxplot.set_xticklabels(["0-24", "24-48", "48-72", "72-96"])
  1076. boxplot.set_ylabel(axis_label, fontsize=15)
  1077. plt.locator_params(axis='y', nbins=4)
  1078. boxplot.xaxis.set_tick_params(labelsize=12)
  1079. boxplot.yaxis.set_tick_params(labelsize=12)
  1081. handles, previous_labels = boxplot.get_legend_handles_labels()
  1083. lgd = plt.legend(bbox_to_anchor=(1.2, 0.5), loc='center right', borderaxespad=0, title = "Player",
  1084. fontsize="12",handles=handles, labels=["1","4", "5", "8", "9"])
  1085. lgd.get_title().set_fontsize('15')
  1086. boxplot.figure.savefig("../Figures/" + metric + ".png", dpi = 300, bbox_inches = "tight")
  1087. plt.cla()
  1089. data_slices.to_csv(metric + "slices.tsv", sep = "\t")
  1090. (stat, p) = stats.friedmanchisquare(data_slices["Day of Impact"],
  1091. data_slices["1 day after"], data_slices["2 days after"],
  1092. data_slices["3 days after"])
  1093. nemenyi_statistic = posthoc_nemenyi_friedman_q(data_slices[["Day of Impact", "1 day after", "2 days after", "3 days after"]].values)
  1095. with open("../acute_results/" + metric + ".txt", "w") as f:
  1096. print(stat, file= f)
  1097. print(p, file = f)
  1098. print(scikit_posthocs.posthoc_nemenyi_friedman(data_slices.loc[:,["Day of Impact", "1 day after", "2 days after", "3 days after"]]), file = f)
  1099. print("Nemenyi Test Statistic:", nemenyi_statistic)
  1100. def slice_data(df, metric):
  1102. for i in range(1, (df.shape[0]-3)):
  1104. if (df.at[df.index.values[i], "Sustained_Impact"] == 1) & (df.at[df.index.values[i-1], "Sustained_Impact"] == 0):
  1105. #check if days are continuous
  1106. range_dates = [*range(int(df.index.values[i-1].split(".")[1]), int(df.index.values[i-1].split(".")[1])+5 )]
  1107. #print(range_dates)
  1108. actual_dates = [int(k.split(".")[1]) for k in df.index.values[i-1:i+4]]
  1109. #print(actual_dates)
  1110. if np.array_equal(range_dates, actual_dates):
  1111. #if True:
  1113. #check for no impacts in i+1, i+2, i+3, i+4, i+5 days
  1114. if max(df.loc[df.index.values[i+1:i+4], 'Sustained_Impact'] ==0):
  1116. Player = df.index.values[0].split(".")[0]
  1117. #values_to_add = [ df.at[df.index.values[i-1], "simpson"],
  1118. # (df.at[df.index.values[i+1], "simpson"] +df.at[df.index.values[i+2], "simpson"])/2,
  1119. # ((df.at[df.index.values[i+3], "simpson"] +df.at[df.index.values[i+4], "simpson"])/2),
  1120. # df.at[df.index.values[i+5], "simpson"] ]
  1121. values_to_add = [ df.at[df.index.values[i], metric],
  1122. df.at[df.index.values[i+1], metric],df.at[df.index.values[i+2], metric],
  1123. df.at[df.index.values[i+3], metric] , Player]
  1124. data_slices.loc[len(data_slices.index)] = values_to_add
  1126. #print(alpha_master_data)
  1127. acute_changes(diversity_master_no_9, metric="BC_Dissimilarity", axis_label="Bray Curtis Dissimilarity")
  1128. acute_changes(diversity_master_no_9, metric="faith_pd", axis_label="Faith's Phylogenetic Diversity")
  1129. #acute_changes(alpha_master_data, metric="simpson", axis_label="Simpson's Diversity")
  1130. #acute_changes(alpha_master_data, metric="UniFrac_Distance", axis_label="Unifrac Distance")
  1133. # acute_changes(df=order_master_no_9, metric = "Bacteroidales", axis_label = "Bacteroidales")
  1134. # acute_changes(df=order_master_no_9, metric = "Bifidobacteriales", axis_label = "Bifidobacteriales")
  1135. # acute_changes(df=order_master_no_9, metric = "Lactobacillales", axis_label = "Lactobacillales")
  1136. # acute_changes(df=order_master_no_9, metric = "Burkholderiales", axis_label ="Burkholderiales")
  1137. # acute_changes(df=order_master_no_9, metric = "Enterobacterales", axis_label= "Enterobacterales")
  1138. # acute_changes(df=order_master_no_9, metric = "Pseudomonadales", axis_label = "Pseudomonadales")
  1139. # acute_changes(df=order_master_no_9, metric = "Campylobacterales", axis_label = "Campylobacterales")
  1140. # acute_changes(df=order_master_no_9, metric = "Clostridiales", axis_label = "Clostridiales")
  1141. # acute_changes(df=order_master_no_9, metric = "Verrucomicrobiales", axis_label = "Verrucomicrobiales")
  1144. # acute_changes(df=family_master_no_9, metric = "Micrococcaceae", axis_label = "Micrococcaceae")
  1145. # acute_changes(df=family_master_no_9, metric = "Prevotellaceae", axis_label = "Prevotellaceae")
  1146. # acute_changes(df=family_master_no_9, metric = "Leuconostocaceae", axis_label = "Leuconostocaceae")
  1147. # acute_changes(df=family_master_no_9, metric = "Streptococcaceae", axis_label = "Streptococcaceae")
  1148. # acute_changes(df=family_master_no_9, metric = "Peptococcaceae", axis_label = "Peptococcaceae")
  1149. # acute_changes(df=family_master_no_9, metric = "Lachnospiraceae", axis_label= "Lachnospiraceae")
  1150. # acute_changes(df=family_master_no_9, metric = "Bacteroidaceae", axis_label = "Bacteroidaceae")
  1151. # acute_changes(df=family_master_no_9, metric = "Rikenellaceae", axis_label = "Rikenellaceae")
  1152. # acute_changes(df=family_master_no_9, metric = "Ruminococcaceae", axis_label = "Ruminococcaceae")
  1153. # acute_changes(df=family_master_no_9, metric = "Clostridiaceae", axis_label = "Clostridiaceae")
  1154. # acute_changes(df=family_master_no_9, metric = "Bifidobacteriaceae", axis_label ="Bifidobacteriaceae")
  1156. # acute_changes(df=genus_master_no_9, metric = "Prevotella", axis_label = "Prevotella")
  1157. # acute_changes(df=genus_master_no_9, metric = "Weissella", axis_label = "Weissella")
  1158. # acute_changes(df=genus_master_no_9, metric = "Lactococcus", axis_label = "Lactococcus")
  1159. # acute_changes(df=genus_master_no_9, metric = "Anaerostipes", axis_label = "Anaerostipes")
  1160. # acute_changes(df=genus_master_no_9, metric = "Lactobacillus", axis_label = "Lactobacillus")
  1161. # acute_changes(df=genus_master_no_9, metric = "Ruminococcus", axis_label = "Ruminococcus")
  1162. # acute_changes(df=genus_master_no_9, metric = "Bifidobacterium", axis_label = "Bifidobacterium")
  1163. # acute_changes(df=genus_master_no_9, metric = "Bacteroides", axis_label = "Bacteroides")
  1165. # acute_changes(df=species_master_no_9, metric = "Anaerostipes_hadrus", axis_label = "Anaerostipes_hadrus")
  1168. # %% [markdown]
  1169. # # Pre Mid Post Analysis
  1171. # %%
  1172. # Pre Mid Post Analysis
  1173. def perform_pmp(dataset, colname, axislabel, skip_first = True):
  1175. global pre_mid_post
  1176. pre_mid_post = pd.DataFrame(columns=["Pre", "Mid", "Post", "Player"])
  1177. dataset.loc[:,[colname, "Player"] ].groupby("Player").apply(lambda x: find_pre_mid_post(x, colname, skip_first) )
  1179. pre_mid_post = pre_mid_post.loc[3:,["Pre", "Mid", "Post", "Player"]]
  1180. print(pre_mid_post)
  1181. pre_mid_post.to_csv(colname + "pre_mid_post.tsv", sep = "\t")
  1183. pre_mid_post_long = pre_mid_post
  1184. pre_mid_post_long["ID"] = pre_mid_post_long.index.values
  1185. pre_mid_post_long = pd.melt(pre_mid_post_long, value_vars=["Pre","Mid","Post"],
  1186. var_name="Time", id_vars= ["ID", "Player"], value_name="Value")
  1187. print(pre_mid_post_long)
  1188. pre_mid_post_long['Player'] = pre_mid_post_long['Player'].astype(str)
  1189. boxplot = sns.boxplot(x='Time',y='Value',data=pre_mid_post_long, palette="Greys", width=0.5, showfliers=False)
  1190. boxplot = sns.stripplot(x='Time',y='Value',hue='Player',data=pre_mid_post_long, s= 6, jitter =False)
  1191. boxplot = sns.lineplot(
  1192. data=pre_mid_post_long, x="Time", y="Value", units="ID",sort =False,
  1193. color=".7", estimator=None,alpha=0.7
  1194. )
  1195. boxplot.set_xlabel("Collection Period",fontsize=15)
  1197. boxplot.set_xticklabels(["Early", "Middle", "Late"])
  1199. boxplot.set_ylabel(axislabel, fontsize=15)
  1200. plt.locator_params(axis='y', nbins=4)
  1201. boxplot.xaxis.set_tick_params(labelsize=12)
  1202. boxplot.yaxis.set_tick_params(labelsize=12)
  1203. handles, previous_labels = boxplot.get_legend_handles_labels()
  1205. lgd = plt.legend(bbox_to_anchor=(1.2, 0.5), loc='center right', borderaxespad=0, title = "Player",
  1206. fontsize="12", handles = handles, labels = ["1","4","5","8","9","16"])
  1207. lgd.get_title().set_fontsize('15')
  1208. boxplot.figure.savefig("../PreMidPost_Figures/" + colname + "_pre_mid_post.png", dpi =300, bbox_inches = "tight")
  1209. plt.cla()
  1210. #(stat, p) = stats.kruskal(pre_mid_post_BC["Pre"], pre_mid_post_BC["Mid"],pre_mid_post_BC["Post"])
  1211. #print(stat)
  1212. #print(p)
  1213. nemenyi_statistic = posthoc_nemenyi_friedman_q(pre_mid_post.loc[:, ["Pre", "Mid", "Post"]])
  1214. with open("../PreMidPost/" + colname + ".txt", "w") as f:
  1215. (stat, p) = stats.friedmanchisquare( pre_mid_post["Pre"], pre_mid_post["Mid"],pre_mid_post["Post"])
  1216. print(stat, file =f)
  1217. print(p, file = f)
  1219. print(scikit_posthocs.posthoc_nemenyi_friedman(pre_mid_post.loc[:, ["Pre", "Mid", "Post"]]), file =f)
  1220. print("Nemenyi Test Statistic:", nemenyi_statistic)
  1223. def find_pre_mid_post(df, colname, skip_first):
  1224. nrows = df.shape[0]
  1225. center = math.floor(nrows/2)
  1226. #pre_mid_post.loc[len(pre_mid_post)] = [df.loc[df.index.values[1:4], colname].mean(), df.loc[df.index.values[center-1:center+2], colname].mean(),df.loc[df.index.values[nrows-3:nrows], colname].mean()]
  1227. #pre_mid_post_BC.loc[len(pre_mid_post_BC)] = [df.at[df.index.values[0], colname], df.at[df.index.values[center-1], colname],df.at[df.index.values[nrows-4], colname]]
  1228. Player = df.index.values[0].split(".")[0]
  1229. if skip_first:
  1230. if (df.index.values[0].startswith("01.")| df.index.values[0].startswith("04.")):
  1231. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[1], colname], df.at[df.index.values[center-1], colname],df.at[df.index.values[nrows-5], colname], Player]
  1232. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[2], colname], df.at[df.index.values[center], colname],df.at[df.index.values[nrows-4], colname], Player]
  1233. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[3], colname], df.at[df.index.values[center+1], colname],df.at[df.index.values[nrows-3], colname], Player]
  1235. else:
  1236. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[1], colname], df.at[df.index.values[center-1], colname],df.at[df.index.values[nrows-3], colname], Player]
  1237. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[2], colname], df.at[df.index.values[center], colname],df.at[df.index.values[nrows-2], colname], Player]
  1238. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[3], colname], df.at[df.index.values[center+1], colname],df.at[df.index.values[nrows-1], colname], Player]
  1239. else:
  1240. if (df.index.values[0].startswith("01.")| df.index.values[0].startswith("04.")):
  1241. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[0], colname], df.at[df.index.values[center-1], colname],df.at[df.index.values[nrows-5], colname], Player]
  1242. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[1], colname], df.at[df.index.values[center], colname],df.at[df.index.values[nrows-4], colname], Player]
  1243. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[2], colname], df.at[df.index.values[center+1], colname],df.at[df.index.values[nrows-3], colname], Player]
  1245. else:
  1246. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[0], colname], df.at[df.index.values[center-1], colname],df.at[df.index.values[nrows-3], colname], Player]
  1247. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[1], colname], df.at[df.index.values[center], colname],df.at[df.index.values[nrows-2], colname], Player]
  1248. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[2], colname], df.at[df.index.values[center+1], colname],df.at[df.index.values[nrows-1], colname], Player]
  1250. perform_pmp(diversity_master_no_9, "BC_Dissimilarity", "Bray Curtis Dissimilarity", True)
  1252. perform_pmp(diversity_master_no_9, "faith_pd", "Faith's Phylogenetic Diversity", True)
  1255. # acute_changes(df=order_master_no_9, metric = "Bacteroidales", axis_label = "Bacteroidales")
  1256. # acute_changes(df=order_master_no_9, metric = "Bifidobacteriales", axis_label = "Bifidobacteriales")
  1257. # acute_changes(df=order_master_no_9, metric = "Lactobacillales", axis_label = "Lactobacillales")
  1258. ##perform_pmp(order_master_no_9, "Coriobacteriales", "Coriobacteriales", False)
  1260. ##perform_pmp(order_master_no_9, "Burkholderiales", "Burkholderiales", False)
  1261. # acute_changes(df=order_master_no_9, metric = "Enterobacterales", axis_label= "Enterobacterales")
  1262. # acute_changes(df=order_master_no_9, metric = "Pseudomonadales", axis_label = "Pseudomonadales")
  1263. # acute_changes(df=order_master_no_9, metric = "Campylobacterales", axis_label = "Campylobacterales")
  1264. # acute_changes(df=order_master_no_9, metric = "Clostridiales", axis_label = "Clostridiales")
  1265. # acute_changes(df=order_master_no_9, metric = "Verrucomicrobiales", axis_label = "Verrucomicrobiales")
  1268. # acute_changes(df=family_master_no_9, metric = "Micrococcaceae", axis_label = "Micrococcaceae")
  1269. ##perform_pmp(family_master_no_9, "Prevotellaceae", "Prevotellaceae", False)
  1270. # acute_changes(df=family_master_no_9, metric = "Leuconostocaceae", axis_label = "Leuconostocaceae")
  1271. # acute_changes(df=family_master_no_9, metric = "Streptococcaceae", axis_label = "Streptococcaceae")
  1272. ##perform_pmp(family_master_no_9, "Peptococcaceae", "Peptococcaceae", False)
  1273. ##perform_pmp(family_master_no_9, "Lachnospiraceae", "Lachnospiraceae", False)
  1274. # acute_changes(df=family_master_no_9, metric = "Bacteroidaceae", axis_label = "Bacteroidaceae")
  1275. # acute_changes(df=family_master_no_9, metric = "Rikenellaceae", axis_label = "Rikenellaceae")
  1276. ##perform_pmp(family_master_no_9, "Ruminococcaceae", "Ruminococcaceae")
  1277. # acute_changes(df=family_master_no_9, metric = "Clostridiaceae", axis_label = "Clostridiaceae")
  1278. # acute_changes(df=family_master_no_9, metric = "Bifidobacteriaceae", axis_label ="Bifidobacteriaceae")
  1280. ##perform_pmp(genus_master_no_9, "Prevotella", "Prevotella", False)
  1281. ##perform_pmp(genus_master_no_9, "Weissella", "Weissella", False)
  1282. # acute_changes(df=genus_master_no_9, metric = "Lactococcus", axis_label = "Lactococcus")
  1283. # acute_changes(df=genus_master_no_9, metric = "Anaerostipes", axis_label = "Anaerostipes")
  1284. # acute_changes(df=genus_master_no_9, metric = "Lactobacillus", axis_label = "Lactobacillus")
  1285. ##perform_pmp(genus_master_no_9, "Ruminococcus", "Ruminococcus", False)
  1286. # acute_changes(df=genus_master_no_9, metric = "Bifidobacterium", axis_label = "Bifidobacterium")
  1287. # acute_changes(df=genus_master_no_9, metric = "Bacteroides", axis_label = "Bacteroides")
  1290. # %% [markdown]
  1291. # # Phyla Level Analysis
  1293. # %%
  1294. phyla_table = pd.read_csv("phyla-table.tsv", sep = "\t", skiprows=[0], index_col = 0)
  1295. phyla_table = phyla_table.T
  1296. phyla_table.index.names = ["Sample_ID"]
  1297. #phyla_table.drop(phyla_table.tail(5).index,inplace=True) # drop last n rows
  1298. for i in range(phyla_table.shape[1]):
  1299. column_name = phyla_table.columns[i]
  1300. if "p__" in column_name:
  1301. phyla_table.rename(columns = {column_name:column_name.split("p__", 1)[1]}, inplace= True)
  1303. phyla_table.loc["09.0808"] = (phyla_table.loc["09.0808.1"] + phyla_table.loc["09.0808.2"]) /2
  1304. phyla_table.rename(columns={'d__Bacteria;__': 'Unclassified'}, inplace=True)
  1305. phyla_table.to_csv("../phyla_processed_with_control.tsv", sep ="\t")
  1309. # %%
  1310. phyla = sorted(phyla_table.columns.values.tolist())
  1311. print(phyla)
  1312. phyla_master_no_9 = pd.read_csv("../phyla_master_without_9.tsv", sep ="\t", dtype={"Sample_ID":"str"})
  1313. phyla_master_no_9 = phyla_master_no_9.set_index(phyla_master_no_9.columns[0])
  1315. pre_phyla = pd.DataFrame(columns=np.append(phyla, ["Player"]))
  1317. mid_phyla = pd.DataFrame(columns=np.append(phyla, ["Player"]))
  1318. post_phyla = pd.DataFrame(columns=np.append(phyla, ["Player"]))
  1319. def find_pre_mid_post_phyla(df, list_of_phyla):
  1320. global pre_phyla
  1321. global mid_phyla
  1322. global post_phyla
  1323. nrows = df.shape[0]
  1324. center = math.floor(nrows/2)
  1325. Player = df.index.values[0].split(".")[0]
  1326. if (df.index.values[0].startswith("01.")| df.index.values[0].startswith("04.")):
  1327. pre_phyla = pre_phyla.append(pd.DataFrame(df.loc[df.index.values[1:4],phyla].mean()).T.assign(Player = Player), ignore_index=True)
  1328. mid_phyla = mid_phyla.append(pd.DataFrame(df.loc[df.index.values[center-1:center+2],phyla].mean()).T.assign(Player=Player), ignore_index=True)
  1329. post_phyla = post_phyla.append(pd.DataFrame(df.loc[df.index.values[nrows-5:nrows-2],phyla].mean()).T.assign(Player = Player), ignore_index=True)
  1330. else:
  1331. pre_phyla = pre_phyla.append(pd.DataFrame(df.loc[df.index.values[1:4],phyla].mean()).T.assign(Player = Player), ignore_index=True)
  1332. mid_phyla = mid_phyla.append(pd.DataFrame(df.loc[df.index.values[center-1:center+2],phyla].mean()).T.assign(Player=Player), ignore_index=True)
  1333. post_phyla = post_phyla.append(pd.DataFrame(df.loc[df.index.values[nrows-3:nrows],phyla].mean()).T.assign(Player = Player), ignore_index=True)
  1336. phyla_master_no_9.groupby("Player").apply(lambda x: find_pre_mid_post_phyla(x, phyla) )
  1338. find_pre_mid_post_phyla(phyla_master_no_9, phyla)
  1340. pre_phyla.drop(index = pre_phyla.head(1).index,inplace=True)
  1341. pre_phyla.drop(index = pre_phyla.tail(1).index,inplace=True)
  1342. mid_phyla.drop(index=mid_phyla.head(1).index,inplace=True)
  1343. mid_phyla.drop(index = mid_phyla.tail(1).index,inplace=True)
  1344. post_phyla.drop(index =post_phyla.head(1).index,inplace=True)
  1345. post_phyla.drop(index =post_phyla.tail(1).index,inplace=True)
  1347. fig, ax = plt.subplots(3,6)
  1349. for i,time in enumerate(["Pre", "Mid", "Post"]):
  1350. player_title = ["1", "4", "5", "8", "9", "16"]
  1351. time_label = ["Early", "Middle", "Late"]
  1352. for j, Player in enumerate(["01", "04", "05", "08", "09", "16"]):
  1354. if time == "Pre":
  1355. fracs = pre_phyla.loc[pre_phyla['Player'] == Player, phyla]
  1356. elif time == "Mid":
  1357. fracs = mid_phyla.loc[mid_phyla['Player'] == Player, phyla]
  1358. else:
  1359. fracs = post_phyla.loc[post_phyla['Player'] == Player, phyla]
  1360. ax[i, j].text(0.5, -0.2, player_title[j], transform=ax[i, j].transAxes,
  1361. horizontalalignment='center', verticalalignment='center',
  1362. fontsize =15)
  1363. #ax[i,j].set_title(player_title[j], y=-0.01)
  1364. ax[i,j].pie(fracs, labels=None,
  1365. autopct=None, shadow=False, startangle=90,
  1366. colors=["blue","orange","red", "cyan","pink","purple","green","grey","brown"])
  1367. ax[i,j].set_aspect('equal')
  1368. # label y axis
  1369. if ax[i,j].is_first_col():
  1370. ax[i,j].set_ylabel(time_label[i], fontsize = 15, rotation = 0,labelpad=18)
  1373. fig.text(0.55, -0.01, "Player", ha="center", fontsize=15)
  1375. plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
  1377. lgd = fig.legend(phyla, loc = 'center right',bbox_to_anchor=(1.3, 0.5))
  1378. fig.savefig("../PreMidPost_Figures/PreMidPost_phyla.png", dpi = 300, bbox_inches = "tight")
  1379. plt.close()
  1381. ### Stacked Bar Plot
  1382. fig, ax = plt.subplots(3,1, sharex = True)
  1384. for i,time in enumerate(["Early", "Middle", "Late"]):
  1385. if time == "Early":
  1386. fracs = pre_phyla
  1387. elif time == "Middle":
  1388. fracs = mid_phyla
  1389. else:
  1390. fracs = post_phyla
  1391. fracs.plot.bar(x = "Player", stacked = True, ax=ax[i], legend = False,
  1392. color=["blue","orange","red", "cyan","pink","purple","green","grey","brown"])
  1394. ax[i].set_title(time, fontsize = 15, loc='left')
  1395. ax[i].yaxis.set_tick_params(labelsize=12)
  1396. ax[1].set_ylabel('Relative Abundance', fontsize = 15)
  1398. ax[2].set_xlabel("Player", fontsize =15)
  1399. ax[2].set_xticklabels(["1", "4", "5", "8", "9", "16"], rotation = 0, fontsize =12)
  1400. plt.tight_layout(pad=0.4, w_pad=0.6, h_pad=1.0)
  1401. plt.subplots_adjust(wspace=0, hspace=0.5)
  1403. lgd = fig.legend(phyla, loc = 'center right',bbox_to_anchor=(1.3, 0.5), borderaxespad = 0., title= 'Phyla')
  1404. lgd.get_title().set_fontsize('15')
  1405. fig.savefig("../PreMidPost_Figures/PreMidPost_phyla_bar.png", dpi = 300, bbox_inches = "tight")
  1406. plt.close()
  1410. # Adding metadata for PDF file
  1413. # %%
  1414. phyla_master_no_9_control = pd.read_csv("../phyla_master_without_9_with_control.tsv", sep ="\t", dtype={"Sample_ID":"str"})
  1415. phyla_master_no_9_control = phyla_master_no_9_control.set_index(phyla_master_no_9_control.columns[0])
  1417. phyla_master_no_9_control_players = phyla_master_no_9_control.drop(
  1418. phyla_master_no_9_control.tail(5).index) # drop last n rows
  1420. phyla_master_no_9_control_players["Player"] = [i.split(".")[0] for i in phyla_master_no_9_control_players.index.tolist()]
  1423. first_data_points = phyla_master_no_9_control_players.groupby("Player").first().reset_index().loc[:,np.append(phyla, ["Player"])]
  1424. controls = phyla_master_no_9_control.tail(5)
  1425. controls.sort_index(axis=1, inplace=True)
  1426. controls["Player"] = controls.index.values
  1427. print(controls)
  1429. fig, ax = plt.subplots(1,2, sharey = True)
  1432. first_data_points.plot.bar(x = "Player", stacked = True, ax=ax[0], legend = False,
  1433. color=["blue","orange","red", "brown","pink","purple","green","grey","cyan"])
  1434. controls.plot.bar(x = "Player", stacked = True, ax=ax[1], legend = False,
  1435. color=["blue","orange","red", "brown","pink","purple","green","grey","cyan"])
  1436. ax[0].set_ylabel('Relative Abundance', fontsize = 15)
  1437. ax[0].set_xlabel("Player", fontsize =15)
  1438. ax[1].set_xlabel("Control Samples", fontsize =15)
  1439. ax[0].set_xticklabels(["1", "4", "5", "8", "9", "16"], rotation = 0, fontsize =12)
  1440. ax[1].set_xticklabels(["Mock", "N1", "N2", "P1", "P2"], rotation = 0, fontsize =12)
  1441. ax[0].yaxis.set_tick_params(labelsize=12)
  1442. plt.tight_layout(pad=0.4, w_pad=0.6, h_pad=1.0)
  1443. plt.subplots_adjust(wspace=0.1, hspace=0)
  1444. plt.locator_params(axis='y', nbins=3)
  1445. lgd = fig.legend(phyla, loc = 'center right',bbox_to_anchor=(1.3, 0.5), borderaxespad = 0., title= 'Phyla')
  1446. lgd.get_title().set_fontsize('15')
  1447. fig.savefig("../PreMidPost_Figures/Control_Bar.png", dpi = 300, bbox_inches = "tight")
  1448. plt.close()
  1452. # %%
  1453. phyla = phyla_table.columns
  1454. phyla_master_no_9 = pd.read_csv("../phyla_master_without_9.tsv", sep ="\t", dtype={"Sample_ID":"str"})
  1455. phyla_master_no_9 = phyla_master_no_9.set_index(phyla_master_no_9.columns[0])
  1457. pre_mid_post_phyla = pd.DataFrame(columns=["Pre", "Mid", "Post", "Player", "Phyla"])
  1458. phyla_master_no_9.loc[:,[phyla, "Player"] ].groupby("Player").apply(lambda x: find_pre_mid_post(x, colname, skip_first) )
  1460. pre_mid_post = pre_mid_post.loc[3:,["Pre", "Mid", "Post", "Player", "Phyla"]]
  1462. pre_mid_post_long = pre_mid_post
  1463. pre_mid_post_long["ID"] = pre_mid_post_long.index.values
  1464. pre_mid_post_long = pd.melt(pre_mid_post_long, value_vars=["Pre","Mid","Post"],
  1465. var_name="Time", id_vars= ["ID", "Player"], value_name="Value")
  1466. print(pre_mid_post_long)
  1467. pre_mid_post_long['Player'] = pre_mid_post_long['Player'].astype(str)
  1472. def find_pre_mid_post(df, colname, skip_first):
  1473. nrows = df.shape[0]
  1474. center = math.floor(nrows/2)
  1475. #pre_mid_post.loc[len(pre_mid_post)] = [df.loc[df.index.values[1:4], colname].mean(), df.loc[df.index.values[center-1:center+2], colname].mean(),df.loc[df.index.values[nrows-3:nrows], colname].mean()]
  1476. #pre_mid_post_BC.loc[len(pre_mid_post_BC)] = [df.at[df.index.values[0], colname], df.at[df.index.values[center-1], colname],df.at[df.index.values[nrows-4], colname]]
  1477. Player = df.index.values[0].split(".")[0]
  1480. if (df.index.values[0].startswith("01.")| df.index.values[0].startswith("04.")):
  1481. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[0], colname], df.at[df.index.values[center-1], colname],df.at[df.index.values[nrows-5], colname], Player]
  1482. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[1], colname], df.at[df.index.values[center], colname],df.at[df.index.values[nrows-4], colname], Player]
  1483. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[2], colname], df.at[df.index.values[center+1], colname],df.at[df.index.values[nrows-3], colname], Player]
  1485. else:
  1486. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[0], colname], df.at[df.index.values[center-1], colname],df.at[df.index.values[nrows-3], colname], Player]
  1487. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[1], colname], df.at[df.index.values[center], colname],df.at[df.index.values[nrows-2], colname], Player]
  1488. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[2], colname], df.at[df.index.values[center+1], colname],df.at[df.index.values[nrows-1], colname], Player]
  1491. # %%
  1492. pre_mid_post = pd.DataFrame(columns=["Pre", "Mid", "Post"])
  1494. def find_pre_mid_post(df, colname):
  1495. global pre_mid_post
  1496. nrows = df.shape[0]
  1497. center = math.floor(nrows/2)
  1498. #pre_mid_post.loc[len(pre_mid_post)] = [df.loc[df.index.values[1:4], colname].mean(), df.loc[df.index.values[center-1:center+2], colname].mean(),df.loc[df.index.values[nrows-3:nrows], colname].mean()]
  1499. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[0], colname], df.at[df.index.values[center-1], colname],df.at[df.index.values[nrows-4], colname]]
  1500. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[1], colname], df.at[df.index.values[center], colname],df.at[df.index.values[nrows-3], colname]]
  1501. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[2], colname], df.at[df.index.values[center+1], colname],df.at[df.index.values[nrows-2], colname]]
  1502. pre_mid_post.loc[len(pre_mid_post)] = [df.at[df.index.values[3], colname], df.at[df.index.values[center+2], colname],df.at[df.index.values[nrows-1], colname]]
  1504. #print(alpha_diversity_master_dummy)
  1505. alpha_diversity_master_dummy.loc[:,["faith_pd", "Player"] ].groupby("Player").apply(lambda x: find_pre_mid_post(x, "faith_pd") )
  1507. pre_mid_post = pre_mid_post.loc[4:,["Pre", "Mid", "Post"]]
  1508. print(pre_mid_post)
  1510. pre_mid_post.boxplot().figure.savefig("alpha_pre_mid_post.png", dpi =300, bbox_inches="tight")
  1511. (stat, p) = stats.f_oneway(pre_mid_post["Pre"], pre_mid_post["Mid"],pre_mid_post["Post"])
  1512. print(stat)
  1513. print(p)
  1515. (stat, p) = stats.friedmanchisquare( pre_mid_post["Pre"], pre_mid_post["Mid"],pre_mid_post["Post"])
  1516. print(stat)
  1517. print(p)
  1519. print(stats.wilcoxon(pre_mid_post["Pre"], pre_mid_post["Post"]))
  1522. # %% [markdown]
  1523. # # Metabolomics Analysis
  1525. # %%
  1526. %%bash
  1527. ../FAPROTAX_1.2.6/collapse_table.py \
  1528. -i species-table/feature-table.biom \
  1529. -o func_table.tsv \
  1530. -g ../FAPROTAX_1.2.6/FAPROTAX.txt -v
  1533. # %%
  1534. func_table = pd.read_csv("func_table.tsv", sep = "\t", index_col = 0)
  1535. func_table = func_table.T
  1536. func_table.index.names = ["Sample_ID"]
  1537. func_table.drop(func_table.tail(5).index,inplace=True) # drop last n rows
  1538. func_table.loc["09.0808"] = (func_table.loc["09.0808.1"] + func_table.loc["09.0808.2"]) /2
  1540. func_table.to_csv("../func_processed.tsv",sep ="\t")
  1542. # %% [markdown]
  1543. # # PCoA Code
  1545. # %%
  1547. df_new = species_master_no_9
  1549. threshold = df_new["Impact.load.sustained.0.24.hours.prior"].quantile(0.75)
  1550. df_new["Sustained_Impact"] = np.where(df_new["Impact.load.sustained.0.24.hours.prior"] > threshold, 1, 0)
  1552. species = species_master_no_9.columns.tolist()
  1553. species = species[0:202]
  1555. #print(species)
  1556. sliced_df = pd.DataFrame(columns=np.append(species, ["Player"]))
  1559. def slice_data_for_PCoA(group_df):
  1560. result_df = pd.DataFrame(columns=sliced_df.columns)
  1561. for i in range(1, (group_df.shape[0] - 3)):
  1562. if (group_df.at[group_df.index.values[i], "Sustained_Impact"] == 1) & (group_df.at[group_df.index.values[i - 1], "Sustained_Impact"] == 0):
  1563. # check if days are continuous
  1564. range_dates = [*range(int(group_df.index.values[i - 1].split(".")[1]), int(group_df.index.values[i - 1].split(".")[1]) + 5)]
  1565. actual_dates = [int(k.split(".")[1]) for k in group_df.index.values[i - 1:i + 4]]
  1566. if np.array_equal(range_dates, actual_dates):
  1567. # check for no impacts in i+1, i+2, i+3, i+4, i+5 days
  1568. if max(group_df.loc[group_df.index.values[i + 1:i + 4], 'Sustained_Impact'] == 0):
  1569. Player = group_df.index.values[0].split(".")[0]
  1570. result_df = result_df.append(group_df.loc[group_df.index.values[i - 1:i + 4], species + ["Player"]].assign(Player=Player), ignore_index = True)
  1571. #print(result_df)
  1572. return result_df
  1574. grouped = df_new.groupby("Player")
  1575. sliced_df = pd.concat([slice_data_for_PCoA(group_df) for _, group_df in grouped], ignore_index = True)
  1579. #sliced_df=sliced_df.loc[:,species]
  1580. print(sliced_df)
  1581. #print(pdist(sliced_df))
  1582. matrix = squareform(pdist(sliced_df, metric = "euclidean"))
  1585. print(matrix)
  1586. pcoa_results = pcoa(matrix)
  1587. pcoa_df = pcoa_results.samples[['PC1', 'PC2']]
  1588. print(pcoa_df.shape)
  1589. a = np.array(["Day Pre", "Day of Impact", "1 day after", "2 days after", "3 days after"])
  1590. pcoa_df["Day"] = np.tile(a,int(sliced_df.shape[0]/5))
  1591. plt.scatter(pcoa_df.loc[pcoa_df["Day"] == "Day Pre", "PC1"],
  1592. pcoa_df.loc[pcoa_df["Day"] == "Day Pre", "PC2"],c = 'black')
  1594. plt.scatter(pcoa_df.loc[pcoa_df["Day"] == "Day of Impact", "PC1"],
  1595. pcoa_df.loc[pcoa_df["Day"] == "Day of Impact", "PC2"],c = 'blue')
  1597. plt.scatter(pcoa_df.loc[pcoa_df["Day"] == "1 day after", "PC1"],
  1598. pcoa_df.loc[pcoa_df["Day"] == "1 day after", "PC2"],c = 'red')
  1601. plt.scatter(pcoa_df.loc[pcoa_df["Day"] == "2 days after", "PC1"],
  1602. pcoa_df.loc[pcoa_df["Day"] == "2 days after", "PC2"],c = 'green')
  1604. plt.scatter(pcoa_df.loc[pcoa_df["Day"] == "3 days after", "PC1"],
  1605. pcoa_df.loc[pcoa_df["Day"] == "3 days after", "PC2"],c = 'purple')
  1607. def encircle2(x,y, ax=None, **kw):
  1608. if not ax: ax=plt.gca()
  1609. p = np.c_[x,y]
  1610. mean = np.mean(p, axis=0)
  1611. d = p-mean
  1612. r = np.max(np.sqrt(d[:,0]**2+d[:,1]**2 ))
  1613. circ = plt.Circle(mean, radius=1.05*r,**kw)
  1614. ax.add_patch(circ)
  1616. encircle2(pcoa_df.loc[pcoa_df["Day"] == "Day Pre", "PC1"],
  1617. pcoa_df.loc[pcoa_df["Day"] == "Day Pre", "PC2"], ec="black", fc="none")
  1618. encircle2( pcoa_df.loc[pcoa_df["Day"] == "Day of Impact", "PC1"],
  1619. pcoa_df.loc[pcoa_df["Day"] == "Day of Impact", "PC2"], ec = 'blue',fc="none")
  1620. encircle2(pcoa_df.loc[pcoa_df["Day"] == "1 day after", "PC1"],
  1621. pcoa_df.loc[pcoa_df["Day"] == "1 day after", "PC2"],ec = 'red', fc = "none")
  1622. encircle2(pcoa_df.loc[pcoa_df["Day"] == "2 days after", "PC1"],
  1623. pcoa_df.loc[pcoa_df["Day"] == "2 days after", "PC2"],ec = 'green', fc = "none")
  1624. encircle2(pcoa_df.loc[pcoa_df["Day"] == "3 days after", "PC1"],
  1625. pcoa_df.loc[pcoa_df["Day"] == "3 days after", "PC2"],ec = 'purple', fc="none")
  1627. plt.gca().relim()
  1628. plt.gca().autoscale_view()
  1629. plt.xlabel("PC1")
  1630. plt.ylabel("PC2")
  1631. plt.legend(["Day Pre", "Day of Impact", "1 day after", "2 days after", "3 days after"])
  1632. plt.savefig("PcoA_01.png", dpi = 300, bbox_inches="tight")

TBI Microbiome.ipynb at commit d13d3d7, no license · at the source

Overview

Authors: Zachary J. Pelland1, Aziz Zafar2,3, Ahmet A. Ay2,3, Kenneth Douglas Belanger2
ORCID iDs: Zachary J. Pelland, Kenneth Douglas Belanger
  1. Program in Neuroscience, Colgate University, Hamilton, New York, United States of America
  2. Department of Biology, Colgate University, Hamilton, New York, United States of America
  3. Department of Mathematics, Colgate University, Hamilton, New York, United States of America
Institutions: Colgate University (United States)
Journal: PloS one, volume 21, issue 5, article e0345651
Dates: received 7 October 2025; accepted 9 March 2026; published online 6 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pone.0345651 · PMID 42090386 · PMCID PMC13148679 · OpenAlex W7160403358
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), traumatic brain injury (population), clinical / translational (subfield)
Methods: Statistics
MeSH: Brain Concussion*, Football*, Gastrointestinal Microbiome*, Athletes, Humans, Male, United States, Young Adult (* major topic)
Journal subjects: Biology and Life Sciences, Microbiology, Medical Microbiology, Microbiome, Genetics, Genomics, Microbial Genomics, Psychology, Behavior, Recreation, Sports, Social Sciences, Sports Science, Immunology, Immune Response, Inflammation, Medicine and Health Sciences, Clinical Medicine, Signs and Symptoms, Pharmacology, Drugs, Analgesics, NSAIDs, Pain management, Critical Care and Emergency Medicine, Trauma Medicine, Traumatic Injury, Neurotrauma, Traumatic Brain Injury, Organisms, Bacteria, Gut Bacteria, Ruminococcus, Anatomy, Head, Physiology, Physiological Processes, Sleep
Topic: Traumatic Brain Injury Research (Epidemiology, Medicine), according to OpenAlex
Funding: Colgate University (N/A)
Citations: cited by 1 paper (Europe PMC); 70 references in the paper

Abstract

Non-concussive head impacts (NHIs) are a significant health concern among at-risk groups, including athletes and military personnel. NHIs are hits to the head or head acceleration events (HAEs) that do not generate clinically detectable symptoms and are unlikely to meet diagnostic criteria for mild traumatic brain injury (mTBI). The composition of the gut microbiota influences many aspects of health and wellness and can be altered by TBIs and by brain-related diseases and disorders; however, microbiome alterations have not previously been linked to NHIs. We investigated whether NHIs in a cohort of American football players correlate with acute and long-term changes in the gut microbiome. This study monitored head impact exposure, gut microbiome composition, and a breadth of clinical and behavioral factors in a cohort of collegiate American football players across a competition season. Both short- and long-term changes in the microbiome were analyzed for correlation with head impact events and mathematical modeling was used to examine the contribution of NHIs and other clinical factors to these changes. We observe that NHI exposure correlates with changes in microbial diversity and composition three days following a head impact event. Furthermore, the athletes’ gut microbiomes change significantly across the season, with evidence from mixed-effects modeling indicating that the cumulative effects of NHIs contribute to this change. Our results provide strong evidence for a link between NHIs and changes in the diversity and composition of the gut microbiome. The outcomes of this study emphasize the importance of careful monitoring of head impacts, including those that do not generate clinical symptoms.

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 12 matches between paragraphs and lines of code.

aziz-zafar/TBI-Microbiome

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: d13d3d7c12fa1833180f138096ba1848f150c08e, 27 January 2026
Languages: Jupyter (2), R (2), Python (1)
Size: 8 files, 5 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, 4 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (3 files), NumPy (3 files), pandas (3 files), SciPy (3 files), car (2 files), ggplot2 (2 files), igraph (2 files), lme4 (2 files), lmerTest (2 files), scikit-learn (2 files), seaborn (2 files), tidyverse (2 files), scikit-posthocs (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
6 files
At the source: github.com/aziz-zafar/TBI-Microbiome

The paper's code and data availability statement is in the Data section.

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;
  • 5 scripts, each with its path and the digest of its content;
  • 12 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

Data links

Data Availability

Raw 16S rRNA sequences generated in this study are available via the NCBI BioProject database under accession number PRJNA1111907 (http://www.ncbi.nlm.nih.gov/bioproject/1111907). All metadata and the code used for data processing and statistical analyses are available on this GitHub repository: https://github.com/aziz-zafar/TBI-Microbiome.

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, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 8 MeSH terms, 1 funder, 68 references.

Cite

This paper

Pelland, Z. J., Zafar, A., Ay, A. A., & Belanger, K. D. (2026). Non-concussive head impacts sustained during American football correlate with changes in gut microbiome diversity and composition. PloS one, 21(5), e0345651. https://doi.org/10.1371/journal.pone.0345651

BibTeX

@article{pelland2026non,
author = {Pelland, Zachary J. and Zafar, Aziz and Ay, Ahmet A. and Belanger, Kenneth Douglas},
title = {{Non-concussive head impacts sustained during American football correlate with changes in gut microbiome diversity and composition}},
journal = {PloS one},
year = {2026},
month = may,
volume = {21},
number = {5},
pages = {e0345651},
publisher = {PLOS},
issn = {1932-6203},
doi = {10.1371/journal.pone.0345651},
url = {https://doi.org/10.1371/journal.pone.0345651},
pmid = {42090386},
pmcid = {PMC13148679}
}

RIS

TY - JOUR
AU - Pelland, Zachary J.
AU - Zafar, Aziz
AU - Ay, Ahmet A.
AU - Belanger, Kenneth Douglas
TI - Non-concussive head impacts sustained during American football correlate with changes in gut microbiome diversity and composition
T2 - PloS one
J2 - PLoS One
PY - 2026
DA - 2026/05/06
VL - 21
IS - 5
SP - e0345651
SN - 1932-6203
PB - PLOS
DO - 10.1371/journal.pone.0345651
UR - https://doi.org/10.1371/journal.pone.0345651
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pone.0345651",
"type": "article-journal",
"title": "Non-concussive head impacts sustained during American football correlate with changes in gut microbiome diversity and composition",
"container-title": "PloS one",
"author": [
{
"family": "Pelland",
"given": "Zachary J."
},
{
"family": "Zafar",
"given": "Aziz"
},
{
"family": "Ay",
"given": "Ahmet A."
},
{
"family": "Belanger",
"given": "Kenneth Douglas"
}
],
"container-title-short": "PLoS One",
"volume": "21",
"issue": "5",
"page": "e0345651",
"DOI": "10.1371/journal.pone.0345651",
"PMID": "42090386",
"PMCID": "PMC13148679",
"ISSN": "1932-6203",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pone.0345651",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
6
]
]
}
}

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.1093/braincomms/fcag176 [code]
Tau topography subtypes account for clinical heterogeneity and longitudinal trajectories in early-onset Alzheimer's disease.
Journal: Brain communications
In common: car, lmerTest, lme4, 9 other tools, clinical / translational
[2] doi:10.34133/csbj.0042 [code]
Using Steady-State Visual Evoked Potentials to Characterize Wide-Ranging Retinopathy Linked to <i>CRB1</i>: Implications for Clinical Trials.
Journal: Computational and structural biotechnology journal
In common: car, lmerTest, lme4, 9 other tools, clinical / translational
[3] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: car, igraph, lmerTest, 9 other tools
[4] doi:10.1016/j.nicl.2026.104012 [code]
Structural-functional multilayer brain network properties and outcome of combined repetitive transcranial magnetic stimulation and psychotherapy for obsessive-compulsive disorder.
Journal: NeuroImage. Clinical
In common: car, lmerTest, lme4, 9 other tools
[5] doi:10.1162/imag.a.105 [code]
Right posterior theta reflects human parahippocampal phase resetting by salient cues during goal-directed navigation
Journal: n/a
In common: car, lmerTest, lme4, 8 other tools
[6] doi:10.3389/fnhum.2026.1839961 [code]
Developmental stability of task-rest neural efficiency in youth using a threat and cognitive control task.
Journal: Frontiers in human neuroscience
In common: car, lmerTest, lme4, 8 other tools
[7] doi:10.1093/jbmrpl/ziag077 [code]
Sex- and site-specific reference data for size-invariant properties using multi-stack HRpQCT.
Journal: JBMR plus
In common: car, lmerTest, lme4, 8 other tools
[8] doi:10.64898/2026.03.09.710596 [code]
Infant gut microbiomes contribute to metabolic states that impact brain function
Journal: bioRxiv (preprint)
In common: car, lmerTest, lme4, 8 other tools
[9] doi:10.1126/sciadv.adz6517 [code]
Corticosterone-linked microglial activity underpins sexually dimorphic neuroplasticity after ketamine anesthesia.
Journal: Science advances
In common: scikit-posthocs, lme4, statsmodels, 8 other tools
[10] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: car, igraph, lme4, 8 other tools

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.