Gene Prediction Commend 01
若使用本方法，请引用此文，谢谢！
JoF | Free Full-Text | Genome Re-Annotation and Transcriptome Analyses of Sanghuangporus sanghuang
如何对基因组序列进行注释_基因组注释_徐洲更hoptop的博客-CSDN博客
使用BRAKER2进行基因组注释(v 2.1.5版)_徐洲更hoptop的博客-CSDN博客
基于PASA进行基因预测 - 简书
PASA的使用 | 陈连福的生信博客
使用Trinity进行转录组组装 - 简书
Trinity的安装与使用 | 陈连福的生信博客
Stringtie的使用说明 - 简书
Genome Guided Trinity Transcriptome Assembly · trinityrnaseq/trinityrnaseq Wiki · GitHub
真核生物基因组的基因分析和预测_基因结构预测_wangyunpeng_bio的博客-CSDN博客
使用MAKER进行注释: 训练SNAP基因模型_徐洲更hoptop的博客-CSDN博客
03 从头(ab initio)预测基因
breaker2
mkdir 10.gene_prediction && cd 10.gene_prediction
mkdir breaker2 && cd breaker2首先使用hisat2根据屏蔽后的参考序列建立索引，进行比对。
# 制作软连接
ln /media/aa/DATA/SZQ2/bj/my/genome/S20191019_23/04genome_feature_analysis_primary/repeat_analysis/repeatMasker/genome.fasta.masked genome.fa
# 新建文件夹
mkdir index
# 进入conda环境
conda activate rnaseq
# 建立索引
hisat2-build genome.fa index/masked
# 新建文件夹
mkdir barker2
# 3个转录组
hisat2 -p 20 -x index/masked -1 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210A/5.1210A_1.clean.fq.gz -2 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210A/5.1210A_2.clean.fq.gz | samtools sort -@ 10 > barker2/A_masked.bam
hisat2 -p 20 -x index/masked -1 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210B/5.1210B_1.clean.fq.gz -2 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210B/5.1210B_2.clean.fq.gz | samtools sort -@ 10 > barker2/B_masked.bam
hisat2 -p 20 -x index/masked -1 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210C/5.1210C_1.clean.fq.gz -2 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210C/5.1210C_2.clean.fq.gz | samtools sort -@ 10 > barker2/C_masked.bam
然后，以未屏蔽重复序列的参考序列和BAM文件作为输入，让BRAKER2（安装会稍显麻烦，因为依赖许多软件）进行预测。
sudo su
密码
conda activate training
# 配置
export AUGUSTUS_CONFIG_PATH=/root/anaconda3/envs/training/config/
export AUGUSTUS_BIN_PHTH=/root/anaconda3/envs/training/bin/
export AUGUSTUS_SCRIPTS_PATH=/root/anaconda3/envs/training/bin/
export BAMTOOLS_PATH=/root/anaconda3/envs/training/bin/
export SAMTOOLS_PATH=/root/anaconda3/envs/training/bin/
export ALIGNMENT_TOOL_PATH=/media/aa/DATA1/software/gth-1.7.3-Linux_x86_64-64bit/bin/
export BLAST_PATH=/media/aa/DATA1/software/rmblast-2.10.0/bin/
export GENEMARK_PATH=/media/aa/DATA/SZQ2/gmes_linux_64_4
export DIAMOND_PATH=/root/anaconda3/envs/training/bin
export PYTHON3_PATH=/root/anaconda3/envs/training/bin
export CDBTOOLS_PATH=/root/anaconda3/envs/training/bin
# 测试
braker.pl --checkSoftware
# 运行
braker.pl --gff3 --cores 20 --species=5.1210 --genome=/media/aa/DATA/SZQ2/bj/my/genome/5.1210/03plion_primary/pilon02.fasta --bam=barker2/A_masked.bam,barker2/B_masked.bam,barker2/C_masked.bam
# 解锁
sudo chmod -R 777 *
exit
# 退出
cd ../..
04 RNA-seq的两种使用策略
# 新建文件夹并进入
mkdir 08.RNA-seq_analysis && cd 08.RNA-seq_analysis1）Trinity de novo
# 新建文件夹
mkdir denovo && cd denovo
conda activate rnaseq
# 运行
Trinity --seqType fq --CPU 20 --max_memory 64G --jaccard_clip --left /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210A/5.1210A_1.fq.gz,/media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210B/5.1210B_1.fq.gz,/media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210C/5.1210C_1.fq.gz --right /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210A/5.1210A_2.fq.gz,/media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210B/5.1210B_2.fq.gz,/media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210C/5.1210C_2.fq.gz &> trinity_denovo.log最终得到Trinity.fasta文件
2）基于HISAT2 + StringTie
①首先，使用HISAT2将RNA-seq数据比对到参考基因组, 这一步和之前相似，但是要增加一个参数--dta，使得StingTie能更好的利用双端信息
# 新建文件夹
mkdir index
# 进入conda环境
conda activate rnaseq
# 建立索引
hisat2-build /media/aa/DATA/SZQ2/bj/my/genome/5.1210/10.gene_prediction/genome.fa index/masked
# 新建文件夹
mkdir hisat2
# 3个转录组
hisat2 --dta -p 20 -x index/masked -1 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210A/5.1210A_1.fq.gz -2 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210A/5.1210A_2.fq.gz | samtools sort -@ 10 > hisat2/A.bam
hisat2 --dta -p 20 -x index/masked -1 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210B/5.1210B_1.fq.gz -2 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210B/5.1210B_2.fq.gz | samtools sort -@ 10 > hisat2/B.bam
hisat2 --dta -p 20 -x index/masked -1 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210C/5.1210C_1.fq.gz -2 /media/aa/DATA/SZQ2/bj/nuohe/RNA-seq/5.1210C/5.1210C_2.fq.gz | samtools sort -@ 10 > hisat2/C.bam
# merge
samtools merge -@ 10 hisat2/merged.bam hisat2/A.bam hisat2/B.bam hisat2/C.bam
②然后用StringTie进行转录本预测
# 进入conda环境
conda activate rnaseq
# 运行
stringtie -p 10 -o hisat2/merged.gtf hisat2/merged.bam
得到merged.gtf
③对于后续的EvidenceModeler而言，它不需要UTR信息，只需要编码区CDS，需要用TransDecoder进行编码区预测
# 从GTF文件中提取FASTA序列
gtf_genome_to_cdna_fasta.pl hisat2/merged.gtf /media/aa/DATA/SZQ2/bj/my/genome/5.1210/04genome_feature_analysis_primary/repeat_analysis/repeatMasker/genome.fasta.masked > hisat2/yc.transcripts.fasta得到yc.transcripts.fasta
3）Trinity genome-guided
Trinity --seqType fq --max_memory 64G --CPU 20 --jaccard_clip --genome_guided_bam hisat2/merged.bam --genome_guided_max_intron 10000 --output trinity_genomeGuide --bflyCalculateCPU &> trinity_genomeGuide.log最终得到Trinity-GG.fasta文件
4）PASA构建综合转录组数据库
# 新建文件夹
mkdir zonghe && cd zonghe
cp ../denovo/trinity_out_dir/Trinity.fasta Trinity.fasta
cp ../trinity_genomeGuide/Trinity-GG.fasta Trinity-GG.fasta
cp ../hisat2/yc.transcripts.fasta yc.transcripts.fasta
# 合并文件夹
$ cat Trinity.fasta Trinity-GG.fasta yc.transcripts.fasta > transcripts.fasta
$ perl -e 'while (<>) { print "$1\n" if />(\S+)/ } ' Trinity.fasta > tdn.accs
# 建立软连接
ln -s /media/aa/DATA/SZQ2/bj/my/genome/5.1210/03plion_primary/pilon02.fasta genome.fasta
对 transcripts 序列进行 end-trimming (vector, adaptor, primer, polyA/T tails)
/media/aa/DATA/SZQ2/PASApipeline-v2.5.2/bin/seqclean transcripts.fasta （使用Univec数据库会报错）统计结果在seqcl_transcripts.fasta.log
构建MYSQL数据库和表-修改PASA的配置文件
# 文件conf.txt
cp /media/aa/DATA/SZQ2/bj/my/genome/3.15188A/08.RNA-seq_analysis/zonghe/conf.txt conf.txt
# 已改动过，不用再改
vim conf.txt
# 文件alignAssembly.config
cp /media/aa/DATA/SZQ2/PASApipeline-v2.5.2/pasa_conf/alignAssembly.config alignAssembly.config
# 需要改动
vim alignAssembly.config
修改 DATABASE=5_1210_db
# 进入conda环境
conda activate rnaseq
# 建立软连接（出现报错，需要做）
ln -s /var/run/mysqld/mysqld.sock /tmp/mysql.sock
# 构建相应数据库（若选择做这一步，则下一步不要输入-C，否则会造成数据库重复）
/media/aa/DATA/SZQ2/PASApipeline-v2.5.2/scripts/create_mysql_cdnaassembly_db.dbi -r -c alignAssembly.config -S /media/aa/DATA/SZQ2/PASApipeline-v2.5.2/schema/cdna_alignment_mysqlschema
# 运行PASA（在/zonghe下，仅进行无参综合转录本的Trinity组装结果;只用blat,用gmap报错）
/media/aa/DATA/SZQ2/PASApipeline-v2.5.2/Launch_PASA_pipeline.pl -c alignAssembly.config -R -g genome.fasta -T -t transcripts.fasta.clean -u transcripts.fasta --ALIGNERS blat --CPU 4 --stringent_alignment_overlap 30.0 --TDN tdn.accs --MAX_INTRON_LENGTH 20000 --TRANSDECODER &> pasa.log构建综合转录组数据库
/media/aa/DATA/SZQ2/PASApipeline-v2.5.2/scripts/build_comprehensive_transcriptome.dbi -c alignAssembly.config -t transcripts.fasta.clean从PASA组装中提取ORF
/media/aa/DATA/SZQ2/PASApipeline-v2.5.2/scripts/pasa_asmbls_to_training_set.dbi --pasa_transcripts_fasta 5_1210_db.assemblies.fasta --pasa_transcripts_gff3 5_1210_db.pasa_assemblies.gff3结果：5_1210_db.assemblies.fasta.transdecoder.genome.gff3
提取蛋白序列长度大于等于300的基因预测结果
/media/aa/DATA/SZQ2/PASApipeline-v2.5.2/scripts/pasa_asmbls_to_training_set.extract_reference_orfs.pl 5_1210_db.assemblies.fasta.transdecoder.genome.gff3 300 > best_candidates.gff3提取完好的基因组模型(两遍)
geneModels2AugusutsTrainingInput --cpu 8 --min_evalue 1e-9 --min_identity 0.7 --min_coverage_ratio 0.7 --min_cds_num 1 --min_cds_length 450 --min_cds_exon_ratio 0.4 --keep_ratio_for_excluding_too_long_gene 0.99 best_candidates.gff3 genome.fasta
# 粘贴统计结果
touch geneModels2AugusutsTrainingInput.result.txt
cd .. 输出：out.filter1.gff3去冗余后的基因模型，out.filter2.gff3完整的基因模型，可用于ab initio基因预测软件的HMM training。
5）准备同源蛋白文件（已下载）
uniprot_sprot_fungi.dat.gz
# 得到/media/aa/DATA/SZQ2/protein.fa
Gene Prediction Commend 01相关推荐
- gene prediction commend 2
使用MAKER进行基因注释(高级篇之AUGUSTUS模型训练)_徐洲更hoptop的博客-CSDN博客https://www.cnblogs.com/zhanmaomao/p/12359964.htm ...
- gene prediction commend 3
若使用本方法,请引用此文,谢谢! JoF | Free Full-Text | Genome Re-Annotation and Transcriptome Analyses of Sanghuang ...
- 程序员们，阿里、腾讯和百度的公司职级、薪资待遇，你有了解吗？
前言 相信程序员们已经度过了一个非常愉快的5.1假期,假期过后就要投入到工作中了,在这愉快的日子里给大家分享一下,一线大厂阿里.腾讯.百度的互联网公司级别和薪资待遇,希望能够给大家增加一些信心,能够努 ...
- 条件随机场（Conditional random fields，CRFs）文献阅读指南
与最大熵模型相似,条件随机场(Conditional random fields,CRFs)是一种机器学习模型,在自然语言处理的许多领域(如词性标注.中文分词.命名实体识别等)都有比较好的应用效果.条 ...
- Prodigal-原核生物基因预测
文章目录 Prokaryotic gene feature What is Prodigal? Installing Prodigal Installing on Mac OS X Installin ...
- Imitation- Seahorse genome
source Abstract report the sequencing de novo assembly of the genome 比较分析得出进化速率 Comparative genomic ...
- 宏基因组分箱整合工具 DAS Tool从零学起笔记
参考https://github.com/cmks/DAS_Tool DAS: dereplication, aggregation and scoring strategy DAS Tool可以将不 ...
- metaProdigal：宏基因组序列中的基因和翻译起始位点预测
文章目录 metaProdigal:宏基因组序列中的基因和翻译起始位点预测 热心肠日报 摘要 动机 Motivation 结果 Results 可用性 Availability 主要结果 表1. 大肠 ...
- Prodigal：原核基因识别和翻译起始位点鉴定
文章目录 Prodigal:原核基因识别和翻译起始位点鉴定 热心肠日报 摘要 背景 结果 结论 实现 Implementation 图1. Prodigal算法的伪代码方式描述 表1. Prodiga ...
最新文章
- jQuery事件的链式写法
- Android studio 实验过程中遇到的问题之android.support.v7.app.AppCompatActivity不能使用的解决办法
- 帮你躲坑：pip install tensorflow 报错怎么办，import tensorflow 报错怎么办？
- linux内核奇遇记之md源代码解读之二
- android 离散分布控件,Android自定义睡眠质量分布控件
- [bzoj4453]cys就是要拿英魂！
- 原子互换：一统公链江湖的神来之笔
- Mol Cell Proteomics. |廖文丽| 阿尔兹海默症临床前期的脑脊液中突触蛋白的变化先于神经变性标志物...
- 投资中最简单的事--阅读笔记
- 小灰，你出书花了多少钱？
- 联通数科一面+二面+面谈 经验分享 base济南
- 哈工大形式语言与自动机2022期末试题
- 收购游戏手机厂商黑鲨背后，腾讯走了一步好棋？
- Pytorch线性模型初体验
- 经典排序算法之快速排序（二分法排序）
- python制作圣诞贺卡_个性化的圣诞贺卡
- BIM技术为什么广泛应用于地铁工程项目上？
- 【博客118】有趣的谬论—龟兔谬论
- 使用lodop实现web精确套打
- Linux- 系统随你玩之--微服务应用出现极少概率会时断时续，它抽风了吗？