randseq

#!/usr/bin/env python2.7
#
# Neat little generator to generate random reads
# from a nucleotide file
#
# Author: Jordan Texier,
# Created: Jan 2014
# Modified: Oct 2014
# 
# DittmerLab, UNC Chapel Hill
#

def usage():
	"""Prints the correct usage of this program"""

	str="usage :" + sys.argv[0] + " -h "
	print str


def getRandString(lines, stringsize, randseed):
	"""Return a random read of size 'stringsize' from the input file"""

	try:
		linenum=randseed/length
		start=randseed%length

		# If first line, go to second
		if linenum == 0:
			linenum=1
			start=0
		
		# Gather the minimum amount of lines to extract
		Nlines = (stringsize/length)+2
		if Nlines < 2:
			Nlines = 2

		templines = lines[linenum:(linenum+Nlines)]

		# Merge multiple lines together
		bigline = ''.join(templines)
		bigline = bigline.replace('\n','')

		if len(bigline) < stringsize:
			raise IndexError('Probably reached end of file')

		# Extract read from multiline
		# print bigline
		read = bigline[start:(start + stringsize)]

		# Check if line is actually reads (first 2 char)
		if read[0] not in ['A','G','C','T'] or read[1] not in ['A','C','G','T']:
			read = 'N'


		return read.strip()

	except IndexError:
		# print "IndexError"
		return 'N'




def isNstring(string):
	"""Return true if first or last char is 'N' (implying there are repeating N's)"""
	if (string is ''):
		return False
	else:
		return ( string[0] is 'N' or string[len(string) -1] is 'N')




def getTotalNuc(linelist):
	"""Return the total number of nucleotide bases in the file, along with the length of each line"""

	N = (len(linelist) + 1) * len(linelist[1])
	return N, len(lines[1]) 





# ------------------------------ #
# ============ Main ============ #
# ------------------------------ #

if __name__ == "__main__":

	import sys
	import argparse
	import random

	# Create Parser, required and optional args. Then parse
	parser=argparse.ArgumentParser(description='Generate random reads from genome file.', 
					prog='randseq', usage='%(prog)s -i inputfile [options]',
					formatter_class=lambda prog: argparse.HelpFormatter(prog, max_help_position = 50, width = 200),
					add_help=False)

	required = parser.add_argument_group('Required Arguments')
	required.add_argument("-i", "--input", default='', help="input file")

	optional = parser.add_argument_group('Optional Arguments')
	optional.add_argument('-h', '--help', action='help', help='this very helpful text')
	optional.add_argument('-H', '--header', action='store_true',  help='enable header in output fasta')
	optional.add_argument('-p', '--prefix', default='>', help='character prefixing header files (requires -H)')
	optional.add_argument('-o', '--output', default='', help='output file, if empty prints to console')
	optional.add_argument('-l', '--readslength', type=int, default=64 ,help='length of individual reads (default=64)')
	optional.add_argument('-n', '--numreads', type=int, default=1000, help='number of random reads (default=1000)')
	
	args = vars( parser.parse_args() )




	# If required arguments are empty
	# print help and bail out	
	if(args["input"] is ''):
		parser.print_help()
		sys.exit()
	if(args['prefix'] != '>' and not args['header'] ):
		print "cant prefix headers without -H tag"
		parser.print_help()
		sys.exit()



	# Get properties of file
	infile = open( args["input"], 'r')
	lines = infile.readlines()[1:]
	(total, length) = getTotalNuc( lines )
	total = total - args["readslength"]
	

	# Define output
	if( args['output'] is ''):
		out = sys.stdout
	else:
		out = open(args['output'], 'w')

	# Header and no-header loops:
	if( not args['header']):
		i = 0

		while( i < args['numreads']):
			str = 'N'
			while( isNstring(str)):
				r = random.randint(1,total)
				str = getRandString(lines, args['readslength'], r)
			
			out.write(str + '\n')
			i+=1
	else:
		i = 0
		
		while( i < args['numreads']):
			str = 'N'
			# Modify this line to change Header format:
			header = args['prefix'] + "Random_reads-" + args['input'] + "_" + `i`

			while( isNstring(str)):
				r = random.randint(1,total)
				str = getRandString(lines, args['readslength'], r)
			
			out.write(header + '\n')
			out.write(str + '\n')
			i+=1

		out.close()


	# Closing open files
	if (not args['output'] is ''):
		out.close()

	infile.close()