diff --git a/init.R b/init.R
index e5d7235..05c28e4 100644
--- a/init.R
+++ b/init.R
@@ -1,6 +1,7 @@
 source("clustering.R")
 source("correlation.R")
 source("input.R")
+source("neural.R")
 source("util.R")
 
 # Load input data
@@ -56,6 +57,24 @@ correlation_replicative <- run_cached(
     tpe_old_genes
 )
 
+neural_all <- run_cached(
+    "neural_all",
+    process_neural,
+    distances,
+    all_species,
+    all_genes,
+    tpe_old_genes
+)
+
+neural_replicative <- run_cached(
+    "neural_replicative",
+    process_neural,
+    distances,
+    replicative_species,
+    all_genes,
+    tpe_old_genes
+)
+
 # Merge processed data as well as gene information.
 
 results_all <- merge(
@@ -72,6 +91,13 @@ results_all <- merge(
     by.y = "gene"
 )
 
+results_all <- merge(
+    results_all,
+    neural_all,
+    by.x = "id",
+    by.y = "gene"
+)
+
 results_replicative <- merge(
     genes,
     clustering_replicative,
@@ -86,6 +112,13 @@ results_replicative <- merge(
     by.y = "gene"
 )
 
+results_replicative <- merge(
+    results_replicative,
+    neural_replicative,
+    by.x = "id",
+    by.y = "gene"
+)
+
 # Rename `id` columns to `gene`.
 
 setnames(results_all, "id", "gene")
diff --git a/neural.R b/neural.R
new file mode 100644
index 0000000..1a1611a
--- /dev/null
+++ b/neural.R
@@ -0,0 +1,110 @@
+library(data.table)
+library(neuralnet)
+
+#' Find genes by training a neural network on reference position data.
+#'
+#' The result will be a data.table with the following columns:
+#'
+#'  - `gene` Gene ID of the processed gene.
+#'  - `neural` Output score given by the neural network.
+#'
+#' @param distances Distance data to use.
+#' @param species_ids Species, whose data should be included.
+#' @param gene_ids Genes to process. This should include the reference genes.
+#' @param reference_gene_ids Genes to compare to.
+#' @param seed A seed to get reproducible results.
+process_neural <- function(distances, species_ids, gene_ids,
+                           reference_gene_ids, seed = 726839) {
+    set.seed(seed)
+    gene_count <- length(gene_ids)
+
+    # Prefilter distances by species.
+    distances <- distances[species %chin% species_ids]
+
+    #' Input data for the network. This contains the gene ID as an identifier
+    #' as well as the per-species gene distances as input variables.
+    data <- data.table(gene = gene_ids)
+
+    #' Buffer to keep track of species included in the computation. Species
+    #' from `species_ids` may be excluded if they don't have enough data.
+    species_ids_included <- NULL
+
+    for (species_id in species_ids) {
+        # Make a column specific to this species.
+
+        species_distances <- distances[species == species_id, .(gene, distance)]
+        setnames(species_distances, "distance", species_id)
+
+        # Only include species with at least 25% known values.
+
+        species_distances <- na.omit(species_distances)
+
+        if (nrow(species_distances) >= 0.25 * gene_count) {
+            species_ids_included <- append(species_ids_included, species_id)
+            data <- merge(data, species_distances, all = TRUE)
+        }
+    }
+
+    # Replace missing data with mean values. The neural network can't handle
+    # NAs in a meaningful way. Choosing extreme values here would result in
+    # heavily biased results. Therefore, the mean value is chosen as a
+    # compromise. However, this will of course lessen the significance of the
+    # results.
+    for (species_id in species_ids_included) {
+        mean_value <- data[, mean(get(species_id), na.rm = TRUE)]
+
+        data[
+            is.na(get(species_id)),
+            eval(quote(species_id)) := round(mean_value)
+        ]
+    }
+
+    # Extract the reference genes.
+
+    reference_data <- data[gene %chin% reference_gene_ids]
+    reference_data[, neural := 1.0]
+
+    # Take out random samples from the remaining genes. This is another
+    # compromise with a negative impact on significance. Because there is no
+    # information on genes with are explicitely *not* TPE-OLD genes, we have to
+    # assume that a random sample of genes has a low probability of including
+    # TPE-OLD genes.
+
+    without_reference_data <- data[!gene %chin% reference_gene_ids]
+
+    reference_samples <- without_reference_data[
+        sample(
+            nrow(without_reference_data),
+            nrow(reference_data)
+        )
+    ]
+
+    reference_samples[, neural := 0.0]
+
+    # Merge training data. The training data includes all reference genes as
+    # well as an equal number of random sample genes.
+    training_data <- rbindlist(list(reference_data, reference_samples))
+
+    # Construct and train the neural network.
+
+    nn_formula <- as.formula(sprintf(
+        "neural~%s",
+        paste(species_ids_included, collapse = "+")
+    ))
+
+    layer1 <- length(species_ids_included) * 0.66
+    layer2 <- layer1 * 0.66
+    layer3 <- layer2 * 0.66
+
+    nn <- neuralnet(
+        nn_formula,
+        training_data,
+        hidden = c(layer1, layer2, layer3),
+        linear.output = FALSE
+    )
+
+    # Return the resulting scores given by applying the neural network.
+
+    data[, neural := compute(nn, data)$net.result]
+    data[, .(gene, neural)]
+}
\ No newline at end of file
diff --git a/server.R b/server.R
index deff16b..a1fce15 100644
--- a/server.R
+++ b/server.R
@@ -30,12 +30,14 @@ server <- function(input, output) {
 
         clusteriness_weight <- input$clusteriness / 100
         correlation_weight <- input$correlation / 100
-        total_weight <- clusteriness_weight + correlation_weight
+        neural_weight <- input$neural / 100
+        total_weight <- clusteriness_weight + correlation_weight + neural_weight
         clusteriness_factor <- clusteriness_weight / total_weight
         correlation_factor <- correlation_weight / total_weight
+        neural_factor <- neural_weight / total_weight
 
         results[, score := clusteriness_factor * clusteriness +
-            correlation_factor * r_mean]
+            correlation_factor * r_mean + neural_factor * neural]
 
         # Apply the cut-off score.
 
@@ -55,6 +57,7 @@ server <- function(input, output) {
                 name,
                 clusteriness,
                 r_mean,
+                neural,
                 score
             )],
             rownames = FALSE,
@@ -64,6 +67,7 @@ server <- function(input, output) {
                 "",
                 "Clusters",
                 "Correlation",
+                "Neural",
                 "Score"
             ),
             style = "bootstrap",
@@ -73,7 +77,11 @@ server <- function(input, output) {
             )
         )
 
-        formatPercentage(dt, c("clusteriness", "r_mean", "score"), digits = 1)
+        formatPercentage(
+            dt,
+            c("clusteriness", "r_mean", "neural", "score"),
+            digits = 1
+        )
     })
 
     output$synposis <- renderText({
diff --git a/ui.R b/ui.R
index 05aec8f..e41886f 100644
--- a/ui.R
+++ b/ui.R
@@ -36,6 +36,15 @@ ui <- fluidPage(
                 step = 1,
                 value = 100
             ),
+            sliderInput(
+                "neural",
+                "Assessment by neural network",
+                post = "%",
+                min = 0,
+                max = 100,
+                step = 1,
+                value = 100
+            ),
             sliderInput(
                 "cutoff",
                 "Cut-off score",