RNA-seq数据分析完全指北-04:创建本地blast库分析物种组成

2021-04-13 14:49:26 浏览数 (1)

啊~~~本来是半个月的专栏不知道到底过了多久才又和大家见面,其中经历不足为外人道也

在上一步骤中,我怀疑最后一个样本可能有一些其他物种的污染。当时没想到应该怎样处理。后来周运来师兄提到了一个解决办法。

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 本地化教程》——不知道怎么取名字【博客园】

0 人点赞