comparison clean_ltr.R @ 4:93d35ae65e1b draft

"planemo upload commit 57a4f4a749b60b4e1d992dc3a879add7bb4bb56b"

3:6ae4a341d1f3

cat("reading gff...")
g <- rtracklayer::import(opt$gff3, format = 'gff3')  # DANTE gff3
cat("done\n")
# testing
if (FALSE) {
  g <- rtracklayer::import("./test_data/sample_ltr_annotation.gff3")
  s <- readDNAStringSet("./test_data/sample_genome.fasta")

  g <- rtracklayer::import("./test_data/DANTE_LTR_Vfaba_chr5.gff3")
  s <- readDNAStringSet("./test_data/211010_Vfaba_chr5.fasta")

  lineage_info <- read.table("databases/lineage_domain_order.csv", sep = "\t", header =
    TRUE, as.is = TRUE)
  source("./R/ltr_utils.R")
}


# get te sequence based on rank

# evaluate best TE -  DLTP grou
s_te <- get_te_sequences(g, s)  # split by 'element quality'

TE_DLTP_info <- analyze_TE(s_te$DLTP, word_size = word_size, ncpus = ncpus)

# TE rank 2:
TE_DLT_plus_DLP_info <- analyze_TE(c(s_te$DLP, s_te$DLT), word_size = word_size, ncpus
  = ncpus)
TE_DLT_plus_DLP_info_DLTP_verified <- compare_TE_datasets(c(s_te$DLT, s_te$DLP), ncpus
  = ncpus,
                                                          TE_DLTP_info$seqs_representative,
                                                          word_size = word_size
)

TE_DLT_plus_DLP_info_multiplicity <- verify_based_on_multiplicity(TE_DLT_plus_DLP_info)
# create additional library from rank 2 verified by multiplicity
id_for_additional_library <- setdiff(
  TE_DLT_plus_DLP_info_multiplicity$id_ok_mp_verified,

  }else {
    seq_representative <- TE_DLTP_info$seqs_representative
  }
}

# TE  rank 3
TE_DL_info_DLTP_verified <- compare_TE_datasets(
  s_te$DL,
  TE_DLTP_info$seqs_representative, min_coverage = 0.98,
  ncpus = ncpus
)


R <- seq_diversity(seq_representative)$richness
SI <- seq_diversity(seq_representative)$shannon_index

4:93d35ae65e1b

cat("reading gff...")
g <- rtracklayer::import(opt$gff3, format = 'gff3')  # DANTE gff3
cat("done\n")
# testing
if (FALSE) {
  g <- rtracklayer::import("./test_data/big_test_data/dante_ltr_unfiltered_t.cacao.gff3")
  s <- readDNAStringSet("./test_data/big_test_data/T_cacao_chromosomes.fasta")

  g <- rtracklayer::import("./test_data/sample_ltr_annotation.gff3")
  s <- readDNAStringSet("./test_data/sample_genome.fasta")

  g <- rtracklayer::import("./test_data/DANTE_LTR_Vfaba_chr5.gff3")
  s <- readDNAStringSet("./test_data/211010_Vfaba_chr5.fasta")

  lineage_info <- read.table("databases/lineage_domain_order.csv", sep = "\t", header =
    TRUE, as.is = TRUE)
  source("./R/ltr_utils.R")
}

## ID in g must be unique - this could be a problem if gff is concatenated from multiple files!
## id ID is renamed - rename parent to!
## add chromosom index to disctinguish same IDs
suffix <- as.numeric(seqnames(g))
g$ID <- ifelse(is.na(g$ID), NA, paste0(g$ID,"_", suffix))
g$Parent <- ifelse(is.na(g$Parent), NA, paste0(g$Parent,"_", suffix))

# get te sequence based on rank

# evaluate best TE -  DLTP grou
s_te <- get_te_sequences(g, s)  # split by 'element quality'

TE_DLTP_info <- analyze_TE(s_te$DLTP, word_size = word_size, ncpus = ncpus)

# TE rank 2:
TE_DLT_plus_DLP_info <- analyze_TE(c(s_te$DLP, s_te$DLT), word_size = word_size, ncpus
  = ncpus)
  TE_DLT_plus_DLP_info_DLTP_verified <- compare_TE_datasets(c(s_te$DLT, s_te$DLP), ncpus
    = ncpus,TE_DLTP_info$seqs_representative, word_size = word_size
  )

TE_DLT_plus_DLP_info_multiplicity <- verify_based_on_multiplicity(TE_DLT_plus_DLP_info)
# create additional library from rank 2 verified by multiplicity
id_for_additional_library <- setdiff(
  TE_DLT_plus_DLP_info_multiplicity$id_ok_mp_verified,

  }else {
    seq_representative <- TE_DLTP_info$seqs_representative
  }
}

# TE  rank 3 - verify agains good DLTP
TE_DL_info_DLTP_verified <- compare_TE_datasets(
  s_te$DL,
  TE_DLTP_info$seqs_representative, min_coverage = 0.95,
  ncpus = ncpus, word_size = word_size
)


R <- seq_diversity(seq_representative)$richness
SI <- seq_diversity(seq_representative)$shannon_index