Phylogeny analysis

#extract the coding region of genome
#./ :including all the coding exon region sequences, like atpA.fasta, atpB.fas .. ycf4.fas.
#fasta :the graft sequences of organell genomes (chloroplast or mitochondria) with only one copy of IR regions and single-line fasta filetype, like Cape_cp.single.fasta .. Toco_cp.single.fasta.
#blastdb: the database of graft sequences of organell genomes bulit by blast+.
#1.pl : a perl script to extract the conding region from graft genome according to the local information.
#1.sh : this file.
#result: the ouput files.
#date : Sat Aug 18 08:26:09 CST 2018
mkir blastdb
cd blastdb
cp -s ../fasta/* ./
for n in *fasta ; do makeblastdb -in $n -dbtype nucl & done
cd ..
rm out* tmp* local* result/*
for p in *fas; do for q in `ls fasta` ; do blastn -query $p -db blastdb/$q -outfmt 6 -out out.${p}.${q}; cat out.${p}.${q} | sort -k 7 -n | awk '{a[NR]=$8; if(a[NR]>=a[NR-1]){print }}' | cut -f 2,9,10 | awk '{if($2<$3){print $1,$2,$3,">"}else{print $1,$3,$2,"<"}}' > local.${p}.${q} ; cat local.${p}.${q} | grep ${q%%.*} > tmp.${p}.${q} ; perl 1.pl fasta/${q%%.*}* tmp.${p}.${q} | awk '{if(/^>/){gsub(">","",$1); m=$1;}else{a[m]=a[m]$1}}END{for(k in a){print k,length(a[k]),a[k]}}' > result/${q/.fasta/}.${p} ; done ; done 
rm out* tmp* local*
cd result
for n in `ls ../*fas | cut -c 4- ` ; do cat *${n} > all.${n/fas/tab} ; done

#alignment and phylogenetic analysis
for n in all*tab ; do cat $n | cut -f 1,3 -d " " | sed 's/^/>/' | sed 's/ /\n/' > ${n/tab/fasta} ; done
rm all*tab
for n in *fasta ; do linsi $n > ${n}.mafft & done
ls *mafft | xargs -I {} -P 10 Gblocks {} -t=d -b5=h -b4=10 -e=.gb

#using SequenceMatrix-Windows-1.7.8 and Mesquite transfer into phy and nex filetype
java -jar ~/software/jmodeltest2-master/jmodeltest-2.1.10/jModelTest.jar -d cp.gene.nex -i -f -g 4 -BIC -AIC -AICc -DT -v -a -w -tr 20 > log.cp.gene
java -jar ~/software/jmodeltest2-master/jmodeltest-2.1.10/jModelTest.jar -d mt.gene.nex -i -f -g 4 -BIC -AIC -AICc -DT -v -a -w -tr 10 > log.mt.gene
#raxml
raxmlHPC-PTHREADS-SSE3 -s mt.gene.phy -n mt.gene -m GTRGAMMAI -T 20 -N 1000 -p 12345 -x 12345 -f a > log.raxml.mt.gene
raxmlHPC-PTHREADS-SSE3 -s cp.gene.phy -n cp.gene -m GTRGAMMAI -T 20 -N 1000 -p 12345 -x 12345 -f a > log.raxml.cp.gene

#align the graft genome using mugsy
mugsy -p cp_genome --directory ~/5.juhua.genome.project/00.1.download.data/02.gene.blast.tab/01.tree/05.mugsy/cp.fasta.files/ -c 100 *fasta
mugsy -p mt_genome --directory ~/5.juhua.genome.project/00.1.download.data/02.gene.blast.tab/01.tree/05.mugsy/mt.fasta.files/ -c 100 *fasta
#alignment and build phylogenies of all mitochondrial gneomes
cat mt_genome.maf |sed '/^a/s/ /_/g' | sed -e 's/score=//' -e 's/label=//' -e 's/mult=/N/' | grep -v ^## | awk '{if(/^a/){m=$1}else{print $0 > "result/"m}}'
cd result
for n in a_*N29 ; do cat $n | sed '/^$/d' | awk '{print ">"$2"\n"$7}' > ${n/a_/b_}.fasta ; done 
rm a_*
for n in *fasta ; do linsi $n > ${n}.mafft & done
ls *mafft | xargs -I {} -P 10 Gblocks {} -t=d -b5=h -b4=10 -e=.gb

#alignment and build phylogenies of all chloroplast gneomes
cat cp_genome.maf |sed '/^a/s/ /_/g' | sed -e 's/score=//' -e 's/label=//' -e 's/mult=/N/' | grep -v ^## | awk '{if(/^a/){m=$1}else{print $0 > "result/"m}}'
cd result
for n in a_*N29; do cat $n | sed '/^$/d' | awk '{print ">"$2"\n"$7}' > ${n/a_/b_}.fasta ; done 
rm a_*
for n in *fasta ; do linsi $n > ${n}.mafft & done
ls *mafft | xargs -I {} -P 10 Gblocks {} -t=d -b5=h -b4=10 -e=.gb
ls *mafft | xargs -I {} -P 10 Gblocks {} -t=d -b5=n -b4=10 -e=.gb


#using SequenceMatrix-Windows-1.7.8 and Mesquite transfer into phy and nex filetype
java -jar ~/software/jmodeltest2-master/jmodeltest-2.1.10/jModelTest.jar -d cp.genome.nex -i -f -g 4 -BIC -AIC -AICc -DT -v -a -w -tr 20 > log.cp.genome
java -jar ~/software/jmodeltest2-master/jmodeltest-2.1.10/jModelTest.jar -d mt.genome.nex -i -f -g 4 -BIC -AIC -AICc -DT -v -a -w -tr 10 > log.mt.genome
raxmlHPC-PTHREADS-SSE3 -s mt.genome.phy -n mt.genome -m GTRGAMMAI -T 20 -N 1000 -p 12345 -x 12345 -f a > log.raxml.mt.genome
raxmlHPC-PTHREADS-SSE3 -s cp.genome.phy -n cp.genome -m GTRGAMMAI -T 20 -N 1000 -p 12345 -x 12345 -f a > log.raxml.cp.genome