啊~~~本来是半个月的专栏不知道到底过了多久才又和大家见面,其中经历不足为外人道也
。
在上一步骤中,我怀疑最后一个样本可能有一些其他物种的污染。当时没想到应该怎样处理。后来周运来师兄提到了一个解决办法。
1、创建本地blast库
平常来说,我们使用NCBI的在线blast功能就能解决绝大多数问题。但是,当我们有大量序列需要对比时,在线查找已经不能满足我们。因此,我们需要构建一个本地blast库。常用的库有nt库和nr库,其中nt库是核酸序列,nr库是氨基酸序列。
1.1、安装ncbi-blast
NCBI的本地blast软件已经构建的非常好了,只需要下载就可以直接使用。
进入网址https://ftp.ncbi.nlm.nih.gov/blast/executables/LATEST/,选择已经编译好的版本。
代码语言:javascript复制## 下载安装包
wget https://ftp.ncbi.nlm.nih.gov/blast/executables/LATEST/ncbi-blast-2.11.0 -x64-linux.tar.gz
tar zxvf ncbi-blast-2.11.0 -x64-linux.tar.gz ## 解压
mv ncbi-blast-2.11.0 ~/biosoft/blast/ ## 移动到自己常用的软件目录并改名
1.2、配置环境变量
代码语言:javascript复制vim ~/.bashrc
## 在最后一行加上如下信息,记得改成自己的路径
export PATH=/home/csw/biosoft/blast/bin:$PATH
## 如果不会使用vim编辑器,那就运行如下命令
echo "export PATH=/home/csw/biosoft/blast/bin:$PATH" >> ~/.bashrc
## 让配置生效
source ~/.bashrc
1.3、设置blast
在根目录下创建一个.ncbirc文件,并添加如下内容
代码语言:javascript复制## 创建.ncbirc文件
vim ~/.ncbirc
## 将下列内容粘贴进去即可,记得改成自己的路径
# Start the section for BLAST configuration
[BLAST]
# Specifies the path where BLAST databases are installed
BLASTDB=/home/csw/reference/linux/blast
# Specifies the data sources to use for automatic resolution
# for sequence identifiers
DATA_LOADERS=blastdb
# Specifies the BLAST database to use resolve protein sequences
BLASTDB_PROT_DATA_LOADER=/home/csw/reference/linux/blast/nr
# Specifies the BLAST database to use resolve protein sequences
BLASTDB_NUCL_DATA_LOADER=/home/csw/reference/linux/blast/nt
BATCH_SIZE=10G
# Windowmasker settings
[WINDOW_MASKER]
WINDOW_MASKER_PATH=/home/csw/reference/linux/blast/windowmasker
# end of file
1.4、使用update_blastdb安装预构建的nt库
代码语言:javascript复制## update_blastdb --showall可以显示所有可以安装的库
nohup update_blastdb --decompress nt & ## 后台下载,在下载中断后会自动断点续传,--decompress表示会在下载后自动解压
## 慢慢等吧,文件挺大的
## 其他参数参考update_blastdb --help
2、本地对比
2.1、从FQ文件中抽取10000条随机序列并转换为FA格式
代码语言:javascript复制ln -s ~/projects/rna-seq/ ./
seqkit sample -n 10000 SRR11178353_1.fastq.gz | seqkit fq2fa - > SRR11178353_1.rd.fa
seqkit sample -n 10000 SRR11178353_2.fastq.gz | seqkit fq2fa - > SRR11178353_2.rd.fa
2.2、使用本地blast库进行对比
代码语言:javascript复制nohup blastn -query SRR11178353_1.rd.fa -out SRR11178353_1.blast -db nt -outfmt "6 std scomname" -evalue 1e-5 -num_threads 16 -qcov_hsp_perc 50.0 -num_alignments 5 &
nohup blastn -query SRR11178353_2.rd.fa -out SRR11178353_2.blast -db nt -outfmt "6 std scomname" -evalue 1e-5 -num_threads 16 -qcov_hsp_perc 50.0 -num_alignments 5 &
2.3、使用R语言进行统计
主要对比到的是灵长类序列,应该是由同源相似性导致的。爷佛了。
硬着头皮继续走流程吧,这个问题我再仔细研究一下。
本文部分代码及思路来自以下文章,在此致谢 《使用本地nt数据库对reads和Trinity结果进行blast》——周小钊【简书】 《Linux系统中NCBI BLAST 本地化教程》——不知道怎么取名字【博客园】