写个zsh脚本(主要是awk)合并gff,去除gff中的重叠区域

问题背景

大致是这样,假如你用一堆注释软件然后用不同物种得到了一堆不同的gff,那这些gff肯定是有重叠的,对于注释信息gff来说只需要保留gene、exon、cds有这些features的区段,并且去除有重叠的区段仅仅保留一个,这里不考虑打分长度什么的(不同软件不同物种也不好比较分数,当然有些软件可以根据权重来合并注释结果,比如RNA注释权重较高然后同源注释次之从头预测的最低),于是脚本就如下(其实用python/R也很简单,但是考虑到性能方面还是用zsh/以awk为主体会快很多),具体思路就是以features为gene的行作为标识符(因为这个gene feature下面所有的记录都是针对这个基因的,所以相当于分隔记录符一样),用sort排序得到对应标识符顺序,用awk读进去成数组,生成“标识符:针对这个gene的所有记录”的哈希表,然后根据排序的数组的标识符顺序来得到整个gff的排序。接下来再判断排序后的gff的"gene features"的首尾区间是否有重叠,同样也是用哈希记录然后最终输出。(简书这个缩进我属实无语懒得改了)

#!/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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,324评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,303评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,192评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,555评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,569评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,566评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,927评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,583评论 0 257
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,827评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,590评论 2 320
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,669评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,365评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,941评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,928评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,159评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,880评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,399评论 2 342

推荐阅读更多精彩内容