RNA_seq植物實戰(zhàn)
Author : yujia
目錄：
- 概述
- salmon工具完成索引建立和生物學定量
- subread工具完成序列比對和定量
- DESeq2差異基因分析
- 總結
一蝇裤、 概述
- 練習數(shù)據(jù)：數(shù)據(jù)來源于擬南芥，共16個樣本粘室，處理分為4組（0day,1day,2day,3day）
- 練習目的：熟悉兩套RNA-seq差異基因表達分析的流程（salmon流程和subread流程）
- 數(shù)據(jù)存放地址：/public/study/mRNAseq/tair/
- 軟件調用地址：/usr/local/bin/miniconda3/bin/
- 實戰(zhàn)項目來源地址：Jimmy學長的生信菜鳥團博客：http://www.bio-info-trainee.com/2809.html （一個植物轉錄組項目的實踐）
二、salmon工具完成索引建立和生物學定量
??salmon是一款不通過序列比對就可以快速完成生物學定量的RNA-seq數(shù)據(jù)分析工具卜范。它的使用流程包括兩步：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
執(zhí)行后舱殿，會得到如下文件：
-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流程到此就結束了嗡害，根據(jù)得到的quant文件焚碌，我們可以在后續(xù)利用DESeq2, edgeR, limma等包進行下游的差異基因分析。現(xiàn)在我們來看subread工具如何完成RNA-seq數(shù)據(jù)的生物學定量霸妹。
三十电、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就是我們比對得到的文件良价，可以發(fā)現(xiàn)寝殴，目錄下還有.vcf和.bed文件，那是因為subjunc可以檢測indel和exon-exon junction
現(xiàn)在我們就可以開始做生物學定量了明垢。
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 -
它的基本語法如下溯革，以雙端測序數(shù)據(jù)為例（具體詳情可見官網(wǎng)）：
#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是我們后續(xù)用作差異基因分析的文件冈闭。
-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為例号显，來進行后續(xù)的差異基因分析臭猜。
目前常見的差異基因分析工具有DESeq2躺酒、limma包等等，此處以DESeq2為工具來進行差異基因的篩選蔑歌。
First step:獲取表達矩陣和分組信息
進行差異基因分析之前羹应，首先要獲取表達矩陣和分組信息。我們的表達矩陣是剛才用featureCounts定量得到的counts_id.txt 次屠，經(jīng)過格式處理之后园匹，是這樣（部分截取）：
<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) #輸入表達矩陣數(shù)據(jù)路徑
us_count<-round(us_count,digits=0) #將輸入數(shù)據(jù)取整
#準備
us_count<-as.matrix(us_count) #將數(shù)據(jù)轉換為矩陣格式
condition<-factor(c("day0","day0","day0","day0","day1","day1","day1","day1")) ## 設置分組信息供汛，建立環(huán)境(8個樣本，2組處理)
coldata<-data.frame(row.names=colnames(us_count),condition)
coldata #展示coldata內容
condition #展示環(huán)境
上面的代碼運行之后涌穆，我們的分組信息就是這樣噠：
<img src="C:\Users\admin\Desktop\分組信息.png" width="600" hegiht="400" align=center />
現(xiàn)在怔昨，就可以正式使用DESeq2做差異基因分析了，總共其實只有三步：
- 構建dds矩陣
- 對dds矩陣進行標準化
- 提取結果并繪制火山圖
代碼如下：
library(DESeq2) #使用library函數(shù)加載DEseq2包
##構建dds矩陣
dds<-DESeqDataSetFromMatrix(us_count,coldata,design=~condition)
head(dds) #查看構建好的矩陣
##進行差異分析
dds<-DESeq(dds) #對原始的dds進行標準化
resultsNames(dds) #查看結果名稱
res<-results(dds) #用results函數(shù)提取結果宿稀，并賦值給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 />
五祝沸、 總結
自己之前處理線蟲數(shù)據(jù)主要是用的RSEM（加bowtie2）來做RNA-seq的序列比對和生物學定量矮烹，但是沒有接觸過salmon和subread。這次使用這兩個工具完成了植物的RNA-seq實戰(zhàn)罩锐，對這些工具有了更深的理解奉狈，以后自己如果要開發(fā)相關軟件，也要多用用唯欣，比較工具之間的優(yōu)劣嘹吨。
