Mercurial > repos > padge > trimal
diff trimal_repo/scripts/get_sequences_gaps_ratio.py @ 0:b15a3147e604 draft
"planemo upload for repository https://github.com/inab/trimal commit cbe1e8577ecb1a46709034a40dff36052e876e7a-dirty"
author | padge |
---|---|
date | Fri, 25 Mar 2022 17:10:43 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/trimal_repo/scripts/get_sequences_gaps_ratio.py Fri Mar 25 17:10:43 2022 +0000 @@ -0,0 +1,78 @@ +#!/usr/bin/python3 + +# +# 'get_sequneces_gaps.py' +# +# Script implemented to obtain the sequences index for those seuqneces +# exceding a minimum gaps' ratio threshold. +# +# [2020] S. Capella-Gutierrez - salvador.capella@bsc.es +# +# this script is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free +# Software Foundation, the last available version. +# +# this script is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for +# more details on <http://www.gnu.org/licenses/> +# +from Bio import AlignIO +import argparse +import sys +import os + +if __name__ == "__main__": + + parser = argparse.ArgumentParser() + + parser.add_argument("-i", "--in", dest = "inFile", required = True, type = \ + str, help = "Input alignment") + + parser.add_argument("-o", "--out", dest = "outFile", default = None, type = \ + str, help = "Set output file. It will be generated into FASTA format") + + parser.add_argument("-f", "--format", dest = "inFormat", default = "fasta", \ + type = str, choices = ["clustal", "fasta-m10", "fasta", "phylip-relaxed", \ + "phylip-sequential", "phylip", "nexus"],help = "Set input alignment format") + + parser.add_argument("-g", "--gap_symbol", dest = "gapSymbol", default = '-', \ + type = str, help = "Define the gap symbol used in the input alignment") + + parser.add_argument("--show_only_index", dest = "showIndexes", default = False, \ + action = "store_true", help = "Show only the indexes of sequences with a " + + "gaps' ratio equal or higher than the established threshold") + + parser.add_argument("--threshold", dest = "gapsThreshold", default = 0.0, \ + type = float, help = "Identify sequences with a minimum of gaps' ratio") + + parser.add_argument("--keep_header", dest = "keepHeader", default = False, + action = "store_true", help = "Keep original alignment sequence IDs indepen" + + "dently of blank spaces on it") + + args = parser.parse_args() + + if not os.path.isfile(args.inFile): + sys.exit(("ERROR: Check input alignment file '%s'") % (args.inFile)) + + index = 0 + indexes = [] + ofile = open(args.outFile, "w") if args.outFile else sys.stdout + for record in AlignIO.read(args.inFile, format = args.inFormat): + sequence_id = record.id if not args.keepHeader else record.description + sequence = str(record.seq) + + length = len(sequence) + valid = len([ps for ps in range(length) if sequence[ps] != args.gapSymbol]) + gaps_ratio = 1 - (valid/length) + + if gaps_ratio >= args.gapsThreshold: + if not args.showIndexes: + print(f'{index:4d}\t{sequence_id:30}\t{gaps_ratio:.4f}', file = ofile) + indexes.append(index) + index += 1 + + if args.showIndexes: + print (','.join(map(str, indexes)), file = ofile) + + ofile.close()