当我们需要做非模式物种的富集时，进行查询后也没有对应的数据，那么需要我们自己准备。那就需要这个物种的基因组数据fasta 和gtf 文件自行构建。
这里我以斑马鱼为例。斑马鱼其实有对应的数据库，只是最近刚好下载了斑马鱼的基因组，其他没有数据库的按照此方法也一样可以。斑马鱼数据库，如下图：
image.png
一、基因组和软件环境准备
1、参考基因组下载：
我比较常用的数据库是NCBI和ensembl，根据需求选择对应的参考基因组版本。我选择的是ensembl
#先下载好该物种的fasta 和gtf文件
wget -c http://ftp.ensembl.org/pub/release-106/fasta/danio_rerio/dna/Danio_rerio.GRCz11.dna.toplevel.fa.gz
wget -c http://ftp.ensembl.org/pub/release-106/gtf/danio_rerio/Danio_rerio.GRCz11.106.gtf.gz
2、使用到的相关软件安装：
从conda安装时最简单的方法，若是安装不上就下载包手动安装。kobas安装
conda install -y samtools gffread ucsc-gtftogenepred transdecoder diamond interproscan kobas hmmer
interproscan安装过程中conda会提示：数据库很大，因此未在此安装中提供。请按照以下命令手动下载并安装数据库。官网
######################################
# First time usage please README !!! #
######################################
The databases are huge and consequently not shiped within this installation.
Please download and install the Databases manually by following the commands below:
!!! /!\ Edit the 2 first lines to match the wished version of the DB /!\ !!!
Commands:
=========
# See here for latest db available: https://github.com/ebi-pf-team/interproscan or http://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/
# Set version
version_major=5.54
version_minor=87.0
# get the md5 of the databases
wget http://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/-/interproscan---64-bit.tar.gz.md5
# get the databases (with core because much faster to download)
wget http://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/-/interproscan---64-bit.tar.gz
# checksum
md5sum -c interproscan---64-bit.tar.gz.md5
# untar gz
tar xvzf interproscan---64-bit.tar.gz
# remove old DB
rm -rf /data/pipeline/scRNA/miniconda/envs/chipseq/share/InterProScan/data/
# copy past the new db
cp -r interproscan--/data /data/pipeline/scRNA/miniconda/envs/chipseq/share/InterProScan/
根据提示进行以下操作：
wget -c http://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/5.54-87.0/interproscan-5.54-87.0-64-bit.tar.gz.md5
wget -c http://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/5.54-87.0/interproscan-5.54-87.0-64-bit.tar.gz
# checksum
md5sum -c interproscan-5.54-87.0-64-bit.tar.gz.md5
# untar gz
tar xvzf interproscan-5.54-87.0-64-bit.tar.gz
# remove old DB
rm -rf /data/pipeline/scRNA/miniconda/envs/chipseq/share/InterProScan/data/
# copy past the new db
cp -r interproscan-5.54-87.0/data /data/pipeline/scRNA/miniconda/envs/chipseq/share/InterProScan/
如果conda安装失败，可以直接手动安装
#解压
tar xvzf interproscan-5.54-87.0-64-bit.tar.gz
# 解压完后，进入目录
cd interproscan-5.54-87.0
#查看是否安装完好，若有用法说明弹出则表示安装成功。
./interproscan.sh
# 进行初始化，此命令大致是准备好HMM模型以供hmmscan使用
python3 initial_setup.py
#测试，第一个测试应该使用默认文件名 test_all_appl.fasta.tsv 创建一个输出文件，然后第二个测试将创建 test_all_appl.fasta_1.tsv（因为默认文件名已经存在）。
./interproscan.sh -i test_all_appl.fasta -f tsv -dp
./interproscan.sh -i test_all_appl.fasta -f tsv
二、从参考基因组提取各种信息
1、使用samtools faidx 建立索引
#如果报错Cannot index files compressed with gzip, please use bgzip，那么把gz文件解压即可
samtools faidx Danio_rerio.GRCz11.dna.toplevel.fa
2、使用Cufflinks里面的工具gffread，将转录本的序列从对应的参考基因组上提取出来
gffread Danio_rerio.GRCz11.106.gtf.gz -g Danio_rerio.GRCz11.dna.toplevel.fa -w Danio_rerio_transcript.fa
3、将GTF转化为GenePred 文件
gtfToGenePred -ignoreGroupsWithoutExons -allErrors Danio_rerio.GRCz11.106.gtf Danio_rerio_Ensemble_flat.txt
gene.seq.fa
5、提取gene序列中最长的ORF
TransDecoder.LongOrfs -t gene.seq.fa
6、通过diamond 进行蛋白查询
diamond blastp --query gene.seq.fa.transdecoder_dir/longest_orfs.pep --db swissport.dmnd --max-target-seqs 1 --outfmt 6 --evalue 1e-5 --threads 8 --out blastp_outfmt6.txt
# dmnd文件的准备
wget -c ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
diamond makedb --in uniprot_sprot.fasta.gz --db swissport
7、结构域注释
#下载pfam数据库，其实interproscan安装包中有pfam数据库，可以直接使用
ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam33.1/Pfam-A.hmm.gz
gzip -d Pfam-A.hmm.gz
# 得到 PFAM 数据库的 HMM 文件。 HMM 文件是文本文件，需要将其变成二进制格式，以加快运算速度，同时进行压缩，并建立成索引数据库。
hmmpress Pfam-A.hmm
#将蛋白结构域的比对结果以表格形式输出到指定的文件中
hmmscan --cpu 8 --domtblout pfam_domtblout.txt -o pfam.log Pfam-A.hmm gene.seq.fa.transdecoder_dir/longest_orfs.pep
#预测可能的编码区
TransDecoder.Predict -t gene.seq.fa --retain_pfam_hits pfam_domtblout.txt --retain_blastp_hits blastp_outfmt6.txt
处理gene.seq.fa.transdecoder.pep文件，得到pep的fatsa文件：
gene.seq.fa.transdecoder.pep
pep.fa
8、使用蛋白序列进行GO 注释，获得gene id 对应的GO号
interproscan.sh -t p -i pep.fa -f TSV -b pep -T go -iprlookup -goterms -pa -appl Pfam,SUPERFAMILY
网友评论