cgeb2001's DokuWiki!

User Tools

Site Tools

Trace:
multi-gene_phylogeny_pipeline

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
multi-gene_phylogeny_pipeline [2017/06/19 08:23] – cgeb2001multi-gene_phylogeny_pipeline [2018/03/10 11:07] (current) – 173.212.69.201
Line 1: Line 1:
-**Multi-gene phylogeny tree using Matt Brown’s Bordor dataset and pipeline**+====== Multi-gene phylogeny tree using Matt Brown’s Bordor dataset and pipeline ====== 
    
 Documentation by Kate Glennon, Sarah Shah, Shelby Williams, and Tommy Harding. Documentation by Kate Glennon, Sarah Shah, Shelby Williams, and Tommy Harding.
    
-The **Bordor** dataset is a set of 351 housekeeping genes that are well-conserved across all eukaryotes. This pipeline uses the gene sequences from //Arabidopsis thaliana// as queries to fish for homologues during the BLAST step.+The **Bordor** dataset is a set of 351 housekeeping genes that are well-conserved across all eukaryotes. This pipeline uses the gene sequences from //Arabidopsis thaliana// as queries to fish for homologues during the BLAST step. Credit goes to Matt Brown & co.: [[https://doi.org/10.1093/gbe/evy014]], [[https://doi.org/10.1093/molbev/msx162]] 
 + 
 +All the original transcriptom/proteome files that are in the Bordor alignment is in **/scratch2/mbrown/PhylogenomicDatabases**
    
 **The Pipeline Overview** **The Pipeline Overview**
Line 78: Line 81:
 /scratch2/mbrown/START.02-09-2017.SHARE.tgz This has more stramenopiles than other versions. /scratch2/mbrown/START.02-09-2017.SHARE.tgz This has more stramenopiles than other versions.
 Additional scripts listed in the following instructions can be copied from: Additional scripts listed in the following instructions can be copied from:
-/scratch2/sarahshah/START.15.06.17ss+/home/sarahshah/AdditionalScripts/Multigene_Phylogeny
  
  
Line 94: Line 97:
    
 Step 3: If you need to add sequences from only one taxon, do the following, otherwise see NOTE at the end of step 3. Edit Bordor.sh and change the input file name to yours, e.g. : Step 3: If you need to add sequences from only one taxon, do the following, otherwise see NOTE at the end of step 3. Edit Bordor.sh and change the input file name to yours, e.g. :
- +<code>
 #!/bin/bash #!/bin/bash
- 
 #$ -S /bin/sh #$ -S /bin/sh
- 
 . /etc/profile . /etc/profile
- 
 #$ -o matt_pipeline.info.txt #$ -o matt_pipeline.info.txt
- 
 #$ -cwd #$ -cwd
- 
 #$ -pe threaded 8 #$ -pe threaded 8
    
-python AddPipeline3.0a.py <short name> <complete species name with “_” in between the genus and species: Genus_species> 1 <AA or NUC> Bordor.351.refdat.txt ./ END.YY-MM-DD yes +python AddPipeline3.0a.py <short name> <complete species name with only one “_” in between the genus and species: Genus_species**> 1 <AA or NUC> Bordor.351.refdat.txt ./ END.YY-MM-DD yes 
-  +</code> 
-The number “1” refers to the standard genetic code. Use “NUC” if your fasta file contains nucleotide sequences, or change it to “AA” for protein sequences. Say “yes” for the last flag at the end of the line if you want your alignments to be trimmed by bmge. Edit the date attached to “END*” to match today’s date. Make sure the AddPipeline3.X.py in your “START*” folder matches the one in this shell script. As of now, AddPipeline3.0a.py is the latest version. Then qsub Bordor.sh+The number “1” refers to the standard genetic code. Use “NUC” if your fasta file contains nucleotide sequences, or change it to “AA” for protein sequences. Say “yes” for the last flag at the end of the line if you want your alignments to be trimmed by bmge. Edit the date attached to “END*” to match today’s date. Make sure the AddPipeline3.X.py in your “START*” folder matches the one in this shell script. As of now, AddPipeline3.0a.py is the latest version. **Make sure your "short name" and "long name" are correct, i.e. 8 characters for the former, and the latter must have one "_"**. Then qsub Bordor.sh. Ensure that you are using a node with enough CPU's available, otherwise your Bordor.sh error file will show that the threadcount is out of range.
    
 NOTE: If you need to add sequences from several taxa: in step 1, instead of renaming the “END*” folder “START*”, rename it with the short name of the first organism you want to add; and in step 2, copy original fasta files for all the organisms of interest (renamed with the organism short names as instructed in step 2 above) upstream of the folder created in step 1. Save the Bordor.sh script at the same location and edit it as follows (for 3 organisms as example): NOTE: If you need to add sequences from several taxa: in step 1, instead of renaming the “END*” folder “START*”, rename it with the short name of the first organism you want to add; and in step 2, copy original fasta files for all the organisms of interest (renamed with the organism short names as instructed in step 2 above) upstream of the folder created in step 1. Save the Bordor.sh script at the same location and edit it as follows (for 3 organisms as example):
- +<code>
 #!/bin/bash #!/bin/bash
- 
 #$ -S /bin/sh #$ -S /bin/sh
- 
 . /etc/profile . /etc/profile
- 
 #$ -o matt_pipeline.info.txt #$ -o matt_pipeline.info.txt
- 
 #$ -cwd #$ -cwd
- 
 #$ -pe threaded 8 #$ -pe threaded 8
    
 cp Genuspec1.fas Genuspec1/ cp Genuspec1.fas Genuspec1/
- 
 cd Genuspec1/ cd Genuspec1/
- 
 python AddPipeline3.0a.py Genuspec1 Genus_species1 1 AA Bordor.351.refdat.txt ./ Genuspec2 no python AddPipeline3.0a.py Genuspec1 Genus_species1 1 AA Bordor.351.refdat.txt ./ Genuspec2 no
    
 mkdir ../Genuspec2 mkdir ../Genuspec2
- 
 cd ../Genuspec2/ cd ../Genuspec2/
- 
 python AddPipeline3.0a.py Genuspec2 Genus_species2 1 AA Bordor.351.refdat.txt ./ Genuspec3 no python AddPipeline3.0a.py Genuspec2 Genus_species2 1 AA Bordor.351.refdat.txt ./ Genuspec3 no
    
 mkdir ../Genuspec3 mkdir ../Genuspec3
- 
 cd ../Genuspec3/ cd ../Genuspec3/
- 
 python AddPipeline3.0a.py Genuspec3 Genus_species3 1 AA Bordor.351.refdat.txt ./ END.YY-MM-DD yes python AddPipeline3.0a.py Genuspec3 Genus_species3 1 AA Bordor.351.refdat.txt ./ END.YY-MM-DD yes
- +</code>
 This will sequentially add the appropriate sequences for all the organisms of interest to the Bordor dataset. Trimming will not occur until the last taxon is added. This will sequentially add the appropriate sequences for all the organisms of interest to the Bordor dataset. Trimming will not occur until the last taxon is added.
- + 
 +NOTE2: If you have alignment files from someone else, and you want to add your own transcriptomes to them, move the alignment files in the folder "old_aln" in your START folder. 
 Step 4: If everything went as expected, there will be a folder named “bmge_trimmed_old” in the “END*” folder. Download a bunch of *.faa (aligned non-trimmed sequences) and *.bmge.fas (trimmed aligned sequences) files to your computer and examine them with a sequence viewer such as AliView. The last line(s) is the sequence from your transcriptome/protein data that was aligned to the other sequences of that particular gene. Make sure they look aligned; for instance, if all other sequences have a “GGG” in a specific location then you should expect your sequence to have the same. The .bmge.fas files are .faa files in which the badly-aligned positions were trimmed away.  Step 4: If everything went as expected, there will be a folder named “bmge_trimmed_old” in the “END*” folder. Download a bunch of *.faa (aligned non-trimmed sequences) and *.bmge.fas (trimmed aligned sequences) files to your computer and examine them with a sequence viewer such as AliView. The last line(s) is the sequence from your transcriptome/protein data that was aligned to the other sequences of that particular gene. Make sure they look aligned; for instance, if all other sequences have a “GGG” in a specific location then you should expect your sequence to have the same. The .bmge.fas files are .faa files in which the badly-aligned positions were trimmed away. 
    
Line 150: Line 139:
    
 Step 6: Run the MakeShell.py by typing in:  Step 6: Run the MakeShell.py by typing in: 
 +<code>
 python MakeShell.py <iqtree, RAxML, or FastTree> python MakeShell.py <iqtree, RAxML, or FastTree>
 +</code>
 You should immediately see new .sh files made for each gene. You should immediately see new .sh files made for each gene.
    
 Step 7: This step until Step 9 are for submitting all those tree shell scripts in an efficient way. Make a simple file containing a list of all the .sh just created by typing:  Step 7: This step until Step 9 are for submitting all those tree shell scripts in an efficient way. Make a simple file containing a list of all the .sh just created by typing: 
 +<code>
 ls iq*sh > list_scripts_iqtree (or Fast*sh or RAxML*sh to the name of another tree program) ls iq*sh > list_scripts_iqtree (or Fast*sh or RAxML*sh to the name of another tree program)
- +</code>
 Step 8: Create a text file (called “max_cpu.txt”) that contain the total number CPUs you would like to use when building your phylogenies and another that specifies the number of seconds to wait between submitting each shell script called “sleep_time.txt”. You can do this by: Step 8: Create a text file (called “max_cpu.txt”) that contain the total number CPUs you would like to use when building your phylogenies and another that specifies the number of seconds to wait between submitting each shell script called “sleep_time.txt”. You can do this by:
-  * echo “40” > max_cpu.txt   #40 is the maximum CPUs you’d like to use +<code> 
-  * echo “2” > sleep_time.txt  #The waiting time between submitting each job is 2 seconds +echo “40” > max_cpu.txt   #40 is the maximum CPUs you’d like to use 
- +echo “2” > sleep_time.txt  #The waiting time between submitting each job is 2 seconds 
 +</code> 
 The above files can be edited in real-time; increase the CPUs and sleep time if jobs are slow in being submitted. A good option for FastTree is 40 CPUs and 2 seconds of sleep time. For IQTree, 40 CPUs and 60 seconds of sleep time seem to work. You do not need to cancel your jobs after editing these files. The above files can be edited in real-time; increase the CPUs and sleep time if jobs are slow in being submitted. A good option for FastTree is 40 CPUs and 2 seconds of sleep time. For IQTree, 40 CPUs and 60 seconds of sleep time seem to work. You do not need to cancel your jobs after editing these files.
    
Line 179: Line 169:
 ii) the data related to these taxa were generated by someone else and are not published yet,       iii) taxa irrelevant to your study. You might want to keep a few organisms as outgroups.  ii) the data related to these taxa were generated by someone else and are not published yet,       iii) taxa irrelevant to your study. You might want to keep a few organisms as outgroups. 
 For example, to keep all the taxa that are publicly available, do: For example, to keep all the taxa that are publicly available, do:
-  * grep “YES” Bordor.Taxa.txt > YESlist +<code> 
-  * awk ‘{print $1}’ YESlist  +grep “YES” Bordor.Taxa.txt > YESlist 
 +awk ‘{print $1}’ YESlist  
 +</code>
 Copy the first column that was printed out by the above command and paste it into a file, let’s call this “yourlist”. Copy the first column that was printed out by the above command and paste it into a file, let’s call this “yourlist”.
    
-Step 13: Copy mvtaxatrimmedfas.py, sepalignmask.py, and BMGE.jar from the main START folder into this new folder. Do python taxon_deletion.py yourlist. Note, this script only recognizes a list of the short names in a column format. This will make *taxatrimmed.fas files. Move all of them into a new folder. +Step 13: Copy mvtaxatrimmedfas.py, sepalignmask.py, and BMGE.jar from the main START folder into this new folder. Do: 
 +<code> 
 +python taxon_deletion.py yourlist 
 +</code> 
 +Note, this script only recognizes a list of the short names in a column format. This will make *taxatrimmed.fas files. Move all of them into a new folder. 
    
 Then do: Then do:
-  * python mvtaxatrimmedfas.py. +<code> 
 +python mvtaxatrimmedfas.py  
 +</code>
 This will rename the taxatrimmed.fas files to .fas files.  This will rename the taxatrimmed.fas files to .fas files. 
    
 Next do: Next do:
-  * python sepalignmask.py numberofthreads. +<code> 
 +python sepalignmask.py numberofthreads 
 +</code>
 This will re-align (using MAFFT) and trim (using BMGE) the sequences in these files, producing .faa and .bmge files. This will re-align (using MAFFT) and trim (using BMGE) the sequences in these files, producing .faa and .bmge files.
    
 You can also make shell scripts for all the .py scripts above and qsub them, especially for ones that take some time, such as the sepalignmask.py. You can also make shell scripts for all the .py scripts above and qsub them, especially for ones that take some time, such as the sepalignmask.py.
    
-Step 14: Move all the .bmge.fas files to a new folder. Copy alvert_septable.py and seqtools.py into this folder from the main START folder. Do python alvert_septable.py -c bmge.fas outputnamehere.dat. This will generate a gene table text and a log file along with a phylogenomic supermatrix outputnamehere.dat. Rename the log file before running IQTree as IQTree also generates an output with a log extension with a similar name. The log file lists the number of genes missing from each taxon.+Step 14: Move all the .bmge.fas files to a new folder. Copy alvert_septable.py and seqtools.py into this folder from the main START folder. Do: 
 +<code> 
 +python alvert_septable.py -c bmge.fas outputnamehere.dat 
 +</code> 
 +This will generate a gene table text and a log file along with a phylogenomic supermatrix outputnamehere.dat. Rename the log file before running IQTree as IQTree also generates an output with a log extension with a similar name. The log file lists the number of genes missing from each taxon.
    
 Step 15: Run IQTree on the outputnamehere.dat file. Make sure to rename your outputnamehere.dat.log file before doing this, as the shell will create a file with the same name. You may need to qsub this to a 256G or a higher RAM node, as it needs a lot of memory to be able to run. Step 15: Run IQTree on the outputnamehere.dat file. Make sure to rename your outputnamehere.dat.log file before doing this, as the shell will create a file with the same name. You may need to qsub this to a 256G or a higher RAM node, as it needs a lot of memory to be able to run.
 Here is an example of a shell script for this: Here is an example of a shell script for this:
 +<code>
 #!/bin/sh #!/bin/sh
- 
 #$ -S /bin/sh #$ -S /bin/sh
- 
 . /etc/profile . /etc/profile
- 
 #$ -cwd #$ -cwd
- 
 #$ #$
- 
 #$ -pe threaded 4 #$ -pe threaded 4
  
 iqtree-omp -bb 1000 -wbt -m LG4X -s outputnamehere.dat -nt 4 iqtree-omp -bb 1000 -wbt -m LG4X -s outputnamehere.dat -nt 4
-  +</code> 
-Step 16: Play around with your final tree! Reroot the tree - the common practice is to choose a point between the amorphea and the rest of the clades. Rename the short species names to their long format (Tommy wrote a script called Rename.py. It is located in /scratch2/sarahshah/START.24.05.2017ss/fas_only_badremoved/taxatrimmedfas/bmge), export the tree as a PDF and then edit it using Adobe Illustrator.+Step 16: Play around with your final tree! Reroot the tree - the common practice is to choose a point between the amorphea and the rest of the clades. Rename the short species names to their long format (Tommy wrote a script called RenameTree.py. It is located in /home/sarahshah/AdditionalScripts/Multigene_Phylogeny), export the tree as a PDF and then edit it using Adobe Illustrator.
  
 ---- ----
Line 222: Line 220:
    
 Bordor.sh Bordor.sh
 +
 Author: Yana Eglit Author: Yana Eglit
 +
 This submits the AddPipeline.py. This submits the AddPipeline.py.
    
Line 237: Line 237:
    
 MakeShell.py MakeShell.py
 +
 Author: Tommy Harding Author: Tommy Harding
 +
 Makes shell scripts for tree-making. Makes shell scripts for tree-making.
    
Line 246: Line 248:
    
 cpu_limit.sh cpu_limit.sh
 +
 Author: Tommy Harding & Gordon Lax Author: Tommy Harding & Gordon Lax
 +
 Uses Laura’s Perl script below to control how each tree shell scripts from a list of scripts are submitted to perun. Uses Laura’s Perl script below to control how each tree shell scripts from a list of scripts are submitted to perun.
    
 submit_cpulimit2.pl submit_cpulimit2.pl
 +
 Author: Laura Eme Author: Laura Eme
 +
 A Perl script to control the number of CPUs and waiting time for shell scripts to be submitted to perun. A Perl script to control the number of CPUs and waiting time for shell scripts to be submitted to perun.
 +
 +
 +----
 +
 +**Adding more hits to the pipeline:**
 +
 +This pipeline only takes in the best hit from the local blastp step against the //A. thaliana// genes. If you're interested in looking at the positions of the other hits (the blastp step lists out 10 hits each gene), I have an edited version of the pipeline in /home/sarahshah/AdditionalScripts/Multigene_Phylogeny. The new scripts are: AddPipeline3.1ss.py, addhits3.1ss.py (currently it takes up to 5 hits, you can go in and change that), and cleanhits.sh.
  
multi-gene_phylogeny_pipeline.1497871418.txt.gz · Last modified: by cgeb2001