Mercurial > repos > peterjc > fasta_filter_by_id

annotate tools/fasta_tools/fasta_filter_by_id.py @ 2:5b552b3005f2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Migrated tool version 0.0.4 from old tool shed archive to new tool shed repository
author peterjc
date Tue, 07 Jun 2011 17:23:07 -0400
parents 5cd569750e85
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
1 #!/usr/bin/env python
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
2 """Filter a FASTA file with IDs from a tabular file, e.g. from BLAST.
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
3
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
4 Takes five command line options, tabular filename, ID column numbers
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
5 (comma separated list using one based counting), input FASTA filename, and
1
5cd569750e85 Migrated tool version 0.0.3 from old tool shed archive to new tool shed repository
peterjc
parents: 0
diff changeset
6 two output FASTA filenames (for records with and without the given IDs).
5cd569750e85 Migrated tool version 0.0.3 from old tool shed archive to new tool shed repository
peterjc
parents: 0
diff changeset
7
5cd569750e85 Migrated tool version 0.0.3 from old tool shed archive to new tool shed repository
peterjc
parents: 0
diff changeset
8 If either output filename is just a minus sign, that file is not created.
0
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
9 This is intended to allow output for just the matched (or just the non-matched)
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
10 records.
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
11
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
12 Note in the default NCBI BLAST+ tabular output, the query sequence ID is
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
13 in column one, and the ID of the match from the database is in column two.
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
14 Here sensible values for the column numbers would therefore be "1" or "2".
2
5b552b3005f2 Migrated tool version 0.0.4 from old tool shed archive to new tool shed repository
peterjc
parents: 1
diff changeset
15
5b552b3005f2 Migrated tool version 0.0.4 from old tool shed archive to new tool shed repository
peterjc
parents: 1
diff changeset
16 This script is copyright 2010 by Peter Cock, SCRI, UK. All rights reserved.
5b552b3005f2 Migrated tool version 0.0.4 from old tool shed archive to new tool shed repository
peterjc
parents: 1
diff changeset
17 See accompanying text file for licence details (MIT/BSD style).
5b552b3005f2 Migrated tool version 0.0.4 from old tool shed archive to new tool shed repository
peterjc
parents: 1
diff changeset
18
5b552b3005f2 Migrated tool version 0.0.4 from old tool shed archive to new tool shed repository
peterjc
parents: 1
diff changeset
19 This is version 0.0.3 of the script.
0
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
20 """
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
21 import sys
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
22 from galaxy_utils.sequence.fasta import fastaReader, fastaWriter
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
23
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
24 def stop_err( msg ):
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
25 sys.stderr.write( msg )
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
26 sys.exit()
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
27
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
28 #Parse Command Line
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
29 try:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
30 tabular_file, cols_arg, in_file, out_positive_file, out_negative_file = sys.argv[1:]
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
31 except ValueError:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
32 stop_err("Expected five arguments, got %i:\n%s" % (len(sys.argv)-1, " ".join(sys.argv)))
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
33 try:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
34 columns = [int(arg)-1 for arg in cols_arg.split(",")]
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
35 except ValueError:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
36 stop_err("Expected list of columns (comma separated integers), got %s" % cols_arg)
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
37
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
38 #Read tabular file and record all specified identifiers
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
39 ids = set()
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
40 handle = open(tabular_file, "rU")
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
41 if len(columns)>1:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
42 #General case of many columns
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
43 for line in handle:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
44 if line.startswith("#"):
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
45 #Ignore comments
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
46 continue
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
47 parts = line.rstrip("\n").split("\t")
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
48 for col in columns:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
49 ids.add(parts[col])
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
50 print "Using %i IDs from %i columns of tabular file" % (len(ids), len(columns))
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
51 else:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
52 #Single column, special case speed up
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
53 col = columns[0]
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
54 for line in handle:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
55 if not line.startswith("#"):
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
56 ids.add(line.rstrip("\n").split("\t")[col])
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
57 print "Using %i IDs from tabular file" % (len(ids))
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
58 handle.close()
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
59
1
5cd569750e85 Migrated tool version 0.0.3 from old tool shed archive to new tool shed repository
peterjc
parents: 0
diff changeset
60 #Write filtered FASTA file based on IDs from tabular file
0
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
61 reader = fastaReader(open(in_file, "rU"))
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
62 if out_positive_file != "-" and out_negative_file != "-":
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
63 print "Generating two FASTA files"
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
64 positive_writer = fastaWriter(open(out_positive_file, "w"))
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
65 negative_writer = fastaWriter(open(out_negative_file, "w"))
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
66 for record in reader:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
67 #The [1:] is because the fastaReader leaves the > on the identifer.
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
68 if record.identifier and record.identifier.split()[0][1:] in ids:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
69 positive_writer.write(record)
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
70 else:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
71 negative_writer.write(record)
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
72 positive_writer.close()
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
73 negative_writer.close()
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
74 elif out_positive_file != "-":
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
75 print "Generating matching FASTA file"
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
76 positive_writer = fastaWriter(open(out_positive_file, "w"))
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
77 for record in reader:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
78 #The [1:] is because the fastaReader leaves the > on the identifer.
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
79 if record.identifier and record.identifier.split()[0][1:] in ids:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
80 positive_writer.write(record)
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
81 positive_writer.close()
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
82 elif out_negative_file != "-":
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
83 print "Generating non-matching FASTA file"
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
84 negative_writer = fastaWriter(open(out_negative_file, "w"))
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
85 for record in reader:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
86 #The [1:] is because the fastaReader leaves the > on the identifer.
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
87 if not record.identifier or record.identifier.split()[0][1:] not in ids:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
88 negative_writer.write(record)
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
89 negative_writer.close()
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
90 else:
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
91 stop_err("Neither output file requested")
2e5f8ad1a096 Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
peterjc
parents:
diff changeset
92 reader.close()