Introduction to the iimi package

1.1. Installation

First, let’s install necessary packages. You may skip this step if you have installed the packages before.

# install iimi
install.packages(c("iimi", "httr"))

# install Biostrings
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("Biostrings")

1.2. Loading packages

We will load necessary packages before we start any analysis.

library(iimi)

library(Biostrings)

library(httr)

1.3. Pre-processing

To get started creating coverage profiles and feature-extracted data frames, we need mapping results. We used Bowtie2 to map the samples (paired-end or single-end) against the official Virtool virus data base (ver. 1.4.0). We recommend Bowtie2 or minimap2 since we have tried both and they yield similar result with minimap2 having a slight decrease. We let both software to report all alignments (-a mode for Bowtie2, --secondary=yes for minimap2). You can also use other mapping tools. Make sure that you convert your mapping results to indexed and sorted BAM files using Samtools.

1.4. Downloading example BAM files

We provide three example BAM files to demonstrate how to use iimi. These files are sourced from the dataset used in the VirHunter paper (Sukhorukov et al. 2022), which we also utilized for external validation in our manuscript (Ning et al. 2024). You can download these files directly from the Recherche Data Gouv website (Candresse, Marais-Colombel, and Brault 2022). We recommend storing all BAM files in a single folder for ease of access. Here, we provide a short tutorial to guide you through using iimi functions to make predictions on the real data.

2.1. State the path to the folder of your BAM files

If you already have coverage profiles in run-length encoding (RLE) format, go to section 2.2.

path_to_bamfiles <- list.files(
  path = "path/to/your/BAM/files/folder",
  pattern = "bam$",
  full.names = TRUE,
  include.dirs = TRUE
)

2.2. Create a data frame that contains the coverage profiles.

example_cov <- convert_bam_to_rle(bam_file = "path_to_bamfiles")

df <-
  convert_rle_to_df(covs = example_cov, unreliable_region_df = unreliable_regions)

3.1 Using pre-trained models and no customization

If you wish to use provided training model, only input your data to newdata and choose a method of your wish using predict_iimi().

There are three methods that you may choose from: xgb, en, and rf, which stand for pre-trained XGBoost, elastic net, and random forest models. The example below uses the pre-trained XGBoost model.

prediction_default <- predict_iimi(newdata = df, method = "xgb")

The detection of your plant sample(s) is finished. The prediction is TRUE if virus infected the sample, FALSE if virus did not infect the sample.

3.2. Customizing your own model

If you would like to train your own model and use this model to test your data, you can use the codes below to train a new model with your own data.

Ideally, the number of the samples used to train the model should be bigger than 100. However, we are only providing a tutorial on how to use the train_iimi() function, only two samples are used to train the model since example_cov() only contains three in-house data’s coverage information.

First, we need to prepare our training data. We are using a 80/20 random split to split the three samples. This means that two samples are used as the training data, and one sample is used as the testing data. If you are training your own data, treat the training data as your data that you want to train the model on; treat the testing data as your data that you would like to have a prediction on.

If you have separate training and testing data, you may only need to convert your data from mapping results (BAM files) to fetaure-extracted data frame and name your data as the following variables.

Here are some definitions/explanation of the objects to input in train_iimi():

1. train_x: the feature-extracted data frame of plant samples that you would like to train iimi model on. Make sure that you have mapped the samples to the virus database and converted the mapping result to sorted and indexes BAM files.

2. train_y: the known truth or labels for your train_x data. Please label the data to make sure that it has a detection label for virus segments as well.

3. test_x: the feature-extracted data frame of plant samples that you would like to predict using your trained iimi model. Make sure that you have mapped the samples to the virus database and converted the mapping result to sorted and indexes BAM files.

# set seed
set.seed(123)

# spliting into 80-20 train and test data set with the three plant samples
train_names <- sample(levels(as.factor(df$sample_id)),
                      length(unique(df$sample_id)) * 0.8)

# trian data
# train_x is the feature-extracted data frame of your train data
train_x = df[df$sample_id %in% train_names,]

# train_y is the known truth or labels for your train_x data, indicating the presence of specific viruses in the samples
train_y = c()

for (ii in 1:nrow(train_x)) {
  train_y = append(train_y, example_diag[train_x$seg_id[ii],
                                         train_x$sample_id[ii]])
}

# test data
# test_x is the feature-extracted data frame of the data you would like to predict
# here we used the sample that is not in the training set for demonstration purpose
test_x = df[df$sample_id %in% train_names == F,]

Then, we plug in the variables into the train_iimi function with the default XGBoost model.

fit <- train_iimi(train_x = train_x, train_y = train_y)

Now, we have a trained model using the toy data.

Then, the process to detect which viruses infect the plant sample(s) is the same as previously described, except we are using a customized trained model.

prediction_customized <-
  predict_iimi(newdata = test_x,
               trained_model = fit)

The detection of the plant sample(s) is finished. The interpretation is the same as above.

3.3. Customizing unreliable regions

Note: if you would like to create your own unreliable regions, please customize them first, then extract features to build a data frame from section 2.2.2. using customized unreliable regions.

If you would like to create your own mappability profile and high nucleotide content regions besides from using your own training model, you may use create_mappability_profile() and high_nucleotide_regions(). Both functions’ output is a data frame with the start and end position of the unmappable region, the virus that the region is on, and the category that it belongs to.

Mappability profile is a profile of areas on a virus genome that can be mapped to (1) other viruses or (2) host genome. We choose Arabidopsis Thaliana as our host genome.

High nucleotide content regions is a profile of areas on a virus genome that has (1) high GC content and/or (2) high A nucleotide percentage.

Including these two profiles into iimi ensures that there are no false peaks like the ones described in the previous section.

# if you would like to keep unmappable regions that can be mapped to other viruses or the host genome separate into two data frames, you may use the following code:

# input your own path that you would want to store regions on a virus that can be mapped to another virus
# you can customize the name of this type of mappability profile
mappability_profile_virus <-
  create_mappability_profile("path/to/bam/files/folder/virus", category = "Unmappable region (virus)")

# input your own path that you would want to store regions on a virus that can be mapped to the host genome
# you can customize the name of this type of mappability profile
mappability_profile_host <-
  create_mappability_profile("path/to/bam/files/folder/host", category = "Unmappable region (host)")

# if you would like to keep everything in the same data frame, you may use the following code:
mappability_profile <-
  create_mappability_profile("path/to/bam/files/folder/of/both/types/", category = "Unmappable region")

high_nucleotide_regions <-
  create_high_nucleotide_content(gc = 0.6, a = 0.45)