宏基因组分析

之前一直都是自学有关宏基因组分析的流程，现在组里向学校租用了服务器，也开始了真正的实操。以下是我大致的数据预处理分析流程，从数据下载到非冗余基因集构建。数据是从NCBI中下载的。以下基于Linux操作系统（分给我256g内存 24核，1000gb ssd）。

 

 1.数据下载（sratoolkit）

 sratoolkit的下载安装

sra是储存二代测序原始数据的数据库，在NCBI官网上提供专门下载sra数据的工具sratoolkit，附上网址。https://github.com/ncbi/sra-tools/wiki/01.-Downloading-SRA-Toolkit

其中可以根据你的系统选择适合的下载连接，我服务器用的是redhat的，这里下载的是centos （如果是自己电脑装的Ubuntu虚拟机则下载Ubuntu版，根据自己的来）。也可以下载win版，然后用远程连接工具将数据传入到服务器，我用的是WinSCP。

 

①获取下载地址后，在Linux系统中用wget命令下载安装包。

wget"https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/3.0.7/sratoolkit.3.0.7-centos_linux64.tar.gz"

②进入包存放目录，进行解压。

tar xzvf sratoolkit.3.0.7-centos_linux64.tar.gz

③配置环境。

vim ~/.bashrc
在末尾添加：export PATH=“$PATH:/XXXX/software/sratoolkit.3.0.0-ubuntu64/bin”;
source ~/.basrc

安装成功。 

 

 从NCBI下载数据

 找到你想要的sra号，比如我这里的SRR9302957，然后运行以下代码，即可获取数据。

prefetch SRR9302957

将sra格式转换为fastq格式 

 sratoolkit还可以将sra格式的数据转换成fastq格式，以便后续分析。代码如下：

fastq_dump --splite-e SRR9302957.sra
# 这里是双端测序，单端的可查看详细参数说明。

 

2.去接头过滤数据 （trimmomatic）

 

 trimmomaatic下载安装及运行

下载地址：http://www.usadellab.org/cms/index.php?page=trimmomatic 

下载完成，解压，trimmomatic是基于java编写的，确保你系统已安装java。运行代码如下：

## 双端测序数据使用方法
1.  java -jar trimmomatic-0.32.jar PE \
2.		[-threads <threads>] \ # 线程数
3.		[-phred33|-phred64] \ # 设置碱基的质量格式
4.		[-trimlog <trimLogFile>] \ # 生成过滤日志
5.		[<inputFile1> <inputFile2>] \
6.		[<outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>] \
7.		[ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>] \ # 接头序列，在adapters目录下
8.		[SLIDINGWINDOW:<windowSize>:<requiredQuality>] \
9.		[LEADING:<quality>] \
10.		[TRAILING:<quality>] \
11.		[MINLEN:<length>] # 切除后的序列最短长度

# 举例:
1. java -jar trimmomatic-0.32.jar PE \
2.		-threads 24 \
3.		-phred33 \
4.		-trimlog trim.log \
5.		input_forward_R1.fq.gz input_reverse_R2.fq.gz \
6.		output_forward_paired_R1.fq.gz  output_forward_unpaired_R1.fq.gz output_reverse_paired_R2.fq.gz output_reverse_unpaired_R2.fq.gz \
7.		ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true \
8.		SLIDINGWINDOW:5:20 \
9.		LEADING:3 \
10.		TRAILING:3 \
11.		MINLEN:36

 

3.短序列组装（megahit） 

 

megahit的下载安装 

我是用的conda安装，关于conda的下载，不赘述。

conda install -c bioconda megahit # 安装完成即存在miniconda的bin目录下

 使用

# 双端序列组装：
	megahit -1 pe_1.fq -2 pe_2.fq -o outdir
	#-1：pair-end 1序列，-2 pair-end 2序列，-o输出目录
# 单端序列：
	megahit -r single_end.fq -o out
# 交错的双端序列：
	megahit --12 interleaved.fq -o out

###
常用参数
两种设置kmer参数的方法：
1）–k-list：组装的kmer size列表，支持多kmer组装，不同kmer size之间逗号分隔，可设置的范围15-255，相邻kmer size间隔必须小于或等于28，默认为21,29,39,59,79,99,119,141
2）组合参数：–k-min：设置最小的kmer size，应小于255，必须为奇数，默认为21
–k-max：设置最大的kmer size，应小于255，必须为奇数，默认为141
–k-step：多kmer组装的kmer size间隔，应小于等于28必须为偶数，默认为12
###

 

 4.组装结果评估（quast）

 

 quast的下载安装

过程如下：

wget https://sourceforge.net/projects/quast/files/quast-5.2.0.tar.gz
tar xzvf quast-5.2.0.tar.gz

使用

quast是基于python开发的，确保系统已装python。

cd quast-5.2.0 # 进入quast目录
python quast.py -o quast_evaluation ~/data/megahit_out/contigs.fasta

 

5.基因预测（prodigal） 

 

prodigal的下载安装 

我是用conda安装的，直接运行：

conda install -c bioconda prodigal
# prodigal即存在miniconda3下的bin目录里

 软件使用

-a 是输出氨基酸文件
-d 输出预测基因的序列文件
-f 选择输出文件格式，有gbk,gff,和sco格式可供选择
-g 指定密码子，原核为第11套
-i 输入文件，即需要预测的基因组序列文件
-m 屏蔽基因组中的N碱基
-o 输出文件，默认为屏幕输出
-p 选择方式，是单菌还是meta样品
-q 不输错错误信息到屏幕
-t 指定训练集
-s 输出所有潜在基因以及分值到一个文件中

prodigal -a protein.fasta -d nucleotide.fasta -f gff -g 11 -o genes.gff -s poteintial.stat -i contig.fasta

 

6.基因去冗余（cd-hit） 

 

cd-hit下载安装 

 我也是用的conda下载

conda install -c bioconda cd-hit

软件使用 

可先将生成的所有蛋白序列合并为一个序列：

cat a.fasta b.fasta c.fasta > all.fasta

 运行：（如果运行cd-hit-est则是对核酸序列操作）

cd-hit -i all.fasta -o new.fa -c 0.8 -aS 0.8 -d 0

###
-i：输入文件，fasta格式。
-o：输出文件前缀，输出文件有两个，分别为fasta格式序列文件和以.clstr结尾的聚类信息文件。
-c：较短序列比对到长序列的bp与自身bp数的比值超过该数值则聚类为一组，默认为0.9。
-d：聚类信息文件中各个聚类组中序列名的长度，设为0则将取完整序列名。
-T 使用的线程数。
-aL：控制代表序列比对严格程度的参数，默认为0，若设为0.8则表示比对区间要占到代表（长）序列的80%。
-aS：控制短序列比对严格程度的参数，默认为0，若设为0.8则表示比对区间要占到短序列的80%。
###

cd-hit运行完会生成两个文件，一个是含所有去冗余后的序列 ，fasta格式；一个是以.clstr结尾的聚类信息文件。

 到此为止，即构建好了非冗余基因集，可用于接下来的物种、基因功能等分析。