diff consensus_from_alignments.R @ 0:0ccbe1c20fc3 draft default tip

planemo upload for repository https://github.com/ColineRoyaux/Galaxy_tool_projects/tree/main/consensus_from_alignments commit ecc21de8f368c6a95c57d4e6511ed42af9e72a66
author ecology
date Tue, 25 Apr 2023 10:05:29 +0000
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/consensus_from_alignments.R	Tue Apr 25 10:05:29 2023 +0000
@@ -0,0 +1,43 @@
+#Rscript
+
+################################################################################
+##      Extract consensus sequence from aligned forward and reverse fasta     ##
+################################################################################
+
+#####Packages
+library(bioseq, quietly = TRUE)
+
+##Load arguments
+args <- commandArgs(trailingOnly = TRUE)
+
+if (length(args) == 0) {
+    stop("This tool needs at least one argument")
+} else {
+    fasta_f <- args[1]
+    seq_type <- args[2]
+    meth_choice <- args[3]
+    gaps_tf <- as.logical(args[4])
+    out_og <- as.logical(args[5])
+}
+
+## Read input file
+seq_l <- bioseq::read_fasta(fasta_f, type = seq_type)
+
+if(bioseq::seq_nseq(seq_l) < 2){
+    stop("Only one sequence in the file, at least two aligned sequences are needed to compute a consensus")
+}else{
+    if(length(unique(bioseq::seq_nchar(seq_l))) > 1) {stop("Sequences have different lengths, please provide aligned sequences")}
+}
+
+##Consensus sequence
+seq_con <- bioseq::seq_consensus(seq_l, method = meth_choice, gaps = gaps_tf)
+
+if(bioseq::seq_nseq(seq_con) > 1){stop("Consensus hasn't worked for an unknown reason, double-check your input file and the parameters you chose")}
+
+names(seq_con) <- paste0("consensus_", Reduce(PTXQC::LCS, names(seq_l)))
+##Create output
+if(out_og){
+    bioseq::write_fasta(c(seq_con, seq_l), file = "output.fasta", line_length = Inf, block_length = Inf)
+}else{
+    bioseq::write_fasta(seq_con, file = "output.fasta", line_length = Inf, block_length = Inf)
+}