1 of 37

The reproducibility of the Crustacean phylogenetic analyses on open source data

Supervisor:

P. Drozdova1�1 Irkutsk State University

Authors:

A. Belyaeva2, 3, A. Zhuravlev3, 4

2 Lomonosov Moscow State University

3 Bioinformatics Institute

4 B.P. Konstantinov Petersburg Nuclear Physics Institute

2 of 37

Introduction

In phylogenetic studies, the problem of reproducibility of data is a major issue. This problem arises due to several factors, including the complexity of the data, the use of different methods and software packages, and the lack of standardization in reporting results. The use of different methods and software packages can lead to different results, which makes it difficult to compare and reproduce results. Additionally, the lack of standardization in reporting results makes it difficult for other researchers to reproduce the same results.

This project examined the issue of phylogenetic analysis reproducibility in a number of studies on freshwater crustaceans based on open data.

3 of 37

Aims and Objectives

Objectives:

​

  • searching for the vault of data (literature review, blastn, NCBI database search due to data errors, study of additional materials)
  • repeating them with modern tools
  • comparison of the results of using several tools at a specific step
  • justification of the non-reproducibility of some steps of the original pipeline
  • creation of “must have” list with essential for this kind of phylogenetic analysis
  • creation of a list of shortcomings, typos, inconsistencies in the articles reviewed
  • creation of “must have” list with essential for replication checkpoints

Aim: The aim of this project was to repeat the phylogenetic research in certain articles and understand about reproducibility of results.

4 of 37

Data

1. Moskalenko: Specimens were collected by the authors of the article in 2012 and 2013 field seasons on Lake Baikal and preserved in 90% ethanol. A subset of species for sequencing was chosen to represent major lineages within Baikal amphipods, including lineages close to basal nodes. Non-Baikal freshwater palearctic Gammarus species morphologically similar to Baikal amphipods were also considered. Sequences that branched out beyond Gammarus species, namely Hyalella azteca and Dikerogammarus villosus were used as outgroups.

2. Bystřický: Individual lines of G. fossarum in their contact zone in the Western Carpathians are being studied.

There is no source data. But there is an archive already with ready alignments in supporting information: https://aslopubs.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1002%2Flno.12239&file=lno12239-sup-0007-DataS1.zip

Reproducibly extracting data from NCBI with the help of Entrez Utilities (E-Utils) is possible:

​

```

from Bio import Entrez

from Bio import SeqIO

​

Entrez.email = 'Your crustacean email@yandex.ru'

accession_numbers = ["MN005083", "MN005110", "MN005067", "MN148361", "MT110188", ....."XM_018164572", "XM_018168833"] (more details in our github)

​

query = " OR ".join(accession_numbers)

handle = Entrez.efetch(db="nucleotide", id=query, rettype="fasta", retmode="text")

sequences = handle.read()

```

​

However, transcriptome assemblies were sometimes mixed among the sequences, which had to be worked with using blastn.

Moskalenko, V. N.; Neretina, T. V.; Yampolsky, L. Y. To the origin of Lake Baikal endemic gammarid radiations, with description of two new Eulimnogammarus spp. Zootaxa 4766(3), 457-471 (2020). https://doi.org/10.11646/ZOOTAXA.4766.3.5

Bystřický, P.K., Rutová, T., Brož, V., Gajdošová, M., Juračka, P.J., Copilaş-Ciocianu, D. and Petrusek, A. Distribution patterns at different spatial scales reveal reproductive isolation and frequent syntopy among divergent lineages of an amphipod species complex in Western Carpathian streams. Limnol Oceanogr 67, 2796-2808 (2022). https://doi.org/10.1002/lno.12239

​

5 of 37

Data

3. Bukin: Specimens of the amphipod G. fasciatus were collected in the littoral zone of Lake Baikal + previously published in pop-set data.

URL-link for pop-set data: https://www.ncbi.nlm.nih.gov/popset/376403838

The sequences obtained in the study were deposited in the database under the following accession numbers: MG214957, MG214958.

MG214957: https://www.ncbi.nlm.nih.gov/nuccore/MG214957.1/

MG214958: https://www.ncbi.nlm.nih.gov/nuccore/MG214958.1/

4. Wattier: COI sequences for the DNA-barcode region were either newly generated or derived from the literature for 4926 Gammarus fossarum individuals, from 498 sampling sites distributed throughout 19 European countries.

All these sequences were derived into three groups.

  1. All newly produced COI sequences are released in GenBank under accession numbers: MT978656–MT980726. https://www.ncbi.nlm.nih.gov/nuccore/?term=MT978656:MT980726[pacc]
  2. The 463 sequences from Lagrue et al. which were initially unreleased in Genbank are now available under accession numbers: MT411018-MT411480. https://www.ncbi.nlm.nih.gov/nuccore/?term=MT411018:MT411480[pacc]
  3. Sequences and metadata are available in BOLD dataset DS-GFOSCDEU. http://www.boldsystems.org/index.php/Public_SearchTerms?query=DS-GFOSCDEU

Bukin, Y.S., Petunina, J.V. & Sherbakov, D.Y. The Mechanisms for Genetic Diversity of Baikal Endemic Amphipod Gmelinoides fasciatus: Relationships between the Population Processes and Paleoclimatic History of the Lake. Russ J Genet 54, 1059–1068 (2018). https://doi.org/10.1134/S1022795418090053

Wattier, R., Mamos, T., Copilaş-Ciocianu, D. et al. Continental-scale patterns of hyper-cryptic diversity within the freshwater model taxon Gammarus fossarum (Crustacea, Amphipoda). Sci Rep 10, 16536 (2020). https://doi.org/10.1038/s41598-020-73739-0

6 of 37

Methods

Data preparation:

Alignment + trimming:

7 of 37

Methods

Data concatenation and partitions:

Best-fitting evolutionary models:

Building tree + dating:

​

8 of 37

Results

Three checklists have been prepared:

In a better quality they can be found in:

9 of 37

Results

Checklists were also prepared for the report format for github:

In a better quality they can be found in:

10 of 37

Results

The article describes two new species of Baikal gammarids, Eulimnogammarus trengove sp. nov. and Eulimnogammarus tchernykhi sp. nov., with some morphology indicating a basal position within the radiation, but strangely, with further concatenation of the data, we have not to consider them due to the lack of data provided in the article. In addition, the description of datasets that was presented in the methods and materials section does not contain these two types of Baikal gammarids. In addition, a rather interesting fact is that by the type of phylogenetic trees in the article, it seems that the data was still divided into two datasets in a slightly different way (this can be judged due to the difference and the number of species under consideration).

Also to process the data from both the source and the tools obtained during the application, a few small scripts were written to generate the correct input data for the next tools in the queue of our pipeline tools. All of them are presented in the file Datasets.ipynb file. For convenience, they are designed in the form of functions.

Moskalenko

Preparing data for alignment

11 of 37

Results

According to one of the articles, a comparative analysis of the use of tools at different stages of our analysis was carried out.

6 tools for multiple alignment were used, as a result, the tool (MUSCLE) used by the authors of the article did not show particularly great benefits in its use:

Moskalenko

Multiple alignment

Our analysis Multiple_al.ipynb showed that kalign tool wins on all indicators based on our data.

12 of 37

Results

Further, at the stage of data concatenation (Data_concatenation.ipynb), it turned out that a tool like catfasta2phyml is great for getting the output of a data splitting scheme in a file.cfg for PartitionFinder, however, concatenated our data incorrectly, which did not even allow us to trim sequences using the trimAl tool. The MEGA 11 program was chosen for concatenation.

Moskalenko

Data concatenation

Trimming

Dataset 1:

Dataset 2:

All genes without gaps in the data:

​

​

Although this was not stated in the article, but by the look of our alignments, it was decided not to skip the trimming stage. This step of our pipeline is presented in the Trimming.ipynb file.

13 of 37

Results

In this article, this stage is skipped and it is simply written in the methods and materials section that an ML-tree is being built in MEGA7.

Our selection of the evolutionary model was performed using both ModelTest (ModelTest-NG) and ModelFinder in IQ-TREE (Evolution_model.ipynb):

​

In general, we believe that both tools have come to about the same thing. However, in the article, again, this point was omitted and it is not even known which of the tools was used to select the model, or at least which model of evolution they selected and by what criterion (BIC or AIC).

Moskalenko

The selection of the evolutionary model

14 of 37

Results

Apart from the fact that the trees were built in MEGA 7, nothing else was said in the article.

Our ML-trees were built in RAxML-NG and in IQ-TREE ML_trees.ipynb, but the optimal solution seemed to be to use IQ-TREE, since it worked faster.

​

Moskalenko

Building ML-trees

15 of 37

Results

Moskalenko

Dataset 1:

In R (ggtree):

In Python (Phylo):

Building ML-trees

16 of 37

Results

Moskalenko

Building ML-trees

Dataset 2:

In R (ggtree):

In Python (Phylo):

17 of 37

Results

Moskalenko

Building ML-trees

Dataset with all sequences without gaps in the data:

In R (ggtree):

In Python (Phylo):

18 of 37

Results

Moskalenko

Building BA-trees

Dataset 1:

Phylogenies were also reconstructed by Bayesian analysis (BA) in BEAST2:

​

19 of 37

Results

Moskalenko

Building BA-trees

But protein coding sequences (H3) entered as amino acid sequences:

20 of 37

Results

Moskalenko

Conclusions

It is very strange that in the original article, the same authors obtained different results on phylogenetic analysis without adding new data, however, none of their articles indicate either evolutionary models for building trees, or parameters set in some tools.

The main conclusions we have received:

  1. The article described two new species of Baikal gammarids, Eulimnogammarus etingovae sp. nov. and Eulimnogammarus tchernykhi sp. nov., and using the construction of phylogenetic trees, they tried to determine the phylogeny placement of these two new species. However, according to the table and description of datasets provided by the authors of this article, they throw them out of consideration, which is indicated by the type of trees we have built.
  2. The hypothesis that Hyalella azteca is an Outgroup has been confirmed by ML methodology.
  3. According to our analysis, such a unique result was not obtained, as the only species that is consistently attributed to the sister group in relation to the rest of the second Baikal radiation was Pachyschesis branchialis.
  4. The hypothesis was confirmed: The inferred relationship of the two Baikal cades with cosmopolitan palearctic Gammarus species strongly depends on the species and genes included and phylogenetic methodology used.
  5. Finally, ML analysis suggests that Baikal’s first clade isn't sister to the basal palearctic species which contradicts the results obtained in the article.
  6. The lack of agreement between different phylogenetic methods with the low resolution allowed by limited sequence data may be interpreted as an evidence of fast ancient divergence. The hypothesis is confirmed.

​

21 of 37

Results

Bystřický

  • One of the main problems of this article is the lack of source data Dataset.ipynb, we found only alignment files in the appendix to the article (supporting information).
  • Without having the source data, it is quite difficult to judge how the alignment files were obtained, besides, in the future we will face the difficulty that the files are not concatenated due to the fact that they contain different sequences and in different order.
  • It is impossible to reproduce these results, the trimming stage was also not described in the article, so in order to avoid any errors, it was not performed.
  • This cannot be reproduced using the provided alignment files. In addition, according to the type of phylogenetic trees, they may have been analyzed separately.

Preparing data for alignment

Multiple alignment

Trimming

Data concatenation

22 of 37

Results

Bystřický

We also refuse to use the PartitionFinder tool, since it is quite difficult to install and in addition, it was not even indicated which model was chosen and by what criteria. Therefore, we chose an evolutionary model using two other tools (RAxML-NG and IQ-TREE) by BIC (Ev_model.ipynb):

Apart from the last file, the tools predicted similar models.

The selection of the evolutionary model

23 of 37

Results

Bystřický

Next, trees were built using RAxML-NG and in IQ-TREE ML_trees_b.ipynb, and using the python module from Bio import Phylo, they were visualized.

16S_alignment.fas:

IQ-TREE

RAxML-NG

Their result

Building ML-trees

24 of 37

Results

Bystřický

28S_alignment.fas:

IQ-TREE

RAxML-NG

Their result

Building ML-trees

25 of 37

Results

Bystřický

EF1a_alignment_phased.fas:

IQ-TREE

RAxML-NG

Their result

Building ML-trees

26 of 37

Bystřický

Results

For 28S, which is homozygous and lineages exhibit reciprocal monophyly (except CWE D), we applied both a tree-based method.

27 of 37

Bystřický

Results

The EF1a marker is heterozygous and lineages are not reciprocally monophyletic. Therefore, tree- or distance-based species delimitation methods are not applicable. Consequently, we used the haploweb method since it is specifically designed to delimit species based on such markers . The haploweb distinguishes putative species clusters based on mutual allelic exclusivity, following Doyle’s view of species as fields for recombination. It was run on the HaplowebMaker web-server (https://eeg-ebe.github.io/HaplowebMaker/):

Species delimitation based on the haploweb method applied to the phased EF1a alleles. Each circle represents a unique haplotype. Haplotypes co-occurring within the same individual are connected by a colored curved line. Each tick mark represents one substitution:

28 of 37

Bystřický

Results

  1. Maximum likelihood family trees for each analyzed marker.
  2. The same synthetic individuals belonging to different pedigrees were obtained,
  3. The outer group is everywhere Gammarus pulex.
  4. Note that both alleles of the EF1a marker are present in the tree.
  5. The haploweb method, applied on EF1a data, confirmed the lack of hybridization since each lineage harbored unique alleles.

Conclusions

29 of 37

Results

Bukin

There was no information about trimming in the original paper

It was decided to make one with -automated1 setting

30 of 37

Results

Bukin

In the original article the best-fitting GTR model with gamma correction

(GTR + G) was used for phylogenetic analysis.

During the replication was got another model

31 of 37

Results

Bukin

On the our picture we can see the original clades of Central, Northern and southwestern populations

Little difference may be the result of poor given data

32 of 37

Results

Wattier

There was no information about trimming in the original paper

But there was problem with turning on pipeline due to the sample size. So we had to trim four times with -automated1 setting

33 of 37

Results

Wattier

The best-fitting model of substitution was the transversion model (TVM) with gamma-distributed rate heterogeneity (G) and a given proportion of invariable sites (I)

Received model matched with the original

34 of 37

Results

Wattier

It isn’t possible to compare trees, due to the Jawa heap space error, which hasn’t been solved yet

35 of 37

Summary

Two of the four articles should be considered complex and poorly reproducible for the following reasons:

  1. Lack of code availability: Without access to the original code, it becomes difficult for other researchers to understand and replicate the analysis;
  2. The absence of source data in our phylogenetic analysis studies was a significant barrier to reproducibility. Without access to the original source data, other researchers like us cannot replicate the analysis or verify the findings;
  3. Missing or incomplete data: Sometimes, articles do not provide the full dataset used in the analysis, making it impossible to reproduce the results accurately. Also, information about the source data and explanations to them are sometimes located throughout the article, which is why you have to go back through the previous steps of the pipeline;
  4. Incomplete or insufficient method descriptions: Some articles may not provide detailed descriptions of the methods used in the phylogenetic analysis. This lack of information makes it challenging for other researchers to accurately reproduce the analysis. Insufficient method descriptions can include incomplete information about the software versions, parameter settings, and other critical details necessary for replication;
  5. Ambiguities in data preprocessing: Without clear instructions or descriptions, reproducing the exact preprocessing steps becomes challenging and can lead to inconsistencies in results.The authors miss the moments of correcting output and, accordingly, input data during pipeline execution;

The first item on the list of shortcomings and the last two items also apply to the other two articles.

More details can be found in:

36 of 37

Summary

Our pipeline, compiled by the authors of this project

More details can be found in:

37 of 37

Further plans

​

​

​

​

Reproducibility is the main goal that Conda and other environment management tools strive to achieve in order to ensure high—quality code execution not only on a personal PC, but also on the machines of team colleagues. That's why it's so important to be able to recreate the same environment on another computer. To do this, we can create a unique virtual environment and design the pipeline as a Bash script.