trimreads.py

   1 #
   2 #  trimquery.py
   3 #  ENRAGE
   4 #
   5 #  Created by Ali Mortazavi on 8/12/08.
   6 #
   7
   8 import sys
   9 import optparse
  10 from cistematic.core import complement
  11 from commoncode import getConfigParser, getConfigBoolOption, getConfigOption
  12
  13 print "trimreads: version 2.2"
  14
  15 def main(argv=None):
  16     if not argv:
  17         argv = sys.argv
  18
  19     usage = "usage: python %prog length infile outfile [--fastq] [--fromback] [--paired] [--flip] [--filter maxN]"
  20
  21     parser = getParser(usage)
  22     (options, args) = parser.parse_args(argv[1:])
  23
  24     if len(args) < 3:
  25         print usage
  26         print "\t where paired fragments are separated by a : when given the -paired flag"
  27         sys.exit(1)
  28
  29     length = int(args[0])
  30     infile = args[1]
  31     outfile = args[2]
  32
  33     trimreads(length, infile, outfile, options.fastq, options.fromBack, options.paired, options.flipseq, options.maxN)
  34
  35
  36 def getParser(usage):
  37     parser = optparse.OptionParser(usage=usage)
  38     parser.add_option("--fastq", action="store_true", dest="fastq")
  39     parser.add_option("--fromback", action="store_true", dest="fromBack")
  40     parser.add_option("--paired", action="store_true", dest="paired")
  41     parser.add_option("--flip", action="store_true", dest="flipseq")
  42     parser.add_option("--filter", type="int", dest="maxN")
  43
  44     configParser = getConfigParser()
  45     section = "trimreads"
  46     fastq = getConfigBoolOption(configParser, section, "fastq", False)
  47     fromBack = getConfigBoolOption(configParser, section, "fromBack", False)
  48     paired = getConfigBoolOption(configParser, section, "paired", False)
  49     flipseq = getConfigBoolOption(configParser, section, "flipseq", False)
  50     maxN = getConfigOption(configParser, section, "maxN", None)
  51
  52     parser.set_defaults(fastq=fastq, fromBack=fromBack, paired=paired, flipseq=flipseq, maxN=maxN)
  53
  54     return parser
  55
  56
  57 def trimreads(length, inFileName, outFileName, fastq=False, fromBack=False, paired=False, flipseq=False, maxN=None):
  58     infile = open(inFileName)
  59     outfile = open(outFileName, "w")
  60
  61     if paired:
  62         pairedlength = 2 * length
  63     index = 0
  64
  65     if fromBack:
  66         length = -1 * length
  67
  68     filtering = False
  69     if maxN is not None:
  70         filtering = True
  71         print "filtering out reads with more than %d Ns" % maxN
  72     else:
  73         maxN = 2
  74
  75     print "trimming reads from %s to %d bp and saving them in %s" % (inFileName, length, outFileName)
  76
  77     filtered = 0
  78     header = ""
  79     for line in infile:
  80         line = line.strip()
  81         if len(line) == 0:
  82             continue
  83
  84         firstChar = line[0]
  85         if (not fastq and firstChar == ">") or (fastq and firstChar in ["@", "+"]):
  86             header = line + "\n"
  87         else:
  88             if filtering:
  89                 if line.count("N") > maxN:
  90                     filtered += 1
  91                     continue
  92
  93             seq1 = line[length:]
  94             seq2 = line[:length]
  95             if flipseq:
  96                 try:
  97                     tempseq1 = seq1
  98                     seq1 = complement(tempseq1)
  99                 except:
 100                     seq1 = tempseq1
 101
 102                 try:
 103                     tempseq2 = seq2
 104                     seq2 = complement(tempseq2)
 105                 except:
 106                     seq2 = tempseq2
 107
 108             if paired:
 109                 if len(line) < pairedlength:
 110                     continue
 111
 112                 outfile.write("%s%s:%s\n" % (header, seq1, seq2))
 113             else:
 114                 if fromBack:
 115                     outfile.write("%s%s\n" % (header, seq1))
 116                 else:
 117                     outfile.write("%s%s\n" % (header, seq2))
 118
 119             index += 1
 120             if index % 1000000 == 0:
 121                 print ".",
 122
 123             sys.stdout.flush()
 124
 125     outfile.close()
 126     print "returned %d reads" % index
 127     if filtering:
 128         print "%d additional reads filtered" % filtered
 129
 130
 131 if __name__ == "__main__":
 132     main(sys.argv)