trimreads.py

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