From 5b804ddeb65b07df32b4f34c9b80c9b386535807 Mon Sep 17 00:00:00 2001 From: Tomas Fer Date: Mon, 6 Feb 2017 21:11:05 +0100 Subject: [PATCH] v.1.3.3 Add possibility to estimate gene trees using 'by exon' partitioned dataset. A *.part file is prepared when concatenating *.mafft alignments in script 4 and later used in RAxML gene tree building in script 6a. --- HybPhyloMaker4_processpslx.sh | 14 +++++++++--- HybPhyloMaker6a_RAxML_for_selected.sh | 32 ++++++++++++++++++++++----- 2 files changed, 37 insertions(+), 9 deletions(-) diff --git a/HybPhyloMaker4_processpslx.sh b/HybPhyloMaker4_processpslx.sh index d829dde..4a0cb30 100644 --- a/HybPhyloMaker4_processpslx.sh +++ b/HybPhyloMaker4_processpslx.sh @@ -22,7 +22,7 @@ # ******************************************************************************** # * HybPhyloMaker - Pipeline for Hyb-Seq data processing and tree building * # * Script 04 - Process pslx files * -# * v.1.3.2 * +# * v.1.3.3 * # * Tomas Fer, Dept. of Botany, Charles University, Prague, Czech Republic, 2016 * # * tomas.fer@natur.cuni.cz * # * based on Weitemier et al. (2014), Applications in Plant Science 2(9): 1400042* @@ -324,8 +324,10 @@ if [[ $cp =~ "no" ]]; then #Copy script if [[ $location == "1" ]]; then cp -r $source/catfasta2phyml.pl . + cp -r $source/AMAS.py . else cp -r ../$source/catfasta2phyml.pl . + cp -r ../$source/AMAS.py . fi #Modify mafft file names (from, i.e., To_align_Assembly_10372_Contig_1_516.fasta.mafft to To_align_Assembly_10372_*mafft) #(all files starting with "To_align_Assembly_10372_" will be merged) @@ -347,8 +349,14 @@ if [[ $cp =~ "no" ]]; then #transform fasta to phylip format, copy results from scratch to home cat fileForLoop.txt | while read -r a b do - perl catfasta2phyml.pl -f $a > $b.fasta - perl catfasta2phyml.pl $b.fasta > $b.phylip + #concatenate exons from the same locus + python3 AMAS.py concat -i $a -f fasta -d dna -u fasta -t ${b}.fasta >/dev/null + python3 AMAS.py concat -i $a -f fasta -d dna -u phylip -t ${b}.phylip >/dev/null + #modify and rename partitions.txt + sed -i.bak -e 's/^/DNA, /g' -e 's/_To_align//g' partitions.txt + mv partitions.txt ${b}.part + #perl catfasta2phyml.pl -f $a > $b.fasta + #perl catfasta2phyml.pl $b.fasta > $b.phylip if [[ $location == "1" ]]; then cp $b.* $path/$type/70concatenated_exon_alignments else diff --git a/HybPhyloMaker6a_RAxML_for_selected.sh b/HybPhyloMaker6a_RAxML_for_selected.sh index e5b5580..578d6c1 100644 --- a/HybPhyloMaker6a_RAxML_for_selected.sh +++ b/HybPhyloMaker6a_RAxML_for_selected.sh @@ -20,7 +20,7 @@ # ******************************************************************************** # * HybPhyloMaker - Pipeline for Hyb-Seq data processing and tree building * # * Script 06a - RAxML gene tree building * -# * v.1.3.2 * +# * v.1.3.3 * # * Tomas Fer, Dept. of Botany, Charles University, Prague, Czech Republic, 2015 * # * tomas.fer@natur.cuni.cz * # ******************************************************************************** @@ -164,11 +164,14 @@ if [[ $location == "1" || $location == "2" ]]; then echo 'SPECIESPRESENCE='"$SPECIESPRESENCE" >> ${group}.sh echo 'type='"$type" >> ${group}.sh echo 'location='"$location" >> ${group}.sh + echo 'genetreepart='"$genetreepart" >> ${group}.sh + echo 'raxmlpthreads='"$raxmlpthreads" >> ${group}.sh echo 'cp '"$path"'/$type/71selected${MISSINGPERCENT}/'"$group"' .' >> ${group}.sh echo 'cp '"$source"'/catfasta2phyml.pl .' >> ${group}.sh echo 'for i in $(cat '"$group"')' >> ${group}.sh echo 'do' >> ${group}.sh echo ' cp '"$path"'/$type/71selected${MISSINGPERCENT}/deleted_above${MISSINGPERCENT}/${i}_modif${MISSINGPERCENT}.fas .' >> ${group}.sh + echo ' cp $path/$type/70concatenated_exon_alignments/${i}.part .' >> ${group}.sh echo ' #Substitute '"'('"' by '"'_'"' and '"')'"' by nothing ('"'('"' and '"')'"' not allowed in RAxML)' >> ${group}.sh echo ' sed -i '"'s/(/_/g'"' ${i}_modif${MISSINGPERCENT}.fas' >> ${group}.sh echo ' sed -i '"'s/)//g'"' ${i}_modif${MISSINGPERCENT}.fas' >> ${group}.sh @@ -179,13 +182,22 @@ if [[ $location == "1" || $location == "2" ]]; then echo 'done' >> ${group}.sh echo '#Make a list of all phylip files' >> ${group}.sh echo 'ls *.phylip | cut -d"." -f1 > FileForRAxML.txt' >> ${group}.sh - echo 'for file in $(cat FileForRAxML.txt)' >> ${group}.sh - echo 'do' >> ${group}.sh + echo 'for file in $(cat FileForRAxML.txt); do' >> ${group}.sh + echo ' #modify name for partition file (remove '_modif${MISSINGPERCENT}')' >> ${group}.sh + echo ' filepart=$(sed "s/_modif${MISSINGPERCENT}//" <<< $file)' >> ${group}.sh echo ' #RAxML with 100 rapid bootstrap' >> ${group}.sh echo ' if [[ $location == "1" ]]; then' >> ${group}.sh - echo ' raxmlHPC-PTHREADS -T $TORQUE_RESC_TOTAL_PROCS -f a -s $file.phylip -n $file.result -m GTRCAT -p 1234 -x 1234 -N 100' >> ${group}.sh + echo ' if [[ $genetreepart == "yes" ]]; then' >> ${group}.sh + echo ' raxmlHPC-PTHREADS -T $TORQUE_RESC_TOTAL_PROCS -f a -s $file.phylip -q $filepart.part -n $file.result -m GTRCAT -p 1234 -x 1234 -N 100' >> ${group}.sh + echo ' else' >> ${group}.sh + echo ' raxmlHPC-PTHREADS -T $TORQUE_RESC_TOTAL_PROCS -f a -s $file.phylip -n $file.result -m GTRCAT -p 1234 -x 1234 -N 100' >> ${group}.sh + echo ' fi' >> ${group}.sh echo ' elif [[ $location == "2" ]]; then' >> ${group}.sh - echo ' $raxmlpthreads -T $NSLOTS -f a -s $file.phylip -n $file.result -m GTRCAT -p 1234 -x 1234 -N 100' >> ${group}.sh + echo ' if [[ $genetreepart == "yes" ]]; then' >> ${group}.sh + echo ' $raxmlpthreads -T $NSLOTS -f a -s $file.phylip -q $filepart.part -n $file.result -m GTRCAT -p 1234 -x 1234 -N 100' >> ${group}.sh + echo ' else' >> ${group}.sh + echo ' $raxmlpthreads -T $NSLOTS -f a -s $file.phylip -n $file.result -m GTRCAT -p 1234 -x 1234 -N 100' >> ${group}.sh + echo ' fi' >> ${group}.sh echo ' fi' >> ${group}.sh echo ' cp *$file.result '"$path"'/$type/72trees'"${MISSINGPERCENT}_${SPECIESPRESENCE}"'/RAxML' >> ${group}.sh echo 'done' >> ${group}.sh @@ -216,6 +228,7 @@ else for i in $(cat selected_genes_${MISSINGPERCENT}_${SPECIESPRESENCE}.txt) do cp $path/$type/71selected${MISSINGPERCENT}/deleted_above${MISSINGPERCENT}/${i}_modif${MISSINGPERCENT}.fas . + cp $path/$type/70concatenated_exon_alignments/${i}.part . #Substitute '(' by '_' and ')' by nothing ('(' and ')' not allowed in RAxML) sed -i.bak 's/(/_/g' ${i}_modif${MISSINGPERCENT}.fas sed -i.bak 's/)//g' ${i}_modif${MISSINGPERCENT}.fas @@ -232,7 +245,14 @@ else calculating=$((calculating + 1)) echo -e "\nProcessing file: ${file}" >> raxml.log echo -e "Processing file: ${file} ($calculating out of $numbertrees)" - $raxmlpthreads -T $numbcores -f a -s $file.fas -n $file.result -m GTRCAT -p 1234 -x 1234 -N 100 >> raxml.log + #modify name for partition file (remove '_modif${MISSINGPERCENT}' + filepart=$(sed "s/_modif${MISSINGPERCENT}//" <<< $file) + #run RAxML + if [[ $genetreepart == "yes" ]]; then + $raxmlpthreads -T $numbcores -f a -s $file.fas -q $filepart.part -n $file.result -m GTRCAT -p 1234 -x 1234 -N 100 >> raxml.log + else + $raxmlpthreads -T $numbcores -f a -s $file.fas -q -n $file.result -m GTRCAT -p 1234 -x 1234 -N 100 >> raxml.log + fi cp *$file.result $path/$type/72trees${MISSINGPERCENT}_${SPECIESPRESENCE}/RAxML done #Copy raxml.log to home