RNA_seq植物实战
Author : yujia
目录：
- 概述
- salmon工具完成索引建立和生物学定量
- subread工具完成序列比对和定量
- DESeq2差异基因分析
- 总结
一、 概述
- 练习数据：数据来源于拟南芥，共16个样本，处理分为4组（0day,1day,2day,3day）
- 练习目的：熟悉两套RNA-seq差异基因表达分析的流程（salmon流程和subread流程）
- 数据存放地址：/public/study/mRNAseq/tair/
- 软件调用地址：/usr/local/bin/miniconda3/bin/
- 实战项目来源地址：Jimmy学长的生信菜鸟团博客：http://www.bio-info-trainee.com/2809.html （一个植物转录组项目的实践）
二、salmon工具完成索引建立和生物学定量
salmon是一款不通过序列比对就可以快速完成生物学定量的RNA-seq数据分析工具。它的使用流程包括两步：1.建立索引 2.对reads进行生物学定量（quantification）。所以，如果我们使用salmon工具来做生物学定量的话，会非常的快速简洁。以下是代码流程：
First step:使用salmon < index >选项对转录组建立索引
利用salmon建立索引的基本语法是：
salmon index -t athal.fa.gz -i athal_index
#index 代表建立索引
#-t .fa文件的路径
#-i 索引存放路径
所以我们的代码如下：
/usr/local/bin/miniconda3/bin/salmon index -t /public/study/mRNAseq/tair/Arabidopsis_thaliana.TAIR10.28.cdna.all.fa.gz -i /trainee/home/yjxiang/practice/index_file
执行后，会得到如下文件：
-rw-rw-r-- 1 yjxiang yjxiang 14K Aug 13 14:40 duplicate_clusters.tsv
-rw-rw-r-- 1 yjxiang yjxiang 751M Aug 13 14:40 hash.bin
-rw-rw-r-- 1 yjxiang yjxiang 357 Aug 13 14:40 header.json
-rw-rw-r-- 1 yjxiang yjxiang 115 Aug 13 14:40 indexing.log
-rw-rw-r-- 1 yjxiang yjxiang 412 Aug 13 14:40 quasi_index.log
-rw-rw-r-- 1 yjxiang yjxiang 116 Aug 13 14:40 refInfo.json
-rw-rw-r-- 1 yjxiang yjxiang 7.8M Aug 13 14:40 rsd.bin
-rw-rw-r-- 1 yjxiang yjxiang 247M Aug 13 14:40 sa.bin
-rw-rw-r-- 1 yjxiang yjxiang 63M Aug 13 14:40 txpInfo.bin
-rw-rw-r-- 1 yjxiang yjxiang 96 Aug 13 14:40 versionInfo.json
Second step:使用salmon < quant >选项完成生物学定量
salmon中进行生物学定量的基本语法是：
salmon quant -i athal_index -l A -1 samp_1.fastq.gz -2 samp_2.fastq.gz -p 8 -o quants/sample_quant
# quant是salmon中进行生物学定量的选项
# -i The -i argument tells salmon where to find the index
# -l A tells salmon that it should automatically determine the library type of the sequencing reads
#The -1 and -2 arguments tell salmon where to find the left and right reads for this sample
# -p 8 argument tells salmon to make use of 8 threads
# -o argument specifies the directory where salmon’s quantification results sould be written
由于我们一共有16个样本，要一一进行生物学定量，因此编写一个bash脚本完成批量处理：
#脚本地址存储在/trainee/home/yjxiang/practice
#!/bin/bash
index=/trainee/home/yjxiang/practice/index_file
for fn in ERR1698{194..209};
do
samp=`basename ${fn}`
echo "Processing sample ${samp}"
/usr/local/bin/miniconda3/bin/salmon quant -i $index -l A \
-1 ${samp}_1.fastq.gz \
-2 ${samp}_2.fastq.gz \
-p 6 -o /trainee/home/yjxiang/practice/quants/${samp}_quant
done
运行脚本之后，得到如下结果：
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:28 ERR1698194_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:38 ERR1698195_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:39 ERR1698196_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:39 ERR1698197_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:40 ERR1698198_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:40 ERR1698199_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:41 ERR1698200_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:42 ERR1698201_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:43 ERR1698202_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:43 ERR1698203_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:44 ERR1698204_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:44 ERR1698205_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:45 ERR1698206_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:46 ERR1698207_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:46 ERR1698208_quant
drwxrwxr-x 5 yjxiang yjxiang 4.0K Aug 13 15:47 ERR1698209_quant
每个文件夹里都有对应样本的quant结果，以样本ERR1698209为例,文件夹里含有这些文件，quant.sf 文件就是我们得到的定量结果：
drwxrwxr-x 2 yjxiang yjxiang 4.0K Aug 13 15:47 aux_info
-rw-rw-r-- 1 yjxiang yjxiang 307 Aug 13 15:47 cmd_info.json
-rw-rw-r-- 1 yjxiang yjxiang 551 Aug 13 15:47 lib_format_counts.json
drwxrwxr-x 2 yjxiang yjxiang 4.0K Aug 13 15:47 libParams
drwxrwxr-x 2 yjxiang yjxiang 4.0K Aug 13 15:00 logs
-rw-rw-r-- 1 yjxiang yjxiang 1.8M Aug 13 15:47 quant.sf
查看一下quant.sf，常见的TPM值、Numreads都在里面：
$ head -n 5 quant.sf
Name Length EffectiveLength TPM NumReads
ATMG00010.1 462 301.089 0.000000 0.000000
ATMG00030.1 324 166.891 0.000000 0.000000
ATMG00040.1 948 786.477 0.000000 0.000000
ATMG00050.1 396 236.034 0.000000 0.000000
salmon流程到此就结束了，根据得到的quant文件，我们可以在后续利用DESeq2, edgeR, limma等包进行下游的差异基因分析。现在我们来看subread工具如何完成RNA-seq数据的生物学定量。
三、subread工具完成序列比对和生物学定量
subread是一个能快速进行序列比对和转录组定量的工具，我们使用它进行转录组定量一共需要三个步骤：1.建立索引 2.序列比对 3.使用featureCounts进行定量 。代码流程如下：
First step:对基因组建立索引
subread建立索引的基本语法如下：
#Build an index for the reference genome (you may provide a single FASTA file including all the reference sequences)：
subread-buildindex -o my_index chr1.fa chr2.fa ...
所以我们的代码如下：
#先解压基因组文件到我的个人目录里
#gunzip -c /public/study/mRNAseq/tair/Arabidopsis_thaliana.TAIR10.28.dna.genome.fa.gz > /trainee/home/yjxiang/practice/raw_data/Arabidopsis_thaliana.TAIR10.28.dna.genome.fa
#然后使用subread软件的subread-buildindex完成索引的建立
/trainee/home/yjxiang/tools/subread-1.6.2-source/bin/subread-buildindex -o /trainee/home/yjxiang/practice/subread_workflow/Ara_genome_index_file /trainee/home/yjxiang/practice/raw_data/Arabidopsis_thaliana.TAIR10.28.dna.genome.fa
运行代码后，得到的索引文件如下：
-rw-rw-r-- 1 yjxiang yjxiang 29M Aug 16 12:29 Ara_genome_index_file.00.b.array
-rw-rw-r-- 1 yjxiang yjxiang 231M Aug 16 12:29 Ara_genome_index_file.00.b.tab
-rw-rw-r-- 1 yjxiang yjxiang 629 Aug 16 12:29 Ara_genome_index_file.files
-rw-rw-r-- 1 yjxiang yjxiang 363K Aug 16 12:29 Ara_genome_index_file.log
-rw-rw-r-- 1 yjxiang yjxiang 76 Aug 16 12:29 Ara_genome_index_file.reads
Second step:序列比对
我们使用subread里集成的subjunc来完成序列比对，关于subjunc的官方说明如下：
The Subjunc aligner is an RNA-seq read aligner, specialized in detecting exon-exon junctions and performing full alignments for the reads (exon-spanning reads in particular).
subjunc的使用基本语法：
#Report uniquely mapped reads only (by default). Mapping output includes BAM files and exon-exon junctions discovered from the data.
subjunc -T 5 -i my_index -r reads1.txt -o subjunc_results.bam
#Report up to three alignments for each multi-mapping read:
subjunc --multiMapping -B 3 -T 5 -i my_index -r reads1.txt -o subjunc_results.bam
#Detect indel of up to 16bp:
subjunc -I 16 -i my_index -r reads1.txt -o subjunc_results.bam
#Map paired-end reads and discover exon-exon junctions:
subjunc -d 50 -D 600 -i my_index -r reads1.txt -R reads2.txt -o subjunc_results.bam
由于样本量较多，所以我们也写一个脚本，来完成这个序列比对任务：
#bash脚本批量完成：
#! /bin/bash
subjunc="/trainee/home/yjxiang/tools/subread-1.6.2-source/bin/subjunc";
index='/trainee/home/yjxiang/practice/subread_workflow/index_file/Ara_genome_index_file';
for fn in ERR1698{194..209};
do
sample=`basename ${fn}`
echo "Processing sample ${sampe}"
$subjunc -i $index \
-r ${sample}_1.fastq.gz \
-R ${sample}_2.fastq.gz \
-T 5 -o /trainee/home/yjxiang/practice/subread_workflow/mapping_out/${sample}_subjunc.bam
done
序列比对完成之后，我们会得到如下文件（截取部分）：
-rw-rw-r-- 1 yjxiang yjxiang 799M Aug 16 12:57 ERR1698194_subjunc.bam
-rw-rw-r-- 1 yjxiang yjxiang 1.1M Aug 16 12:57 ERR1698194_subjunc.bam.indel.vcf
-rw-rw-r-- 1 yjxiang yjxiang 2.6M Aug 16 12:57 ERR1698194_subjunc.bam.junction.bed
-rw-rw-r-- 1 yjxiang yjxiang 1.1G Aug 16 13:01 ERR1698195_subjunc.bam
-rw-rw-r-- 1 yjxiang yjxiang 1.4M Aug 16 13:01 ERR1698195_subjunc.bam.indel.vcf
-rw-rw-r-- 1 yjxiang yjxiang 3.6M Aug 16 13:01 ERR1698195_subjunc.bam.junction.bed
.bam就是我们比对得到的文件，可以发现，目录下还有.vcf和.bed文件，那是因为subjunc可以检测indel和exon-exon junction
现在我们就可以开始做生物学定量了。
Third step:featureCounts完成生物学定量
featureCounts是一款可以进行生物学定量的软件，它集成在subread的package里了，官方关于它的说明如下：
featureCounts is a highly efficient general-purpose read summarization program that counts mapped reads for genomic features such as genes, exons, promoter, gene bodies, genomic bins and chromosomal locations. It can be used to count both RNA-seq and genomic DNA-seq reads.
我们就使用它来完成定量操作。需要注意的是：使用featureCounts，我们需要提供gtf格式的注释或者SAF格式的注释，大概像这样：
GeneID Chr Start End Strand
497097 chr1 3204563 3207049 -
497097 chr1 3411783 3411982 -
497097 chr1 3660633 3661579 -
它的基本语法如下，以双端测序数据为例（具体详情可见官网）：
#Summarize paired-end reads and count fragments (instead of reads):
featureCounts -p -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_PE.bam
所以，我们的代码如下：
gff3='/public/study/mRNAseq/tair/Arabidopsis_thaliana.TAIR10.28.gff3.gz'
gtf='/public/study/mRNAseq/tair/Arabidopsis_thaliana.TAIR10.28.gtf.gz'
featureCounts='/trainee/home/yjxiang/tools/subread-1.6.2-source/bin/featureCounts'
nohup $featureCounts -T 5 -p -t exon -g gene_id -a $gtf -o /trainee/home/yjxiang/practice/subread_workflow/counts_out/counts_id.txt /trainee/home/yjxiang/practice/subread_workflow/mapping_out/*.bam &
代码运行后，就可以查看我们的定量结果了：
featureCounts会输出两份文件，一个是summary,一个是count reads的结果.counts_id.txt是我们后续用作差异基因分析的文件。
-rw-rw-r-- 1 yjxiang yjxiang 6.3M Aug 16 13:52 counts_id.txt
-rw-rw-r-- 1 yjxiang yjxiang 2.3K Aug 16 13:52 counts_id.txt.summary
四、DESeq2差异基因分析
获得reads counts之后，我们就可以开展差异基因分析了。我们以subread中的featureCounts工具得到的counts_id.txt为例，来进行后续的差异基因分析。
目前常见的差异基因分析工具有DESeq2、limma包等等，此处以DESeq2为工具来进行差异基因的筛选。
First step:获取表达矩阵和分组信息
进行差异基因分析之前，首先要获取表达矩阵和分组信息。我们的表达矩阵是刚才用featureCounts定量得到的counts_id.txt ，经过格式处理之后，是这样（部分截取）：
<img src="C:\Users\admin\Desktop\表达矩阵.png" width="600" hegiht="400" align=center />
第一列是基因ID，之后的列都是样本ID
每一行代表不同的基因在不同样本中的表达量.
我们选day0和day1做比较，
为了方便，分组矩阵的制作我在R里面完成，输入如下代码：
us_count<-read.csv("C:\\Users\\admin\\Desktop\\RNA_seq_Ara_counts_day1_day0.csv",head=T,row.names=1) #输入表达矩阵数据路径
us_count<-round(us_count,digits=0) #将输入数据取整
#准备
us_count<-as.matrix(us_count) #将数据转换为矩阵格式
condition<-factor(c("day0","day0","day0","day0","day1","day1","day1","day1")) ## 设置分组信息，建立环境(8个样本，2组处理)
coldata<-data.frame(row.names=colnames(us_count),condition)
coldata #展示coldata内容
condition #展示环境
上面的代码运行之后，我们的分组信息就是这样哒：
<img src="C:\Users\admin\Desktop\分组信息.png" width="600" hegiht="400" align=center />
现在，就可以正式使用DESeq2做差异基因分析了，总共其实只有三步：
- 构建dds矩阵
- 对dds矩阵进行标准化
- 提取结果并绘制火山图
代码如下：
library(DESeq2) #使用library函数加载DEseq2包
##构建dds矩阵
dds<-DESeqDataSetFromMatrix(us_count,coldata,design=~condition)
head(dds) #查看构建好的矩阵
##进行差异分析
dds<-DESeq(dds) #对原始的dds进行标准化
resultsNames(dds) #查看结果名称
res<-results(dds) #用results函数提取结果，并赋值给res变量
summary(res) #查看结果
plotMA(res,ylim=c(-2,2))
mcols(res,use.names=TRUE)
plot(res$log2FoldChange,res$pvalue) #绘制火山图
#提取差异基因
res <- res[order(res$padj),]
resdata <-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
deseq_res<-data.frame(resdata)
up_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange > 1)) #提取上调差异表达基因
down_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange < -1)) #提取下调差异表达基因
write.csv(up_diff_result,"C:\\Users\\admin\\Desktop\\上调_day0_VS_day1_diff_results.csv") #输出上调基因
write.csv(down_diff_result,"C:\\Users\\admin\\Desktop\\下调_day0_VS_day1_diff_results.csv") #输出下调基因
至此，差异基因就成功提取了，看看火山图：
<img src="C:\Users\admin\Desktop\火山图.png" width="600" hegiht="400" align=center />
五、 总结
自己之前处理线虫数据主要是用的RSEM（加bowtie2）来做RNA-seq的序列比对和生物学定量，但是没有接触过salmon和subread。这次使用这两个工具完成了植物的RNA-seq实战，对这些工具有了更深的理解，以后自己如果要开发相关软件，也要多用用，比较工具之间的优劣。
