phylogenetic建树及自动美化

好久没更新了，来个福利
step1 tree

#!/usr/bin/bash
###blast
makeblastdb -in genome.faa -out genome -dbtype prot
blastp -query all.fa -db genome -out all_to_genome -outfmt 6 -evalue 1e-5 -num_threads 12
#########################################################################################

###every gene blatp out
cat all_to_genome | awk '$3>60 && $11<1e-10'>0-all_to_genome
mkdir genome_out
for a in $(cat falist)
do 
grep $a all_to_genome > genome_out/$a.txt
done
#########################################################################################

###get every gene blastp out list
cd genome_out
mkdir list
for a in $(ls *.txt)
do
    cat $a |cut -f2|sort -u >list/${a%.*}.list
done

#####
###hmm
#cd hmm
#for a in *
#do
#   hmmsearch --cpu 10 --domtblout hmm_out/${a%.*}.out $a /data/genome.faa
#done
#for a in $(ls *.out); do grep -v "#" $a|awk '($7 + 0) < 1E-10'|cut -f1 -d  " "|sort -u > list/${a%.*}.list; done
#comm -12 blastp_result_id.list final.NBS.list > common.list

#########################################################################################

###get fa from list
cd list
mkdir fa
for a in $(ls *.list)
do
    less genome.faa | seqkit grep -f $a > fa/${a%.*}.fa
done
#########################################################################################

###matff for multiple sequence alignment
cd fa
mkdir mafft
for a in $(ls *.fa)
do
    mafft --auto --thread 16 --inputorder --anysymbol $a > mafft/${a%.*}.faa
done
##########################################################################################

###Gblocks for Conserved domains
cd matff
for a in $(ls *.faa)
do
Gblocks $a -t=p -b4=2 -b5=a
done
##########################################################################################

###iqtree to tree
for a in $(ls -Sr *.fa-gb)
do
echo "iqtree -s $a -mset LG,JTT -mrate G,I,G+I -bb 1000 -bnni -nt AUTO -ntmax 4" >>$a.sh
done

for a in *.bash
do
echo "qs -m1 -p1 $a" >> qs.sh
do
bash qs.sh
##########################################################################################

step2 get list

#!/usr/bin/bash
#####
for name in $(ls ../*.faa)
do
for a in $(grep ">" $name |awk -F ' ' '{print $1;}')
do 
    echo "${a#>*}">>${name%.*}.list
done
done

#####
for a in $(ls *.list); do cat $a |sort >${a%.*}.sort; done
rm *.list

#####
for list in $(ls *.sort)
do
    for a in $(cat RNS.txt)
    do
        grep $a $list >> ${list%.*}.A
    done
done


#####
for list in $(ls *.sort)
do
       for a in $(cat AMS.txt)
       do
               grep $a $list >> ${list%.*}.B
       done
done

#####
for list in $(ls *.sort)
do
    cat ${list%.*}.A ${list%.*}.B |sort>${list%.*}.A_B
    comm -13 ${list%.*}.A_B $list >${list%.*}.C && rm ${list%.*}.A_B
done

#####
for list in $(ls *.sort)
do
    for a in $(cat ${list%.*}.A); do echo -e $a"\tA">>${list%.*}.list; done
    for a in $(cat ${list%.*}.B); do echo -e $a"\tB">>${list%.*}.list; done
    for a in $(cat ${list%.*}.C); do echo -e $a"\tC">>${list%.*}.list; done
done

#####
rm *.A *.B *.C *.sort

step3 ggtree

###R ggtree pack
#!/usr/bin/bash
setwd("D:/Desktop/tree")
library("ape")
library("Biostrings")
library("ggplot2")
library("ggtree")
library("treeio") 
#tree=read.tree("PvSIN1_Medtr4g133660.1.faa.treefile") 
tree <- read.newick("LjSST1_Medtr6g086170.1.faa.treefile",node.label = "support")
group_file <- read.table("LjSST1_Medtr6g086170.1.list",header = T,row.names = 1)
groupInfo <- split(row.names(group_file), group_file$Group)
# 将分组信息添加到树中
tree <- groupOTU(tree, groupInfo)
#p<-ggtree(tree)
# 绘制进化树
tregraph=ggtree(tree, layout="rectangular", size=0.3, aes(color=group)) +
  # 树体
  geom_tiplab(size=0.8, aes(color=group), hjust=-0.05) +
  # 枝名
  geom_text2(aes(subset=!isTip,label=support, hjust=-0.5),size=1)+
  # bootstrap
  geom_tippoint(size=0.1, aes(color=group)) +
  # point
  geom_nodepoint(aes(color=group), alpha=1/4, size=0.1) +
  # 末节点
  theme_tree2() +
  # x轴标尺
  xlim(NA, 8)
# x轴宽度

ggsave(filename = "LjSST1_Medtr6g086170.1.pdf",
       plot= tregraph,height= 20,width= 20,units= 'in', dpi= 300)

#pdf("Glyma.YUCCA2.pdf")
#tregraph
#dev.off()
#b=tree[["node.label"]]
#write.csv(b, file ="A.csv", row.names =FALSE)

喜欢的来个赞，感谢