对于单拷贝直系同源基因的比对结果，常规操作是用Gblocks进行处理，但是Gblocks会删的比较彻底，以至于好多基因的文件夹都是空的，所以这次尝试了用TrimAI来进行删除。
1
2
3
4
| # trimal提取保守序列
for i in *.mafft;do
trimal -in $i -out $i.trimal -automated1;
done
|
然后就又遇到一个新的问题，那就是trimal删除以后那些全是gaps的序列会被直接删掉，在序列串联的时候也不能用了，所以这个需要写一个code把缺失的序列补齐。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
| # 经典allspecies.txt囊括所有物种名称
allSpecies=$(cat allspecies.txt)
# 所有的序列排成一列
for i in *.trimal;do
seqkit seq $i -w 0 > $i.seqkit;
done
# 然后就开始补齐了
for i in *.seqkit;do
grep ">" $i > $i.al;
for j in $allSpecies;do
if grep -w "$j" $i.al;
then
echo "yes $j";
else
echo "$j" >> ${i};
len=$(cat ${i} | wc -L); echo $(eval "printf "%0.s-" {1..$len}") >> ${i};
fi
done
done
|
但是这行code有个bug，就是我以最长行的列数作为需要补齐的gaps数量，偶尔会出现有些序列还没有名称长，就导致了补齐的gaps会偏多，所以后面还需要调整成以第二行的列数作为gaps的补齐数量。
暂时用不上，等下次用的时候再进行修改了。
