updated, production version

Author: icwells

AlignmentProcessor.py

@@ -69,7 +69,7 @@ def makeDir(path, outdir, axt, phylip):
     '''Makes all sub-directories used by program.'''
     os.chdir(outdir)
     for i in ["01_splitFastaFiles", "02_rmHeader", "03_checkFrame",
-              "04_countBasesPercent", "05_ReplaceStopCodons"]:
+              "04_countBasesPercent", "05_ReplaceStopCodons", "Logs"]:
         try:
             os.mkdir(i)
         except FileExistsError:
@@ -142,8 +142,12 @@ def countBases(outdir, percent):
 def replaceStop(outdir, retainStops):
     '''Removes stop codons for downstream analysis.'''
     print("Removing stop codons...")
-    rs = Popen(split("python bin/05_ReplaceStopCodons.py " + " "
-                     + outdir + " " + str(retainStops)))
+    if retainStops == False:
+        rs = Popen(split("python bin/05_ReplaceStopCodons.py " + " "
+                     + outdir))
+    elif retainStops == True:
+        rs = Popen(split("python bin/05_ReplaceStopCodons.py " + " "
+                     + outdir + " retainStops"))
     rs.wait()
     if rs.returncode == 0:
         return True
@@ -192,8 +196,7 @@ def phylipConvert(outdir, starttime, codeml):
 def runcodeml(outdir, starttime):
     '''Runs codeml on a directory.'''
     print("Running codeml...")
-    cm = Popen(split("python bin/07_CodeMLonDir.py " + outdir + "codeml.ctl "
-                     + outdir))
+    cm = Popen(split("python bin/07_CodeMLonDir.py " + outdir))
     cm.wait()
     if cm.returncode == 0:
         elapsedtime = datetime.now() - starttime
@@ -207,6 +210,10 @@ def helplist():
     print("    example usage: python AlignmentProcessor.py -% <decimal> \
 --axt/phylip --kaks/codeml -i <input fasta file> -o \
 <path to output directory> -r <reference species>")
+    print("    --printNameList    prints list of genome build names and\
+ associated common names")
+    print("    --addNameToList    add new genome build name and\
+ associated common name to list")
     print("    --axt    Converts files to axt for use in KaKs_Calcuator.")
     print("    --phylip   Converts files to phylip for use in PhyML.")
     print("    --kaks     Runs KaKs_Calcuator if --axt is also specified")
@@ -236,6 +243,7 @@ def main():
     codeml = False
     conv = False
     percent = "0.5"
+    ref = "void"
     # Extract values from command line:
     for i in argv:
         if i == "-h" or i == "--help":
@@ -274,6 +282,11 @@ def main():
         elif i == "--retainStops":
             retainStops = True
     # Check inout commands prior to running:
+    if ref == "void":
+        print()
+        print("\tPlease pecify a reference species.")
+        print()
+        quit()
     checkInput(axt, kaks, phylip, codeml)
     # Save working directory to variable to call other scripts:
     path = os.getcwd()

AlignmentProcessorReadMe.txt

@@ -18,6 +18,8 @@
 #			Python 3 version of Biopython
 #			Perl
 #			PAML 
+#			R
+#			ape R package
 ###############################################################
 
 ### Contents ###
@@ -32,10 +34,16 @@
 # 0. Introduction
 #-------------------------------
 
-AlignmentProcessor is a pipeline mean to quickly convert a multi-fasta 
+AlignmentProcessor is a pipeline meant to quickly convert a multi-fasta 
 alignment file into a format that can be read by KaKs_Calculator or PAML and 
-optionally run those programs. Each python script can be run individually, or 
-you can run the AlignmentProcessor wrapper which will call all of the files.
+optionally run those programs. When running codeml, alignment processor will
+use input control and tree files as templates to create unique control and
+tree files for each gene. It will call the ape R package to dynamically trim
+the given phylogenic tree so that only species which remain in each gene's 
+alignment after trimming are represented in that gene's tree. 
+
+You can run the AlignmentProcessor wrapper which will call all of the python 
+scripts in sequence, or you may call each script individually.
 
 # Installing Biopython
 
@@ -66,12 +74,26 @@ If you plan to use CodeML, you must first download PAML
 (http://abacus.gene.ucl.ac.uk/software/paml.html) and move the folder into the
 AlignmentProcessor directory. Make sure that it is titled "paml".
 
+# Ape
+
+The most straightforward way to install ape, and most R packages, is through 
+Bioconductor. If you do not have Bioconductor installed, open R and paste:
+
+	source("https://bioconductor.org/biocLite.R")
+	biocLite()
+
+To install ape, enter:
+
+	library("BiocInstaller")
+	biocLite("ape")
+
+
 #-------------------------------
 # 1. Obtaining a fasta alignment
 #-------------------------------
 # UCSC Fasta Alignment
 It is possible to download CDS fasta alignments from the UCSC Table browser.
-This does, unfortunately, limit you to the alignments that they have available.
+This does, unfortunately, limit you to currently available alignments.
 If they do have the alignment that you are interested in, however, it will 
 probably be faster to download it rather than generate a new one. Since this
 precludes the use of user-generated alignments, AlignmentProcessor has been 
@@ -125,18 +147,16 @@ run the substitutions quickly, so everything can be run with one command. Each
 script can be run individually if necessary (each script's function and 
 options will be discussed later). The input order for AlignmentProcessor's 
 command line does not matter, but it does matter for the individual scripts.
-AlignmentProcessor also comes packaged with KaKs_Calculator2.0 and PAML4.8 for
-downstream analysis.
 
 To execute the AlignmentProcessor pipeline, you must first change into
 the package directory. Otherwise it will not be able to locate the scripts
 in the bin/ directory.
 
 # Example Usage: 
 
-	python AlignmentProcessor.py -% <decimal> --retainStops \
-		--axt/phylip --kaks/codeml --ucsc -i <input fasta file> \
-		-o <path to output directory> -r <reference species>
+	python AlignmentProcessor.py --ucsc --axt/phylip --kaks/codeml \
+		--retainStops -% <decimal> -r <reference species> -i <input fasta file> \
+		-o <path to output directory> 
 
 # Required Arguments:
 
@@ -148,7 +168,8 @@ in the bin/ directory.
 
 	--ucsc	This will invoke 00_convertHeader.py, which will convert the 
 		headers from UCSC fasta files so they only contain build
-		names and gene IDs. 
+		names and gene IDs. This does not need to be run on Stich Gene
+		Blocks output.
 
 	--axt/phylip	Specifies which format to convert the files into 
 			(axt format for KaKs_Calculator; phylip for CodeML 
@@ -163,13 +184,12 @@ in the bin/ directory.
 			control file for CodeML which must be located in the 
 			output directory you specified with the the -o option.
 			This file must be titled "codeml.ctl" (the default 
-			name given by PAML. The name is hard coded into
-			AlignmentProcessor to simplify the commands). 
+			name given by PAML). 
 
 	--retainStops	This will tell AlignmentProcessor to retain sequences
 			that contain internal stop codons. By default,
 			sequences with internal stop codons will be removed 
-			from the analysis.
+			from the analysis as they may bias the results.
 
 	-%	a decimal value specifying the minimum percentage of reads 
 		that must remain after replacing unknown codons with gaps 
@@ -184,10 +204,11 @@ in the bin/ directory.
 
 	--printNameList	will print the contents of 02_nameList.txt which 
 			contains the list of genome builds and associated 
-			common names.
+			common names. Must be run without any other arguments
 
 	--addNameToList will add an entry to the 02_nameList.txt file.
-e.g. python AlignmentProcessor.py -- addNameToList <build> <common name>
+			e.g. python AlignmentProcessor.py -- addNameToList \
+				<build> <common name>
 
 
 # Genome Builds and Common Names
@@ -201,7 +222,7 @@ e.g. python AlignmentProcessor.py -- addNameToList <build> <common name>
 	new row, followed by a tab, then the desired common name. Make sure 
 	there are no spaces in either name.
 
-# The codeml control file
+# The CodeML control file
 
 	Codeml requires that all of its parameters be specified in one control 
 	file (http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf). Provide a 
@@ -211,7 +232,20 @@ e.g. python AlignmentProcessor.py -- addNameToList <build> <common name>
 	tree file for codeml (see PAML manual above).
 
 	The control file must be titled titled “codeml.ctl”, and it must be 
-	locatedin the output directory.
+	located in the output directory.
+
+# The CodeML tree file
+
+	If your CodeML analysis requires a phylogneic tree, provide your 
+	desired tree, titled "codeml.tree" in the output directory. Specify
+	the tree as you want to appear to CodeML and be sure that the species
+	names are specified as the common names in the 02_nameList.txt file, 
+	but trimmed to ten characters (some programs still set a ten character
+	limit on the length of names, so AlignmentProcessor trims the names). 
+	The 07_CodeMLonDir.py script will save any nodes you have specified 
+	before sending a plain Newick tree to ape (which will not work if 
+	there are PAML node symbols). It will then add any nodes back into the
+	tree after it has been trimmed.
 
 # Invoking the Ka/Ks pipeline with a UCSC alignment:
 
@@ -231,8 +265,8 @@ Each script performs one or two functions on the input file or files, and
 saves the output to a new subdirectory. The location of the working directory
 must be specified, but the subdirectories are hard-coded in the programs. 
 Be sure to include the trailing "/" in the path. The AlignmentProcessor
-wrapper will add it, but, to avoid redundancies, the individual scripts
-will not.
+wrapper will add it, but to avoid redundancies the individual scripts will 
+not.
 
 Remember that the order of the arguments does matter for these scripts.
 
@@ -246,17 +280,17 @@ Remember that the order of the arguments does matter for these scripts.
 # 01_SplitFastaFiles.py
 
 	This script will split the input mult-fasta alignment into one file
-	per gene.
+	per gene. It will produce an output file for a gene if it has at least
+	two sequences. 
 
-	python 01_splitFastaFiles.py <# of species in alignment> \
-		<input fasta alignment> <path to output directory>
+	python 01_splitFastaFiles.py <input fasta alignment> \
+		<path to output directory>
 
 # 02_RemoveHeader.py
 
 	This program will read through a directory that contains 
-	aligned multiple FASTA files and remove FASTA header and all 
-	the trailing information for each gene and change all the 
-	genome build names to common names.
+	aligned multiple FASTA files and replace FASTA headers with each 
+	species' common name.
 
 	python 02_RemoveHeaderOnDir.py 	<path to inut and output directories>
 
@@ -269,21 +303,20 @@ Remember that the order of the arguments does matter for these scripts.
 	frame. It will then replaces codons with missing nucleotides with gaps
 	to remove unknown amino acids from the sequence.
 
-	python 03_CheckFrameOnDir.py <# of species> \
-		<path to inut and output directories> <reference_species>
+	python 03_CheckFrameOnDir.py <path to inut and output directories> \
+		<reference_species>
 
 # 04_CountBases.py
 
 	This program will check a multiple FASTA file to see that each species
  	retains a certain percentage of its nucleotide sequence. If not, it 
 	will remove that sequence.
 
-	Note: the AlignmentProcessor wrapper specifies 50% as adefault value, 
+	Note: the AlignmentProcessor wrapper specifies 50% as a default value, 
 	but the script itself does not, so you must specify one if you invoke 
 	it on its own.
 
-	python 05_CountBasesOnDir.py <number of species> \
-		<threshold percentage as a decimal> \
+	python 05_CountBasesOnDir.py <threshold percentage as a decimal> \
 		<path to inut and output directories>
 
 # 05_ReplaceStopCodons.py
@@ -294,16 +327,18 @@ Remember that the order of the arguments does matter for these scripts.
 	codon.
 
 	Terminal stop codons will be replaced, while sequences with internal 
-	stop codons will have their gene id, as well sequence name and 
-	location of the first occuring stop codon, recorded in the 
-	internalStops.txt file. 
+	stop codons will have their gene id and sequence name recorded in the 
+	internalStops.txt file. If --retainStops is specified, these sequences
+	will be retained. If it is not, these sequences will be removed and 
+	any which does not have at least two remaining sequences will not be 
+	written to file.
 
-	python 05_ReplaceStopCodonsOnDir.py <number of species> \
-		<path to inut and output directories>
+	python 05_ReplaceStopCodonsOnDir.py \
+		<path to inut and output directories> --retainStops(optional)
 
 # 06_FASTAtoAXT.py
 
-	This program executes 06_parseFastaIntoAXT.pl on an entire direcotry,
+	This program executes 06_parseFastaIntoAXT.pl on an entire directory,
 	allowing all of the contents of the directory to be converted to 
 	axt files.
 
@@ -322,7 +357,7 @@ Remember that the order of the arguments does matter for these scripts.
 
 # 07_KaKsonDir.py
 
-	This program executes KaKs_Calculator on every file in a direcotry. 
+	This program executes KaKs_Calculator on every file in a directory. 
 
 	python 07_KaKsonDirectory.py <path to inut and output directories> \
 		<name of refernce species>
@@ -331,10 +366,22 @@ Remember that the order of the arguments does matter for these scripts.
 
 	This script will run codeml on every file in a directory. It requires
 	the codeml.ctl file, and likely a tree file which it will supply to
-	codeml.
+	codeml. It will overwrite the "seqfile", "treefile", "outfile" lines 
+	include the paths to the input phylip file, the output file, and the 
+	tree file. It will also call the ape R package to trim the tree file 
+	so that it only includes species which have not been filtered out.
 
 	python 07_CodeMLonDir.py <path to codeml control file> \
-		<path to input and output directorie>
+		<path to input and output directories>
+
+3 07_pruneTree.R
+
+	This R script will call the ape package to dynamically trim input 
+	trees for CodeML if any sequences have been removed. Species 
+	whose sequences were removed in steps 4 or 5 will be removed from
+	the temporary tree given to CodeML.
+
+	(Caled by 07_CodeMLonDir.py)
 
 # 08_compileKaKs_CSV.py
 
@@ -361,26 +408,37 @@ printed to the output directory that you specified with the -o option.
 
 If --phylip is specified, the program will convert the fasta files to phylip
 after the stop codons have been removed. If you also specified --codeml, 
-AlignmentProcessor will edit and submit the control file to codeml and the 
+AlignmentProcessor will edit and submit the control file to CodeML and the 
 output files will be saved in 07_codeml. Since there are many different things 
 that can be done with the codeml output files, AlignmentProcessor does not 
 attmept to concatenate specific parts from the output files.
 
-If you wish to convert convert the files to both formats, with or without the 
---kaks option, specify both --axt and --phylip the program will convert the 
-fasta files to both formats. You may also run one of the individual scripts on
-the 07_rmStops directory to convert the files in a separate step. 
-AlignmentProcessor will not, however, run KaKs_Calculator and codeml 
-simultaneously, as this could require too much memory.
+If you wish to convert convert the files to both formats, specify both --axt
+and --phylip the program will convert the fasta files to both formats. You may
+also run one of the individual scripts on the 07_rmStops directory to convert
+the files in a separate step. AlignmentProcessor will not, however, run 
+KaKs_Calculator and CodeML simultaneously, as this could require too much 
+memory.
 
 #-------------------------------
 # 5. Run AlignmentProcessor on test data
 #-------------------------------
 
+# To test KaKs_Calculator:
 Change directory into the AlignmentProcessor folder. Paste the followig into
-a terminal (change the output directory to the desired loaction):
+a terminal:
 
 python AlignmentProcessor.py --axt --kaks --ucsc -r Green_anole \
--i test.fa -o test/
+-i kaksTest.fa -o test/
 
 This will return a text file with 11 lines.
+
+# To test CodeML:
+The test directory already contains sample CodeML control and tree files, so
+all you need  to do is change into the AlignmentProcessor direcotry and paste
+the following:
+
+python AlignmentProcessor.py --phylip --codeml --ucsc -r Green_anole \
+-i codemlTest.fa -o test/
+
+There should be 8 .mlc files in the 07_codeml directory.