搜索
查看: 1009|回复: 1

非模式生物blast注释笔记

[复制链接]

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
QQ
发表于 2019-10-10 09:09:50 | 显示全部楼层 |阅读模式
本帖最后由 惠真-市二医院 于 2019-10-10 09:14 编辑

Blast不是diamond最近做单菌的转录组分析,研究了下无参转录组,因此就碰到了不少问题,这里笔记一下
1. 组装后的scaffolds,不要直接使用软件预测,因为测序错误会引入转录终止子,这样,就会把长蛋白序列折短。直接对组装好的scaffolds,blast预测,参考数据库呢,一定要越全越好,ncbi的nr,我这里用的是uniprot的Reviewed (Swiss-Prot)/Unreviewed (TrEMBL)。那问题就是,blast速度超级慢,我这里的办法就是,先通过taxomy将数据库数据进行拆分,然后再blast。2. blast注释过滤,我的策略就是blast时可先简单设置写下e值,然后写脚本对结果过滤,blast各个阈值见图。
挺多啊,哈哈,我采用的策略是,query至少要覆盖target的60%,然后选择最大的bitscore。脚本如下,估计大家也不想看,可脚本着实花了我好久才搞好其中逻辑顺序:
[Shell] 纯文本查看 复制代码
##diamond blastx
##diamond blastx --query-gencode 11 -d ../sprot.dmnd -q trinity_out_dir/kp_21_trinity_unique.fasta -f 6 qseqid qlen qstart qend sseqid slen sstart send length pident mismatch gapopen bitscore qcovhsp evalue btop qseq > kp_21_trinity_prokka.results
##diamond blastx qseq为query fasta sequence,if(start < end):substr($seq,start,end-start+1);if(start >end):substr($seq,end,start-end+1)

##blastx
##blastx -db ../Bacteria/sprot -evalue 1e-5 -query_gencode 11 -query trinity_out_dir/kp_21_trinity_unique.fasta -outfmt "6 qseqid qlen qstart qend sseqid slen sstart send length pident mismatch gapopen bitscore qcovs evalue btop qseq" > kp_21_trinity_prokka_blast.results

##extact the mapping interval from contig fasta files, if(start < end):substr($seq,start,end-start+1); if(start > end):revcomp(substr($seq,end,start-end+1)) equal to the Protein sequences

##querycov filter > 60 query interval sequence mapping at least 60% of the target sequence

##Filter the search result with the criteria:If the collected interval would intersect with the next interval, the script would selcet the one with the higher bitscore
BEGIN{
    FS="\t"
    OFS="\t"
}
{
    if($1 in Contig){
        Tmp[$1]=0
        for(var in Annot){
             #print "length Annot",length(Annot)
            split(Annot[var],T,"\t")
            if(T[1] == $1){
                S[1]=$3;S[2]=$4;S[3]=T[3];S[4]=T[4]
                #print "before",S[1],S[2],S[3],S[4]
                asort(S)
                if( ((S[1]==$3) && (S[2]==$4)) || ((S[1]==$4) && (S[2]==$3)) || ((S[1]==T[3]) && (S[2]==T[4])) || ((S[1]==T[4]) && (S[2]==T[3])) ){
                    Tmp[$1]=Tmp[$1]+1
                    #print "alfter",S[1],S[2],S[3],S[4]
                    #print Tmp[$1],Contig[$1]
                }else{
                    #print "$13",$13,T[13]
                    if($13 > T[13]){
                        #print Annot[var]
                        print Annot[var] > "Omitted_results_by_filter.txt"
                        delete Annot[var]
                        Annot[$0]=$0
                        #Tmp[$1]=0
                        next
                    }else{
                        #print Annot[var]
                        print $0 > "Omitted_results_by_filter.txt"
                        #Tmp[$1]=0
                        next
                    }
                }
            }
        }
            #print Tmp[$1],Contig[$1],$1
            if(Tmp[$1]==Contig[$1]){
                    Annot[$0]=$0
                    Contig[$1]=Contig[$1]+1
                    #print "length",length(Annot)
                    next
            }

    }else{
        Contig[$1]=Contig[$1]+1
        #print "Contig[$1]",Contig[$1]
        Annot[$0]=$0
    }
}


END{
    for(var in Annot){
        print Annot[var]
    }
}


最后,请坚持使用blast,diamond虽然快但是敏感性差的不是一点半点!

8.png
7.png
6.png
5.png
4.png
3.png
2.png
1.png
苛求远离完美
回复

使用道具 举报

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
QQ
 楼主| 发表于 2019-10-10 09:14:52 | 显示全部楼层
按名称数字顺序看图
苛求远离完美
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-7-13 09:36 , Processed in 0.025757 second(s), 32 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.