生信菜鸟团博客2周年精选文章集(6)三个最基础生信软件教程

其实我现在已经不写软件教程了！

fastqc对原始测序reads质控 NCBI的blast++软件使用说明书 SRA工具sratoolkit把原始测序数据转为fastq格式

目录

一：下载安装该软件

二：准备数据

三：运行命令

四：输出文件解读

正文

一：下载安装该软件

在NCBI的ftp站点里面可以找到blast++的下载链接

wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/LATEST/ncbi-blast-2.2.30+-x64-linux.tar.gz

我们一般选择适合我们操作系统的二进制版本，解压即可使用

可以把它们添加到PATH，前提是有root权限，或者把该目录添加到PATH也行。

cp * /home/jmzeng/my-bin/bin/

我把my-bin添加到了我的PATH，所以可以直接使用这些程序了 二：准备数据

只需要fasta文件的数据即可，query和target都可以是该fasta文件，可以随便找两个fa文件做测试

三：运行命令

1，建库，用makeblastdb，标准是

makeblastdb -in db.fasta -dbtype prot -parse_seqids -out dbname

具体参数看help里面的，但是我们一般用这几个就够了的

我的例子 ：对200M的蛋白文件

makeblastdb -in uniprot_sprot.trinotate_v2.0.pep -dbtype prot -parse_seqids -out sprot

输出的文件如下，基本不需要看，反正调用的时候只用sprot这个

对8G的uniref90，

makeblastdb -in uniprot_uniref90.trinotate_v2.0.pep -dbtype prot -parse_seqids -out uniref90

2，比对分为好几种，blastn, blastp,blastx，tblastn，tblastx

• blastp -query seq.fasta -out seq.blast -db dbname -outfmt 6 -evalue 1e-5 -num_descriptions 10 -num_threads 8
• blastn -query seq.fasta -out seq.blast -db dbname -outfmt 6 -evalue 1e-5 -num_descriptions 10 -num_threads 8
• blastx -query seq.fasta -out seq.blast -db dbname -outfmt 6 -evalue 1e-5 -num_descriptions 10 -num_threads 8

参数说明:

• -query： 输入文件路径及文件名
• -out：输出文件路径及文件名
• -db：格式化了的数据库路径及数据库名
• -outfmt：输出文件格式，总共有12种格式，6是tabular格式对应BLAST的m8格式
• -evalue：设置输出结果的e-value值
• -num_descriptions：tabular格式输出结果的条数
• -num_threads：线程数

四：输出文件解读

重点是-outfmt 6，也就是之前版本的m 8格式

结果中从左到右每一列的意义分别是：

• [00] Query id
• [01] Subject id
• [02] % identity
• [03] alignment length
• [04] mismatches
• [05] gap openings
• [06] q. start
• [07] q. end
• [08] s. start
• [09] s. end
• [10] e-value
• [11] bit score

一，下载该软件

wget http://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-ubuntu64.tar.gz

tar xzf sratoolkit.current-centos_linux64.tar.gz

解压直接使用即可，里面有一大堆的软件，针对不同的测序仪，不同的数据

我一般只用/home/jmzeng/down_software/sratoolkit.2.3.5-2-ubuntu64/bin/fastq-dump

/home/jmzeng/down_software/sratoolkit.2.3.5-2-ubuntu64/bin/fastq-dump –split-3 SRR1793917.sra

二：下载数据

首先去NCBI里面搜索并找到你想要的数据的SRA地址，然后写脚本批量下载。

如果文献里面的SRA号，那么可以直接打开NCBI里面的搜索界面下载

如果文献里面是SRP号，那么该SRP会涉及到好几个SRA数据，得一个个开网站下载

三：用命令解压数据

下载之后的数据是

非常简单的命令，就可以把当前文件夹下的所有sra都解压开来！

[shell]

for i in *sra do echo $i /home/jmzeng/bio-soft/sratoolkit.2.3.5-2-ubuntu64/bin/fastq-dump –split-3 $i done

[/shell]

解压的同时它也会显示每个SRA文件的数据量

四：结果文件解读

可以看到，每个SRA文件都产生了两个reads，分别是左右两端测序，说明这个SRA文件是双端测序策略。

随便打开一个fastq文件可以看到，它的读长是300bp

fastqc软件的使用

一：下载安装该软件

具体搜索其地址下载，fastqc是一个java软件，下载后可以直接使用，但是需要自行配置好java环境，具体配置方法，见linux下java配置。

二：准备数据

数据就是我们测序得到的fastq文件的reads，压缩包也可以直接运行

三：运行命令

我习惯了批处理解决问题，脚本如下：

for id in *fastq

do

echo $id

/home/jmzeng/bio-soft/FastQC/fastqc $id

Done

运行过程中会显示以下的提示信息

估计还是要运行很久的，比较这几个RNA-seq文件每个都是16G的

按住ctrl+A+D即可退出该后台，继续去前台执行简单任务

好像二十分钟就跑完了

输出文件如下

四：输出文件解读

可以直接打开那个html网页文件就可以查看每一个图片内容，也可以解压那个zip压缩包具体看每一张图片

下载fastqc跑出来的结果一个个解读

1，简单统计表格

这些英文我就不翻译了，reads均长是100bp，共四千多万条reads

2，测序质量图

这个图其实很容易看，就是100bp长度reads上的1-100的坐标在这四千万条reads里面的测序质量的箱线图，看那个红线均值就可以了，超过Q30就蛮好了，超过Q20也是合格的

3，碱基（A,T,C,G）含量图

这也是100bp长度reads上的1-100的坐标在这四千万条reads里面的A,T,C,G的比例，如果是全基因组全转录组的随机打断，那么就应该A,T,C,G的比例都接近于25%，如果测序是有目的性的，那么比例也就相应的改变了

4，reads的GC含量频数分布图

这是对四千万条reads里面的GC含量值做统计密度曲线，可以看到绝大部分的reads的GC含量都集中在50%附近。极端情况很少。

5，reads长度分布图

可以看到大多reads都是100bp长度，很整齐

6，可能的重复序列表格

可以看到这些重复序列比例很高，高达千分之一，而且被注释了可能的来源，adapter，是需要去除的。