交流是最快的成長方式，最近就自己被坑了。哎不說了寫點有用的吧。分享個基因家族分析比較完整的pipelin

#!/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
cat all_to_fungi | awk '$3>60 && $11<1e-10'>0-fungi10
for a in $(cat falist); do grep $a 0-moregenome > more_genome/$a.txt; done
###list############
mkdir list
for a in $(ls *.txt)
do
    cat $a |cut -f2|sort -u >list/${a%.*}.list
done

#############
###hmm
cd /all_hmm
for a in *
do
    hmmsearch --cpu 10 --domtblout ${a%.*}.out $a genome.faa
done
###list
for a in $(ls *.out); do grep -v "#" $a|awk '($7 + 0) < 1E-10'|cut -f1 -d  " "|sort -u > /data/shiyan/data/2-hmm/list/${a%.*}.list; done
###comm
comm -12 *.list *.list > common.list
###fa###########
mkdir all_fa
for a in $(ls all_list/)
do
    less genome.faa | seqkit grep -f all_list/$a > ${a%.*}.fa
done
###


###matff###########
cd /home/*/all_fa
mkdir mafft
for a in $(ls *.fa)
do
    mafft --auto --thread 16 --inputorder --anysymbol $a > ${a%.*}.faa
done
###保守結(jié)構(gòu)域
Gblocks proteins.fasta -t=p
###iqtree
iqtree -s AglAG12.faa -m MFP -bb 1000 -bnni -nt AUTO -ntmax 10