基因组注释流程

一、使用Regtag将contig挂到染色体上

###这里没有进行补洞，补洞代码 ragtag.py patch  补完后染色体大小会变得一样大，会丢失很多变异信息
cat id1 |while read id
do
ragtag.py scaffold /public/home/fengting/database/rice/9311/9311.fa \ 
/public/home/fengting/task/5.12ragtag/dataNIP/${id}.fasta \
-o data9311/${id}.fa
done

二、使用Repeatmasker进行基因组数据屏蔽：

cat /id|while read id 
do
RepeatMasker -q -nolow -no_is -norna  -e ncbi   \ 
-species rice -div 40 -cutoff 255 -pa 88 -s -gff \ 
-dir ../datazs97/${id}.fa
done

三、基因预测

August使用未屏蔽基因组

• August从头预测：

cat id|while read id
do augustus --species=rice --gff3=on ${id}.fa > aug/${id}.gff
done

• genemark 从头预测

~/soft/gmes_linux_64/gmes_petap.pl --ES --sequence /public/home/fengting/task/5.12anno/masked/${id}.fa.masked --cores 50 

##格式转换
gffread $id.gtf -o gmk/$id.gff

• gth同源预测

cat ../id1|while read id 
do
gth -genomic ../masked/${id}.fa.masked -protein /public/home/fengting/database/rice/ri/Oryza_sativa.IRGSP-1.0.pep.all.fa -intermediate -gff3out > out/${id}.gff
done

使用python脚本进行格式转换

#!/lustre/home/guohan_lab/local/python-3.6/bin/python3
#liyumei
#./change-name_lym.py proteinprediction.gff proteinprediction.gff proteinprediction_r.gff

import sys
import re
fout = open(sys.argv[3],'w')

ref_dict={}
with open(sys.argv[1]) as gff:
        for line in gff:
                line_s = line.strip().split('\t')
                if 'gene' == line_s[2]:
                        ref_g = re.split('=',line_s[8])
                        ref_gene = ref_g[1]
                        ref_dict[ref_gene]=[]
                else:
                        pos = line_s[3]
                        ref_dict[ref_gene].append(pos)

with open(sys.argv[2]) as gff_r:
        for eachline in gff_r:
                i = eachline.strip().split('\t')
                info = re.split('=',i[8])
                name = info[1]
                ref_set = []
                for n in range(0,8):
                        ref_set.append(i[n])
                ref_list = '\t'.join(ref_set)
                ref_set1 = ['CDS' if x == 'exon' else x for x in ref_set]
                ref_list1 = '\t'.join(ref_set1)
                if 'gene' == i[2]:
                        fout.write(eachline)
                elif 'exon' == i[2]:
                        if name in ref_dict:
                                vs = ref_dict[name]
                                len_vs = len(vs)
                                for a in range(0,len_vs):
                                        if vs[a] == i[3]:
                                                b = vs.index(vs[a])
                                                fout.write('%s\tID=%s.t%d;Parent=%s\n'%(ref_list,name,b+1,name))
                                                fout.write('%s\tID=%s.c%d;Parent=%s\n'%(ref_list1,name,b+1,name))
fout.close()

cat /public/home/fengting/task/222anno/id1 |while read id
do
#/public/home/fengting/tools/EVidenceModeler-1.1.1/EvmUtils/misc/genomeThreader_to_evm_gff3.pl ${id}.gff > gth/evm.${id}.gff

grep -v '^#' ${id}.gff |grep -v 'prime_cis_splice_site'|awk -F ';' '{print$1}' > gth/gth.${id}.gff
done

cat /public/home/fengting/task/222anno/id1 |while read id
do
python ../change-gff_lum.py gth.${id}.gff gth.${id}.gff gth/gth.${id}.gff
done

• TransDecoder转录本预测

cat id|while read id
do
hisat2-build /public/home/fengting/task/5.12anno/masked/${id}.fa.masked index/${id}
hisat2 --dta -p 20 -x index/${id} -1 /public/home/lianglunping/work/HXM/RNA/input/HA798-normal-1A_1.fq.gz -2 /public/home/lianglunping/work/HXM/RNA/input/HA798-normal-1A_2.fq.gz -S sam/${id}.1.sam
hisat2 --dta -p 20 -x index/${id} -1 /public/home/lianglunping/work/HXM/RNA/input/HA798-roll-1A_1.fq.gz -2 /public/home/lianglunping/work/HXM/RNA/input/HA798-roll-1A_2.fq.gz -S sam/${id}.2.sam
hisat2 --dta -p 20 -x index/${id} -1 /public/home/lianglunping/work/HXM/RNA/input/NPB-1A_1.fq.gz -2 /public/home/lianglunping/work/HXM/RNA/input/NPB-1A_2.fq.gz  -S sam/${id}.3.sam
hisat2 --dta -p 20 -x index/${id} -1 /public/home/lianglunping/work/HXM/RNA/input/TR-grass-1A_1.fq.gz -2 /public/home/lianglunping/work/HXM/RNA/input/TR-grass-1A_2.fq.gz  -S sam/${id}.4.sam
hisat2 --dta -p 20 -x index/${id} -1 /public/home/lianglunping/work/HXM/RNA/input/TR-normal-1A_1.fq.gz -2 /public/home/lianglunping/work/HXM/RNA/input/TR-normal-2A_1.fq.gz  -S sam/${id}.5.sam
samtools view -bS sam/$id.1.sam -o bam/$id.1.bam
samtools view -bS sam/$id.2.sam -o bam/$id.2.bam
samtools view -bS sam/$id.3.sam -o bam/$id.3.bam
samtools view -bS sam/$id.4.sam -o bam/$id.4.bam
samtools view -bS sam/$id.5.sam -o bam/$id.5.bam
samtools merge -@ 88   bam/$id.merge    bam/$id.1.bam bam/$id.2.bam bam/$id.3.bam bam/$id.4.bam bam/$id.5.bam 

samtools sort -@ 88 bam/$id.merge  -o bam/$id.s.bam

stringtie -p 10 -o sbam/$id.gtf bam/$id.s.bam

/public/home/fengting/tools/TransDecoder-v5.5.0/util/gtf_genome_to_cdna_fasta.pl sbam/$id.gtf /public/home/fengting/task/5.12anno/masked/${id}.fa.masked > tra/transcripts.${id}.fasta
/public/home/fengting/tools/TransDecoder-v5.5.0/util/gtf_to_alignment_gff3.pl sbam/$id.gtf > tra/transcripts.${id}.gff3
TransDecoder.LongOrfs -t tra/transcripts.${id}.fasta
TransDecoder.Predict -t tra/transcripts.${id}.fasta
/public/home/fengting/tools/TransDecoder-v5.5.0/util/cdna_alignment_orf_to_genome_orf.pl \
     transcripts.${id}.fasta.transdecoder.gff3 \
     tra/transcripts.${id}.gff3 \
     tra/transcripts.${id}.fasta > out/trs.${id}.gff3
done

四、evm整合

权重文件
#####weights.txt  
PROTEIN gth 4
TRANSCRIPT  transdecoder    8
ABINITIO_PREDICTION AUGUSTUS    1
ABINITIO_PREDICTION GeneMark.hmm3   1

/public/home/fengting/tools/EVidenceModeler-1.1.1/EvmUtils/partition_EVM_inputs.pl \
    --genome data/H7L33.fa \
    --gene_predictions gpre/H7L33.gff  \
    --protein_alignments gth/gth.H7L33.gff  \
    --transcript_alignments trs/trs.H7L33.gff3 \
    --segmentSize 100000 --overlapSize 10000 --partition_listing out/H7L33/partitions_list.out

/public/home/fengting/tools/EVidenceModeler-1.1.1/EvmUtils/write_EVM_commands.pl \
    --genome data/H7L33.fa \
    --weights `pwd`/weights.txt \
        --gene_predictions gpre/H7L33.gff \
        --protein_alignments gth/gth.H7L33.gff  \
        --transcript_alignments trs/trs.H7L33.gff3 \
        --output_file_name evm.out --partitions out/H7L33/partitions_list.out >  out/H7L33/commands.list
#parallel --jobs 10 < commands.list

sh out/H7L33/commands.list

~/tools/EVidenceModeler-1.1.1/EvmUtils/recombine_EVM_partial_outputs.pl --partitions out/H7L33/partitions_list.out --output_file_name out/H7L33/evm.out

/public/home/fengting/tools/EVidenceModeler-1.1.1/EvmUtils/convert_EVM_outputs_to_GFF3.pl  --partitions out/H7L33/partitions_list.out --output evm.out  \
    --genome data/H7L33.fa
find . -regex ".*evm.out.gff3" -exec cat {} \; | /public/home/fengting/miniconda3/envs/GS/bin/bedtools sort -i - > out/H7L33/H7L33.EVM.all.gff

注释结果

注释结果存在部分假阳性和假阴性的问题

五、ppacp进行PAN注释

/public/home/fengting/demo/ppscp/ppsPCP/bin/make_pan.pl \
--ref /public/home/fengting/task/5.12anno/datazs97/H7L1.fa \
--ref_anno /public/home/fengting/task/3.26pan/new/gmk/H7L1.gff \
--query /public/home/fengting/task/5.12anno/datazs97/H7L26.fa  \
/public/home/fengting/task/5.12anno/datazs97/H7L27.fa \
/public/home/fengting/task/5.12anno/datazs97/H7L28.fa \
/public/home/fengting/task/5.12anno/datazs97/H7L29.fa \
/public/home/fengting/task/5.12anno/datazs97/H7L30.fa \
/public/home/fengting/task/5.12anno/datazs97/H7L31.fa \
/public/home/fengting/task/5.12anno/datazs97/H7L32.fa \
/public/home/fengting/task/5.12anno/datazs97/H7L33.fa \
--query_anno /public/home/fengting/task/3.26pan/new/gmk/H7L26.gff \
/public/home/fengting/task/3.26pan/new/gmk/H7L27.gff \
/public/home/fengting/task/3.26pan/new/gmk/H7L28.gff \
/public/home/fengting/task/3.26pan/new/gmk/H7L29.gff \
/public/home/fengting/task/3.26pan/new/gmk/H7L30.gff \
/public/home/fengting/task/3.26pan/new/gmk/H7L31.gff \
/public/home/fengting/task/3.26pan/new/gmk/H7L32.gff \
/public/home/fengting/task/3.26pan/new/gmk/H7L33.gff \
--tmp rice_tmp2/
--thread 88