✅基因组注释：基因注释，就是在“读完基因组序列”之后，给这些 A、T、G、C 的字母加上“解释”和“标签”的过程。
也就是说，不只是知道这段 DNA 的“字母顺序”，还要弄清楚：
这里是不是一个基因？
从哪儿开始、在哪儿结束？
它能不能翻译成蛋白质？
它跟什么功能、什么通路、什么疾病有关？
通过结构注释（Structural annotation）解决“长什么样、在哪儿”的问题，比如：
基因的位置、长度
外显子 / 内含子边界
是否有可变剪接转录本
各种元素：启动子、UTR、增强子等
本期分享利用Augustus、Barrnap、cmscan、tRNAscan、Infernal预测exon、intron、cds、lncRNA、rRNA、tRNA及一些其他小rRNA。
git clone https://github.com/Gaius-Augustus/Augustus.gitcd Augustusmake augustus#如果想编译 Augustus 的附加工具（utilities / 辅助程序）make auxprogs
mamba create -n genemark \-c conda-forge -c bioconda \perl perl-yaml perl-hash-merge perl-json perl-json-pp \perl-file-which perl-libwww-perl perl-devel-size -ymamba activate genemarkmamba install -c conda-forge perl-parallel-forkmanagermamba install -c conda-forge perl-mce#-----下载GenemarkES的安装包----tar -zxvf gmes_linux_64.tar.gz#----下载key安装包gunzip gm_key_64.gzmv gm_key_64 .gm_key
网址：http://topaz.gatech.edu/GeneMark/license_download.cgi
就可以得到tar.gz文件以及他的key文件（这个也是必不可少的）
git clone https://github.com/tseemann/barrnap.gitcd barrnap/bin./barrnap --help
网址：https://github.com/tseemann/barrnap
git clone https://github.com/UCSC-LoweLab/tRNAscan-SE.gitcd tRNAscan-SE/pwd./configure --prefix=$pwdmake && make installecho 'export PERL5LIB=/software/tRNAscan-SE/lib:$PERL5LIB' >> ~/.bashrcecho 'export PATH=$PATH:/software/tRNAscan-SE/bin' >> ~/.bashrcsource ~/.bashrc#----调用infernal-----wget -c http://eddylab.org/infernal/infernal-1.1.2-linux-intel-gcc.tar.gztar -zxvf infernal-1.1.5-linux-intel-gcc.tar.gzcd infernal-1.1.5-linux-intel-gcc/pwd./configure --prefix=$pwdmake && make installcp ./binaries/* bin#---tRNAscan-SE.conf---修改一些文件夹位置bin_dir: /software/tRNAscan-SE/binlib_dir: /software/tRNAscan-SE/lib/infernal_dir: /software/infernal-1.1.5-linux-intel-gcc/bin
网址：https://github.com/UCSC-LoweLab/tRNAscan-SE
http://eddylab.org/infernal/
#----下载Rfam CM库以及Rfam clanin文件----wget https://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.claninwget https://ftp.ebi.ac.uk/pub/databases/Rfam/CURRENT/Rfam.cm.gzgunzip Rfam.cm.gz#---利用 Infernal里面的 cmpress 为Rfam.cm建立索引./cmpress /Rfam/Rfam.cm
git clone https://github.com/gpertea/gffreadcd gffreadmake release
git clone https://github.com/gpertea/gffcomparecd gffcomparemake release
gffcompare 是一个用于比较、合并、注释和估计一个或多个GFF文件准确性的工具。它可以将新组装的转录本与参考注释进行比较，帮助我们发现新的转录本。
主要输出文件
.gtf文件：主要输出注释文件。
.combined.gtf：当提供多个GTF/GFF文件时，gffcompare会生成一个包含所有样本中所有转录本的“并集”的GTF文件。
.annotated.gtf：如果提供单个GTF/GFF查询文件作为输入，并且没有指定删除“重复的”/“冗余”转录本的特定选项，则gffcompare将输出一个名为*.annotated.gtf的文件。
.loci：报告每个转录本中每个新的转录本位点信息。
.stats：报告了与参考注释相比的输入转录本的准确性（或一致性）相关的各种统计数据。
.tracking：匹配样本之间的转录本。
.refmap：对于每个参考转录本，查询的转录本是完全还是部分匹配到参考转录本上。
.tmap：对于每个转录本最佳的参考转录本。
wget http://cpc2.gao-lab.org/data/CPC2-beta.tar.gztar xzvf CPC2-beta.tar.gzcd CPC2-betaexport CPC_HOME="$PWD"cd libs/libsvmgzip -dc libsvm-3.18.tar.gz | tar xf -cd libsvm-3.18make clean && makecd $CPC_HOME
网址：https://github.com/gao-lab/CPC2_standalone
#利用RepeatModeler和RepeatMasker去除重复序列BuildDatabase -name index/index_chr1 -engine ncbi /DATA/pre_chr1.faRepeatModeler -database index/index_chr1 -engine ncbi -pa 20 &> index_chr1.out#这一步有一个.masked结尾的文件是屏蔽了重复序列的用这个文件后续分析RepeatMasker -lib index/index_chr1-families.fa -e ncbi -pa 10 -dir . /DATA/pre_chr1.fa#--这里的chr1.fa就是去重以后的---#---利用已知训练集预测------augustus --species=chicken chr1.fa > ./chr1.gff#---利用GenemarkES先训练再用Augustus训练，最后用Augustus预测#---这一步得到genemark.gtf用于autoAug.pl训练perl gmes_linux_64/gmes_petap.pl --ES --sequence chr1.fasta --soft_mask auto --cores 56Augustus/scripts/autoAug.pl --genome=chr1.fasta --species=species --trainingset=GenemarkES/genemark.gtf --useexisting#---这条命令一定要有，是为了把自己训练的物种路径指定不然会报错export AUGUSTUS_CONFIG_PATH=/software/Augustus/config~/software/Augustus/bin/augustus --species=silkie chr1.fasta > genemark_auto_Aug_chr1.gff
👉这里AUGUSTUS用的chr1.fa是经过RepeatModeler+RepeatMasker去重复后的序列。
barrnap/bin/barrnap --kingdom euk --threads 20 --outseq rRNA_chr1.fa chr1.fasta >rRNA_chr1.gff3❗ 关键参数：--kingdom
选择数据库（HMM 模型）：
bac→ 细菌arc→ 古菌euk→ 真核生物 ✔mito→ 线粒体（仅用于线粒体基因组）
🔧 常用参数解释
1）--threads
指定线程数（加速）
2）--reject
低于期望长度比例的 rRNA 不输出
默认 0.25，通常足够
大基因组可调高（例如 0.50）避免碎片：
3）--evalue
HMMER 搜索的 E-value cutoff：
4）--lencutoff
定义 partial（部分匹配）的比例，默认就行。
software/tRNAscan-SE-install/bin/tRNAscan-SE -o Data/chr1.gff Data/chr1.fa
❗结果包含染色体号、基因组坐标、tRNA Type（氨基酸对应类型）、Anti-codon（反密码子）、Intron Begin / End
如果 tRNA 有内含子，这里给出内含子坐标、Inf Score（Infernal Score），是预测可信度。
cmscan \--cpu 32 \--cut_ga \--rfam \--noali \--tblout rfam.tbl \--fmt 2 \--clanin /Rfam/Rfam.clanin /Rfam/Rfam.cm > chr1.cmscan
✔--cpu 32
多线程，提高速度（根据你节点核心数改）。
✔ --cut_ga
告诉 Infernal 使用 Rfam 模型内置的 GA cutoff（可信度阈值），避免大量假阳性。
这是做 ncRNA 注释必须开的选项。
✔ --rfam
自动使用 Rfam 模型注释风格（比如假阳性过滤、incl/excl rules）。
做基因组注释时必须开。
✔ --noali
不输出序列比对（alignment），减少输出体积 100 倍。
GFF 注释不需要 alignment。
✔ --tblout rfam.tbl
生成总结表格（非常重要）。
你最后需要的 snoRNA、snRNA 就在这里。
✔ --fmt 2
输出格式 = Rfam 风格，更好解析。
✔--clanin Rfam.cm 是下载并 cmpress 的 CM 模型 以及clanin文件
右边 > rfam.cmscan 是详细结果
这是cmscan结果。包含了基因组坐标及对应的ncRNA类型。
#----对fna文件构建索引---hisat2-build chr1.fa chr1#-----map比对---hisat2 --new-summary -p 2 -x chr1 -1 A_R1.fastq.gz -2 A_R2.fastq.gz -S A.sam --rna-strandness RFsamtools view -@ 4 -b -o A.bam A.samsamtools sort -@ 24 -O bam -o A_sort.bam A.bam#--stringtie-----stringtie -p 10 -G chr1.gff -l A_stringtie -o A__stringtie_no_e.gtf A_sorted.bamstringtie --merge -p 10 -G chr1.gff -o stringtie_merged.gtf list.txt
先对fna文件建立索引→hisat2比对samtools转化为bam文件→stringtie对每个个体组装转录组，因为是从头组装不要加-e参数 →stringtie merge参数把每个个体得到的gtf文件合并一下
gffcompare -r chr.gff -o gffcompare_chr1 stringtie_merged.gtf对上一步得到的stringtie_merged.gtf用gffcompare鉴定一下转录本类型
class code 用于表示输入中的转录本与注释中的转录本的关系。常见的class code包括：
=：内含子链的完全匹配。
c：预测转录本包含在参考转录本中。
j：潜在的新异构体（片段），至少一个剪接连接与参考转录本共享。
e：单外显子转录本与参考外显子和参考内含子至少10 bp重叠，表明可能存在前mRNA片段。
i：完全比对到参考内含子内的转录本。
o：部分外显子与参考转录本重叠。
p：可能的聚合酶片段，在参考转录本2k bp范围内。
r：重复序列。
u：未知，基因间转录本。
x：与参考转录本在相反链上的外显子重叠。
lncRNA一般是
常用做 lncRNA 的 class_code：
u：intergenic（典型 lincRNA）x：antisensei：完全在已知 intron 内o：其他重叠j：新的剪接组合（可选，看你要不要）
awk 'NR>1 && $3 ~ /^[uxioj]$/ {print $5}' \gffcomp_merged.stringtie_merged.gtf.tmap | sort -u > novel_ids.txt#将uxioj类型的得到gtfgrep -F -f <(sed 's/^/transcript_id "/; s/$/"/' novel_ids.txt) \stringtie_merged.gtf > novel_raw.gtf#做长度 ≥ 200 nt + 外显子数 ≥ 2 过滤awk -F'\t' '$3=="exon" {# 从第9列属性里拆出 transcript_idn = split($9, a, ";")tid = ""for (i=1; i<=n; i++) {gsub(/^ +| +$/, "", a[i]) # 去掉前后空格if (a[i] ~ /^transcript_id /) {split(a[i], b, " ")tid = b[2]gsub(/"/, "", tid) # 去掉双引号}}if (tid != "") {exon_count[tid]++if (!(tid in start) || $4 < start[tid]) start[tid]=$4if (!(tid in end) || $5 > end[tid]) end[tid]=$5}}END {for (tid in exon_count) {len = end[tid] - start[tid] + 1if (len >= 200 && exon_count[tid] >= 2) {print tid}}}' novel_raw.gtf > lnc_ids.txt#将lnc的transcriptid得到gtf以便后面得到fastagrep -F -f <(sed 's/^/transcript_id "/; s/$/"/' lnc_ids.txt) \novel_raw.gtf > lncRNA_candidate.gtf
gffread lncRNA_candidate.gtf -g chr1.fasta -w lncRNA_candidate.fapython2.7 CPC2.py -i lncRNA_candidate.fa -o lncRNA_cpc2
#ID transcript 的名字（就是你 fasta 里的 ID，比如 MSTRG.1.2）
transcript_length 转录本长度（nt，核苷酸个数）
peptide_length 预测到的最长 ORF 的氨基酸长度（aa）
Fickett_score Fickett test code score，反映序列像不像编码序列（越大越“像”）
pI 这个 ORF 翻译出来的蛋白的等电点（isoelectric point）
ORF_integrity ORF 完整性：1 = 有完整起止密码子，-1 = 没有完整 ORF :contentReference[oaicite:0]{index=0}
coding_probability 该转录本是“编码”的概率（0~1），CPC2 模型算出来的 越接近 1 越“像编码”，接近 0 越“像非编码”
label 最终分类结果：`coding` 或 `noncoding`
将noncoding的转录本提取出来即为lncRNA。
这样就会得到一个基因组文件的结构注释信息。好啦，本期分享就到这结束啦~仅供参考，感谢大家关注支持，祝大家科研顺利~
