Mercurial > repos > petr-novak > dante_ltr

annotate clean_ltr.R @ 5:0c3111ab729b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit 5a81f1734493a3acf5ff24a573fa2887d92b58a8"
author petr-novak
date Mon, 16 May 2022 07:50:41 +0000
parents 93d35ae65e1b
children b91ca438a1cb
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
1 #!/usr/bin/env Rscript
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
2 initial_options <- commandArgs(trailingOnly = FALSE)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
3 file_arg_name <- "--file="
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
4 script_name <- normalizePath(sub(file_arg_name, "", initial_options[grep
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
5 (file_arg_name,
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
6 initial_options)]))
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
7 script_dir <- dirname(script_name)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
8 library(optparse)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
9
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
10 parser <- OptionParser()
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
11 option_list <- list(
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
12 make_option(c("-g", "--gff3"), action = "store", type = "character",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
13 help = "gff3 with LTR Transposable elements", default = NULL),
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
14 make_option(c("-s", "--reference_sequence"), action = "store", type = "character",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
15 help = "reference sequence as fasta",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
16 default = NULL),
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
17 make_option(c("-o", "--output"), action = "store", type = "character",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
18 help = "output file prefix", default = NULL),
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
19 make_option(c("-c", "--cpu"), type =
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
20 "integer", default = 5, help = "Number of cpu to use [default %default]",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
21 metavar = "number")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
22
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
23 )
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
24 description <- paste(strwrap(""))
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
25
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
26 epilogue <- ""
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
27 parser <- OptionParser(option_list = option_list, epilogue = epilogue, description =
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
28 description, usage = "usage: %prog COMMAND [OPTIONS]")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
29 opt <- parse_args(parser, args = commandArgs(TRUE))
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
30
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
31
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
32 # load packages
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
33 suppressPackageStartupMessages({ library(rtracklayer)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
34 library(Biostrings)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
35 library(BSgenome)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
36 library(parallel)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
37 })
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
38
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
39 # CONFIGURATION
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
40 # load configuration files and functions:
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
41 lineage_file <- paste0(script_dir, "/databases/lineage_domain_order.csv")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
42 ltr_utils_r <- paste0(script_dir, "/R/ltr_utils.R")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
43 if (file.exists(lineage_file)) {
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
44 lineage_info <- read.table(lineage_file, sep = "\t",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
45 header = TRUE,
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
46 as.is = TRUE)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
47 source(ltr_utils_r)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
48 }else {
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
49 lineage_file <- paste0(script_dir, "/.
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
50 ./share/dante_ltr/databases/lineage_domain_order.csv")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
51 ltr_utils_r <- paste0(script_dir, "/.
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
52 ./share/dante_ltr/R/ltr_utils.R")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
53 if (file.exists(lineage_file)) {
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
54 lineage_info <- read.table(lineage_file, sep = "\t",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
55 header = TRUE,
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
56 as.is = TRUE)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
57 source(ltr_utils_r)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
58 }else {
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
59 (stop("configuration files not found"))
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
60 }
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
61 }
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
62
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
63
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
64 ncpus <- opt$cpu
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
65
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
66
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
67 # load data
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
68 cat("reading fasta...")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
69 s <- readDNAStringSet(opt$reference_sequence) # genome assembly
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
70 cat("done\n")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
71 outfile <- opt$output
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
72 # clean sequence names:
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
73 names(s) <- gsub(" .+", "", names(s))
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
74 cat("reading gff...")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
75 g <- rtracklayer::import(opt$gff3, format = 'gff3') # DANTE gff3
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
76 cat("done\n")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
77 # testing
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
78 if (FALSE) {
4
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
79 g <- rtracklayer::import("./test_data/big_test_data/dante_ltr_unfiltered_t.cacao.gff3")
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
80 s <- readDNAStringSet("./test_data/big_test_data/T_cacao_chromosomes.fasta")
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
81
0
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
82 g <- rtracklayer::import("./test_data/sample_ltr_annotation.gff3")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
83 s <- readDNAStringSet("./test_data/sample_genome.fasta")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
84
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
85 g <- rtracklayer::import("./test_data/DANTE_LTR_Vfaba_chr5.gff3")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
86 s <- readDNAStringSet("./test_data/211010_Vfaba_chr5.fasta")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
87 names(s) <- gsub(" .+", "", names(s))
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
88 ncpus <- 10
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
89 lineage_info <- read.table("databases/lineage_domain_order.csv", sep = "\t", header =
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
90 TRUE, as.is = TRUE)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
91 source("./R/ltr_utils.R")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
92 }
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
93
4
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
94 ## ID in g must be unique - this could be a problem if gff is concatenated from multiple files!
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
95 ## id ID is renamed - rename parent to!
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
96 ## add chromosom index to disctinguish same IDs
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
97 suffix <- as.numeric(seqnames(g))
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
98 g$ID <- ifelse(is.na(g$ID), NA, paste0(g$ID,"_", suffix))
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
99 g$Parent <- ifelse(is.na(g$Parent), NA, paste0(g$Parent,"_", suffix))
0
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
100
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
101 # get te sequence based on rank
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
102
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
103 # evaluate best TE - DLTP grou
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
104 s_te <- get_te_sequences(g, s) # split by 'element quality'
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
105 # best quality - split by lineage
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
106 word_size <- 15
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
107 # best TE
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
108 TE_DLTP_info <- analyze_TE(s_te$DLTP, word_size = word_size, ncpus = ncpus)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
109
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
110 # TE rank 2:
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
111 TE_DLT_plus_DLP_info <- analyze_TE(c(s_te$DLP, s_te$DLT), word_size = word_size, ncpus
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
112 = ncpus)
4
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
113 TE_DLT_plus_DLP_info_DLTP_verified <- compare_TE_datasets(c(s_te$DLT, s_te$DLP), ncpus
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
114 = ncpus,TE_DLTP_info$seqs_representative, word_size = word_size
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
115 )
0
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
116
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
117 TE_DLT_plus_DLP_info_multiplicity <- verify_based_on_multiplicity(TE_DLT_plus_DLP_info)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
118 # create additional library from rank 2 verified by multiplicity
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
119 id_for_additional_library <- setdiff(
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
120 TE_DLT_plus_DLP_info_multiplicity$id_ok_mp_verified,
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
121 TE_DLT_plus_DLP_info_DLTP_verified$id_ok_verified)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
122
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
123 if (length(id_for_additional_library) > 1) {
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
124 seqs_for_additional_library <- c(s_te$DLP, s_te$DLT)[names(c(s_te$DLP, s_te$DLT)) %in%
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
125 id_for_additional_library]
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
126 seqs_additional_info <- analyze_TE(seqs_for_additional_library, word_size =
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
127 word_size, ncpus = ncpus)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
128 seq_representative <- c(
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
129 TE_DLTP_info$seqs_representative,
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
130 seqs_additional_info$seqs_representative
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
131 )
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
132 }else {
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
133 if (length(id_for_additional_library) == 1) {
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
134 seq_representative <- c(
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
135 TE_DLTP_info$seqs_representative,
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
136 c(s_te$DLP, s_te$DLT)[names(c(s_te$DLP, s_te$DLT)) %in% id_for_additional_library]
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
137 )
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
138 }else {
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
139 seq_representative <- TE_DLTP_info$seqs_representative
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
140 }
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
141 }
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
142
4
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
143 # TE rank 3 - verify agains good DLTP
0
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
144 TE_DL_info_DLTP_verified <- compare_TE_datasets(
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
145 s_te$DL,
4
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
146 TE_DLTP_info$seqs_representative, min_coverage = 0.95,
93d35ae65e1b "planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"
petr-novak
parents: 3
diff changeset
147 ncpus = ncpus, word_size = word_size
0
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
148 )
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
149
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
150
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
151 R <- seq_diversity(seq_representative)$richness
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
152 SI <- seq_diversity(seq_representative)$shannon_index
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
153
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
154 # final RM library:
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
155 seq_representative_no_ssr <- seq_representative[R > 20]
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
156
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
157 ID <- g$ID[g$type == "transposable_element"]
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
158 names(ID) <- paste0(seqnames(g), "_",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
159 start(g), "_",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
160 end(g), "#",
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
161 g$Final_Classification
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
162 )[g$type %in% "transposable_element"]
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
163
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
164 # create clean gff3
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
165 id_of_good_te <- unique(c(
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
166 TE_DLTP_info$te_conflict_info$ok,
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
167 TE_DLT_plus_DLP_info_DLTP_verified$id_ok_verified,
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
168 TE_DLT_plus_DLP_info_multiplicity$id_ok_mp_verified,
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
169 TE_DL_info_DLTP_verified$id_ok_verified)
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
170 )
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
171
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
172 c1 <- g$ID %in% ID[id_of_good_te]
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
173 c2 <- sapply(g$Parent, function(x)ifelse(length(x) == 0, "", x)) %in% ID[id_of_good_te]
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
174
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
175 gff_out <- g[c1 | c2]
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
176
3
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
177 gff_te <- gff_out[gff_out$type %in% "transposable_element"]
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
178 gff_5ltr <- gff_out[gff_out$LTR %in% "5LTR"]
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
179 gff_3ltr <- gff_out[gff_out$LTR %in% "3LTR"]
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
180 full_te <- getSeqNamed(s, gff_te)
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
181 ltr5 <- getSeqNamed(s, gff_5ltr)
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
182 ltr3 <- getSeqNamed(s, gff_3ltr)
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
183 inc <- gff_te$Rank != "DL"
0
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
184
3
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
185 writeXStringSet(seq_representative, paste0(opt$output, "_RM_lib_non_redundant.fasta"))
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
186 writeXStringSet(full_te, paste0(opt$output, "_RM_lib_full_TE.fasta"))
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
187 writeXStringSet(ltr5, paste0(opt$output, "_RM_lib_5LTR.fasta"))
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
188 writeXStringSet(ltr3, paste0(opt$output, "_RM_lib_3LTR.fasta"))
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
189
0
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
190 export(gff_out, paste0(opt$output, "_clean.gff3"), format = "gff3")
7b0bbe7477c4 "planemo upload commit 92c684dff3b377c8c08654c7f3d46a133385e3e0-dirty"
petr-novak
parents:
diff changeset
191
3
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
192 lv <- sort(unique(gff_te$Final_Classification))
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
193 te_count <- table(factor(gff_te$Final_Classification, levels=lv))
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
194
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
195 pdf(paste0(opt$output, "_summary.pdf"), width = 13, height=8, pointsize = 10)
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
196 par(mfrow=c(1,2), mar=c(5,7,2,0), las=1)
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
197 boxplot(width(gff_te) ~ factor(gff_te$Final_Classification, levels=lv),
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
198 horizontal = TRUE, xlab="length[bp]", ylab="",
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
199 names = paste0(gsub("^.+[|]", "", lv), " (", te_count, ")"),
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
200 main = 'Full TE', at = seq_along(lv)*4
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
201 )
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
202 boxplot(width(gff_te[inc]) ~ factor(gff_te$Final_Classification[inc], levels=lv),
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
203 horizontal = TRUE, xlab="length[bp]", ylab="",
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
204 names = rep("", length(lv)),
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
205 main = 'Full TE', at = seq_along(lv)*4-1, add=TRUE, col=2
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
206 )
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
207 par(mar=c(5,0,2,7))
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
208 boxplot(width(gff_5ltr) ~ factor(gff_5ltr$Final_Classification, levels=lv),
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
209 horizontal = TRUE, xlab="length[bp]", ylab="",
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
210 names = rep("", length(lv)),
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
211 main = "5'LTR", at = seq_along(lv)*4
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
212 )
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
213 boxplot(width(gff_5ltr[inc]) ~ factor(gff_5ltr$Final_Classification[inc], levels=lv),
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
214 horizontal = TRUE, xlab="length[bp]", ylab="",
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
215 names = rep("", length(lv)),
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
216 main = "5'LTR", at = seq_along(lv)*4-1, add=TRUE, col=2
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
217 )
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
218 legend('bottomright', col=c("grey","2"), legend=c("All TE", "TE with PBS/TSD"), pch=15)
6ae4a341d1f3 "planemo upload commit 8bd6029a4de4a8f5031a5cc71303bb06217cc88a"
petr-novak
parents: 0
diff changeset
219 dev.off()