linux脚本 | 线粒体组基因提取合并

介绍

基于线粒体基因组进行构树时，一般仅以其中的编码基因(需要的话也可以提取rRNA)完成。之就需要对每个物种的线粒体基因组的每个基因序列进行提取，合并。无论是自己拼装的还是NCBI发表的线粒体全序，想要快速的统一格式，最简单的办法还是直接去MITOS(http://mitos2.bioinf.uni-leipzig.de/index.py)进行注释，这样所有基因的名称即可统一，且规避掉起始位置不统一、序列需要单个提取这样的麻烦。
MITOS注释完后，我们直接用MITOS提供的fas序列（MITOS提供的核苷酸序列，并非自己的原始序列），对多物种序列进行整理即可。

文件准备

• 核酸序列的fas文件XXX.fas，XXX无需任何含义
• 需要提取的基因名列表，单列，储存为name.list.txt
• 脚本，run_MitoGene.sh

#  run_MitoGene.sh
for fn in *.fas
do
        # 提取物种名
        name=`cat ${fn} | grep ">" | awk -F, 'BEGIN{FS=";"}{print $1}' | awk -F, 'BEGIN{FS=">"}{print $2}' | uniq`
        # 简化基因名，仅剩下物种名和基因名
        cat ${fn} | awk -F, 'BEGIN{FS="("}{print $1}' | awk  -F, 'BEGIN{FS=";"}{print $1,$4}'| awk '{print $1,$2}' > ${name}_sn.fas
done
#去掉序列中的物种名
for fn in *_sn.fas
do
        echo $fn
        name=${fn%_*}
        echo $name
        sed 's/>'"${name}"' />/1g' < ${fn} > ${name}_rn.fas
done
#将需要的序列通过name list 提取
cat name.list.txt | while read line
do
        #seqtk subseq 提取的基因名需要一个单独文件
        echo ${line} > this.name.txt
        #按基因提取，从每个物种中提取相同基因并合并
        for fn in *_rn.fas
        do
                sp_name=${fn%_*}
                seqtk subseq $fn this.name.txt > gene.tmp
                #提取到tmp文件中，重新加上物种名而后合并
                sed 's/>/>'"${sp_name}"'|/g' < gene.tmp >> ${line}_tmp.fas
        done
        #去掉所有空格
        sed 's/ //g' < ${line}_tmp.fas > ${line}_GeneSet.fas
        #删除基因集tmp文件
        rm ${line}_tmp.fas
done
rm this.name.txt
rm gene.tmp
#
#去掉基因集里的基因名，仅保留统一的物种名
for fn in *_GeneSet.fas
do
        name=${fn%_*}
        sed 's/|/]>/g' ${fn} | awk -F, 'BEGIN{FS="]"}{print $2,$1}' | awk '{print $2,$1}' | awk '{print $1}' > ${name}_Del.fas
done
#
#整理所有文件
mkdir simple_name
mkdir single_name
mkdir Gene_set
mkdir GS_Del
mv *_sn.fas ./simple_name
mv *_rn.fas ./single_name
mv *_GeneSet.fas ./Gene_set
mv *_Del.fas ./GS_Del
mkdir raw_data
mv *.fas ./raw_data

结果

原始数据被打包仅raw_data中
最终的结果被放在GS_Del文件中，文件按基因名称命名，每个基因的序列集中仅保留了各物种的物种名，从而方便多基因直接联合。
后序可以仍在linux中进行比对、删减等处理。物种不多的话，直接导入PhyloSuite即可。