交流是最快的成长方式，最近就自己被坑了。哎不说了写点有用的吧。分享个基因家族分析比较完整的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
###保守结构域
Gblocks proteins.fasta -t=p
###iqtree
iqtree -s AglAG12.faa -m MFP -bb 1000 -bnni -nt AUTO -ntmax 10
