基因组GTF优化全流程

收集数据

NCBI下载相同物种的参考FASTA和GTF注释文件

PelSin_1.0 (Wang Z et al., 2014) —只有scaffold水平的基因组，暂无genome水平

我们这里有三个下机cyclone数据

将数据挂载到任务下

/data/input/Files/cyclone_2/TB20002087/TB20002087-202401121501020_read.fq.gz

/data/input/Files/cyclone_2/TB2000209F/TB2000209F-202401121501061_read.fq.gz

/data/input/Files/cyclone_2/TB200033B9/TB200033B9-202401121501520_read.fq.gz

参考基因组路径

/data/work/turtle/NCBI_Pelodiscus_sinensis/GCF_000230535.1_PelSin_1.0_genomic.fna（FASTA）

/data/work/turtle/NCBI_Pelodiscus_sinensis/genomic.gtf（GTF）

利用Isoquant进行数据比对和转录本预测

考虑到基因组较大，采用离线脚本任务投递

脚本投递示例

isoquant.py --reference /data/work/turtle/NCBI_Pelodiscus_sinensis/GCF_000230535.1_PelSin_1.0_genomic.fna --complete_genedb --genedb /data/work/turtle/NCBI_Pelodiscus_sinensis/genomic.gtf --fastq /data/input/Files/cyclone_2/TB20002087/TB20002087-202401121501020_read.fq.gz --data_type nanopore -o /data/work/turtle/isoquant_output

三组测序分别输出到三个文件夹isoquant_output111/222/333

离线任务日志示例

。。。。。。

结果路径文件一览

Output222

Output333

1. 利用gffcompare比较/合并注释文件

比较对象路径

/data/work/turtle/turtle.20231106.cdy.gtf

用picbio三代注释（chen）和我们构建的转录本比较

脚本投递示例

gffcompare -r /data/work/turtle/turtle.20231106.cdy.gtf -o /data/work/turtle/gffcompare_cdyto111 /data/work/turtle/isoquant_output111/OUT/OUT.extended_annotation.gtf

单独比较结果文件

合并原有注释和新构建注释

gffcompare -r /data/work/turtle/turtle.20231106.cdy.gtf -o /data/work/turtle/cdy_111 -D /data/work/turtle/isoquant_output111/OUT/OUT.extended_annotation.gtf /data/work/turtle/turtle.20231106.cdy.gtf

结果文件一览

1. 利用Transdecoder预测基因结构

1）从上一步重新生成的GTF中提取FASTA序列

gtf_genome_to_cdna_fasta.pl /data/work/turtle/combined_cdy_111/cdy_111.combined.gtf /data/work/turtle/NCBI_Pelodiscus_sinensis/GCF_000230535.1_PelSin_1.0_genomic.fna > transcripts_111.fasta

2）将GTF格式转换为GFF3格式（以备后续使用）

import sys

inFile = open(sys.argv[1],'r')

for line in inFile:
  #skip comment lines that start with the '#' character
  if line[0] != '#':
    #split line into columns by tab
    data = line.strip().split('\t')

    #parse the transcript/gene ID. I suck at using regex, so I usually just do a series of splits.
    #transcript_id和gene_ids是GTF第九列肯定存在的
    transcriptID = data[-1].split('transcript_id')[-1].split(';')[0].strip()[1:-1]
    geneID = data[-1].split('gene_id')[-1].split(';')[0].strip()[1:-1]

    #replace the last column with a GFF formatted attributes columns
    #I added a GID attribute just to conserve all the GTF data
    data[-1] = "ID=" + transcriptID + ";GID=" + geneID

    #print out this new GFF line
    print('\t'.join(data))

3）预测转录本中的长开放阅读框（need some time）

TransDecoder.LongOrfs -t /data/work/turtle/transdecoder_111/transcripts_111.fasta

（该步骤跑完的结果路径一览）

4）使⽤DIAMOND对上⼀步输出的transcripts.fasta.transdecoder.pep在蛋⽩数据库中进⾏搜索，寻找同源证据⽀持

下载数据库并解压

https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/u

niprot_sprot.fasta.gz

gzip -d /data/work/turtle/uniprot_sprot.fasta.gz

建立索引

diamond makedb --in /data/work/turtle/transdecoder_111/uniprot_sprot.fasta --db swissprot

BLASTP⽐对

diamond blastp -d /data/work/turtle/transdecoder_111/swissprot.dmnd -q /data/work/turtle/transdecoder_111/transcripts_111.fasta.transdecoder_dir/longest_orfs.pep --evalue 1e-5 --max-target-seqs 1 > blastp.outfmt6

预测可能的编码区

TransDecoder.Predict -t /data/work/turtle/transdecoder_111/transcripts_111.fasta --retain_blastp_hits /data/work/turtle/transdecoder_111/blastp.outfmt6

（注意这一步要把之前预测ORF时产生的transcripts_111.fasta.transdecoder_dir文件夹放在当前路径下，因为结果文件会被继续添加进这个文件夹）

一些过程报错

5）生成基于新转录本的编码区注释文件

cdna_alignment_orf_to_genome_orf.pl /data/work/turtle/transdecoder_111/transcripts_111.fasta.transdecoder.gff3 /data/work/turtle/transdecoder_111/cdy_111.combined.gff3 /data/work/turtle/transdecoder_111/transcripts_111.fasta > transcripts.fasta.transdecoder.genome.gff3

结果文件保存在

/data/work/turtle/transdecoder_111/transcripts.fasta.transdecoder.genome.gff3

less /data/work/turtle/transdecoder_111/transcripts.fasta.transdecoder.genome.gff3

6）基因重命名

less /data/work/turtle/transdecoder_111/transcripts.fasta.transdecoder.genome.gff3 | awk '{OFS="\t";if($3=="gene" && $9~"%7"){split($9,a,"%7C");split($9,b,"=");$9=b[1]"="b[2]"="a[3]"#"a[2];print;} else if($3=="gene" && $9!~"%7"){split($9,a,";");split(a[1],b,"=");$9=a[1]";Name="b[2];print;}else{print;}}' >transcripts.fasta.transdecoder.genome.renamegene.gff3

基因组注释文件结构

BUSCO评价构建转录本的准确性和完整性

首先官网下载对应的同源拷贝数据库到本地，并解压

Index of /v5/data/lineages/

（Eukaryota真核生物，或者更具体的Metazoa动物）

tar -zxvf /data/work/turtle/BUSCO/eudicots_odb10.2024-01-08.tar.gz (-C output_path)

# 提取最长转录本

perl /path/to/trinity/opt/trinity-2.4.0/util/misc/get_longest_isoform_seq_per_trinity_gene.pl /path/to/transcripts.fasta > longest_transcripts.fasta

最后运行busco代码如下：

busco -i transcripts_111.fasta -l ./eudicots_odb10 -o buscops -m tran --cpu 8 --offline &

（最好都放在当前路径下，避免绝对路径）

（以及若不能联网，要加上--offline &离线运行参数，否则报错如下）

结果

使用基因浏览器IGV进行可视化

IGV使用