問題背景

大致是這樣，假如你用一堆注釋軟件然后用不同物種得到了一堆不同的gff属划，那這些gff肯定是有重疊的恬叹，對(duì)于注釋信息gff來說只需要保留gene、exon同眯、cds有這些features的區(qū)段绽昼，并且去除有重疊的區(qū)段僅僅保留一個(gè)，這里不考慮打分長(zhǎng)度什么的（不同軟件不同物種也不好比較分?jǐn)?shù)须蜗，當(dāng)然有些軟件可以根據(jù)權(quán)重來合并注釋結(jié)果硅确，比如RNA注釋權(quán)重較高然后同源注釋次之從頭預(yù)測(cè)的最低），于是腳本就如下（其實(shí)用python/R也很簡(jiǎn)單明肮，但是考慮到性能方面還是用zsh/以awk為主體會(huì)快很多）菱农，具體思路就是以features為gene的行作為標(biāo)識(shí)符（因?yàn)檫@個(gè)gene feature下面所有的記錄都是針對(duì)這個(gè)基因的，所以相當(dāng)于分隔記錄符一樣）柿估，用sort排序得到對(duì)應(yīng)標(biāo)識(shí)符順序循未，用awk讀進(jìn)去成數(shù)組，生成“標(biāo)識(shí)符：針對(duì)這個(gè)gene的所有記錄”的哈希表秫舌，然后根據(jù)排序的數(shù)組的標(biāo)識(shí)符順序來得到整個(gè)gff的排序只厘。接下來再判斷排序后的gff的"gene features"的首尾區(qū)間是否有重疊，同樣也是用哈希記錄然后最終輸出舅巷。(簡(jiǎn)書這個(gè)縮進(jìn)我屬實(shí)無語(yǔ)懶得改了)

#!/usr/bin/zsh
# -*- coding: utf-8 -*-
### ------------------------------------
### merge gff results generated by annotation pipeline
### ------------------------------------

# Get gffdir from input
while {getopts d: arg} {
        case $arg {
                (d)
                gffdir=$OPTARG
                d=$arg
                # Test if the gffdir exists.
                if [[ -d $gffdir ]] {
                        echo "Your gffdir is $gffdir"
                } else {
                           echo "The directory that you specified does not exist, please specify the correct path."
                           exit
                }
                ;;
                (?)
                echo "Wrong option!"
                ;;
        }
}

#  Test if the -d option is provided.
if [[ -z $d ]] {
        echo "You must use -d to specify the directory that contains all the gff files that you want to merge."
        exit
}

# Merge all gff to a big gff
if [[ -f $gffdir/merged.gff ]] {
        rm -rf $gffdir/merged.gff
}
if [[ -f $gffdir/sortedmerged.gff ]] {
        rm -rf $gffdir/sortedmerged.gff
}
if [[ -f $gffdir/filtermerged.gff ]] {
        rm -rf $gffdir/filtermerged.gff
}
ls $gffdir/*gff | while read gff
do
        print "Proccessing $gff..."
        grep -v "^#" $gff >> $gffdir/merged.gff
done && print "Successfully merged gff files!"

# sort gff
## Get a sorted id array that contains only the rows whose feature is gene
id_sorted=`gawk '$3=="gene"{print $0}' $gffdir/merged.gff |
        sort -t $'\t' -k 1,1 -k 7,7 -k 4n,4 -k 5n,5`
## Use awk and srotedIDarray to get a sorted gff file
gawk -v arr=$id_sorted '
BEGIN{
split(arr,id_sorted,"\n")
RS="\n"
FS="\t"
OFS="\t"
ORS="\n"
}
$3=="gene"{
        recs[id]=id"\n"lines
        id=$0
        lines=""
}
$3!="gene"{
        lines=lines"\n"$0
}
END{
recs[id]=id+"\n"+lines
for (id in id_sorted){
        # print recs[id]
        print recs[id_sorted[id]]
}
}
' $gffdir/merged.gff | sed '/^\s*$/d' > $gffdir/sortedmerged.gff


# Exclude overlap records of sorted gff
gawk 'BEGIN{
RS="\n"
FS="\t"
flag=1
}
$3=="gene"{
if(!start[$1","$7]){
        start[$1","$7]=$4
        end[$1","$7]=$5
        flag=1
        print $0
}
else{
        s=start[$1","$7]
        e=end[$1","$7]
        if((s<=$4&&$4<=e)||(s<=$5&&$5<=$e)){
                flag=0
        }
        else{
                start[$1","$7]=$4
                end[$1","$7]=$5
                flag=1
                print $0
        }
}
}
$3!="gene"&&flag==1{print $0}
' $gffdir/sortedmerged.gff > $gffdir/filtermerged.gff