收集数据
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
-
利用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
结果文件一览
-
利用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评价构建转录本的准确性和完整性
首先官网下载对应的同源拷贝数据库到本地,并解压
(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使用