今天分享的學(xué)習(xí)筆記是一套轉(zhuǎn)錄組分析簡(jiǎn)單流程劫拗，適用于初學(xué)者入門閱讀间校，從原始測(cè)序數(shù)據(jù)開始，經(jīng)過質(zhì)控页慷、序列比對(duì)憔足、定量表達(dá)胁附、差異表達(dá)、功能富集等一系列分析步驟滓彰，最終獲得基因表達(dá)信息控妻。本文所有數(shù)據(jù)都經(jīng)過特殊修改，僅供學(xué)習(xí)參考使用揭绑。

轉(zhuǎn)錄組是在特定時(shí)空條件下細(xì)胞中基因轉(zhuǎn)錄表達(dá)產(chǎn)物弓候，廣義的轉(zhuǎn)錄組包括信使RNA，核糖體RNA他匪，轉(zhuǎn)運(yùn)RNA及非編碼RNA菇存，狹義上是指所有mRNA的集合，轉(zhuǎn)錄組分析能夠獲得不同基因的表達(dá)情況诚纸。

1. 數(shù)據(jù)來源

假設(shè)有兩個(gè)不同組織（PR和SR）撰筷，每個(gè)組織各區(qū)三個(gè)樣本，一共六個(gè)樣本畦徘，利用illumina平臺(tái)進(jìn)行轉(zhuǎn)錄組測(cè)序毕籽，得到雙端測(cè)序數(shù)據(jù)。數(shù)據(jù)原始格式為.fq,共有12條測(cè)序數(shù)據(jù)文件（每個(gè)樣本產(chǎn)生兩條）

PR1_2.fq  PR2_2.fq  PR3_2.fq  SR1_2.fq  SR2_2.fq  SR3_2.fq
PR2_1.fq  PR3_1.fq  SR1_1.fq  SR2_1.fq  SR3_1.fq

創(chuàng)建工作文件夾井辆，并新建子步驟目錄：00_trainingRawdata 01_trimmomaticFiltering 02_hisat2Mapping 03_featurecountsQuatification 04_DESeq2DEGanalysis

mkdir RNA-Seq_Practice 
cd RNA-Seq_Practice
mkdir 00_trainingRawdata 01_trimmomaticFiltering 02_hisat2Mapping 03_featurecountsQuatification 04_DESeq2DEGanalysis   
#批量創(chuàng)建工作文件夾

2. 測(cè)序數(shù)據(jù)質(zhì)量評(píng)估

利用fastQC軟件對(duì)獲得的fastq序列文件進(jìn)行質(zhì)量分析关筒，生成html格式的結(jié)果報(bào)告，其中含有各項(xiàng)指標(biāo)杯缺，以下用PR1樣本為例蒸播。

fastqc *.fq        #批量對(duì)fq后綴文件運(yùn)行fastqc程序
輸出結(jié)果：PR1_1_fastqc.html  
Filename    PR1_1.fq  
File type   Conventional base calls
Encoding    Illumina 1.5
Total Sequences 105300       #總序列數(shù)
Sequences flagged as poor quality   0
Sequence length 90              #序列長(zhǎng)度
%GC 52                #GC堿基含量

質(zhì)量評(píng)估報(bào)告，使用瀏覽器打開：

3. 過濾低質(zhì)量序列

利用Trimmomatic軟件除去序列文件中的接頭（adapter）萍肆，并對(duì)堿基進(jìn)行合適的修改袍榆，然后對(duì)堿基進(jìn)行修剪，對(duì)低質(zhì)量的序列進(jìn)行過濾塘揣。

time java -jar trimmomatic-0.39.jar PE  #trimmomatic是依賴java環(huán)境運(yùn)行
-threads 1 
-phred33 PR1_1.fq PR1_2.fq 
-summary ../01_trimmomaticFiltering/PR1.summary 
-baseout ../01_trimmomaticFiltering/PR1 LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 HEADCROP:13 MINLEN:36

運(yùn)行程序?qū)π蛄形募械唾|(zhì)量的序列進(jìn)行過濾包雀，將輸出結(jié)果存儲(chǔ)到01_trimmomaticFiltering，每一個(gè)樣本序列會(huì)生成如下輸出文件亲铡，summary文件包含過濾的總結(jié)信息才写。

-rw-r--r-- 1  19M Nov 28 14:29 PR1_1P
-rw-r--r-- 1    0 Nov 28 14:29 PR1_1U
-rw-r--r-- 1  19M Nov 28 14:29 PR1_2P
-rw-r--r-- 1    0 Nov 28 14:29 PR1_2U
-rw-r--r-- 1  282 Nov 28 14:29 PR1.summary

打開其中一個(gè)PR1.summary文件進(jìn)行查看，其中Input Read Pairs表示過濾之前的reads數(shù)奖蔓，Both Surviving Reads表示過濾之后的reads數(shù)赞草。

$ cat PR1.summary  #查看總結(jié)文件
Input Read Pairs: 102300
Both Surviving Reads: 102300
Both Surviving Read Percent: 100.00
Forward Only Surviving Reads: 0
Forward Only Surviving Read Percent: 0.00
Reverse Only Surviving Reads: 0
Reverse Only Surviving Read Percent: 0.00
Dropped Reads: 0
Dropped Read Percent: 0.00

4. 比對(duì)到參考基因組

利用hisat2軟件，將fasta序列比對(duì)到參考基因組吆鹤。首先需要構(gòu)建索引文件厨疙，下載或者拷貝參考基因組序列文件Ref和基因組注釋文件，通過hisat2-build命令構(gòu)建索引文件檀头。

cd xx/RNA-Seq_Practice
cp -r xx/Ref .
cd Ref
gunzip -c chr.fa.gz > chr.fa  #解壓參考基因組
time hisat2-build -p 1 chr.fa Chr   #建立索引文件

完成上述步驟后轰异，將過濾得到的reads比對(duì)到參考基因組上岖沛。輸入文件為兩個(gè)fasta序列文件暑始，將運(yùn)行過程中的輸出提示重定向到log文件搭独。

• -S設(shè)置輸出文件為sam
• -x設(shè)置參考基因組
• -p設(shè)置運(yùn)行的線程數(shù)量

time hisat2 -p 1 #線程數(shù)為1
-x Ref/Chr #參考基因組文件路徑
-1 01_trimmomaticFiltering/PR1_1P #輸入文件
-2 01_trimmomaticFiltering/PR1_2P #輸入文件
-S 02_hisat2Mapping/PR1.sam       #輸出文件sam格式
--new-summary 1>02_hisat2Mapping/PR1_hisat2Mapping.log 2>&1  
#輸出日志

得到的日志文件中包含比對(duì)成功的reads數(shù)量和比對(duì)率等信息：

HISAT2 summary stats:
        Total pairs: 102300
   Aligned concordantly or discordantly 0 time: 94209 (92.09%)
   Aligned concordantly 1 time: 7613 (7.44%)
   Aligned concordantly >1 times: 293 (0.29%)
   Aligned discordantly 1 time: 185 (0.18%)
        Total unpaired reads: 188418
   Aligned 0 time: 186684 (99.08%)
   Aligned 1 time: 1575 (0.84%)
   Aligned >1 times: 159 (0.08%)
        Overall alignment rate: 8.76%   #總比對(duì)率8.76%

上述步驟完成后，會(huì)在當(dāng)前目錄下生成多個(gè)sam格式文件廊镜，SAM（sequence Alignment/mapping)格式是高通量測(cè)序中存放比對(duì)數(shù)據(jù)的標(biāo)準(zhǔn)格式牙肝。另外，bam是sam的二進(jìn)制格式嗤朴，可以減小sam文件的大小配椭，因此利用samtools對(duì)sam進(jìn)一步處理，得到bam文件雹姊，以下用PR1為例股缸，其他樣本按相同方式處理。

time samtools view -bS 
     02_hisat2Mapping/PR1.sam > 02_hisat2Mapping/PR1.bam
time samtools sort 
     02_hisat2Mapping/PR1.bam > 02_hisat2Mapping/PR1.srt.bam

5. 基因表達(dá)定量分析

利用featuresCounts軟件吱雏，對(duì)reads進(jìn)行精確的計(jì)數(shù)敦姻，最后將所有樣本的reads數(shù)合并為一個(gè)文件，將RNA-Seq_Practice_countstable文件導(dǎo)出歧杏，根據(jù)FPKM和TPM的計(jì)算公式定量 每個(gè)基因在每個(gè)樣本中的表達(dá)量镰惦。利用網(wǎng)站（https://www.ncbi.nlm.nih.gov/gene）將基因ID轉(zhuǎn)化為Gene description，從而了解其功能和相關(guān)信息犬绒。

time featureCounts 
-a RefMaize_AGPv4/Zea_mays.AGPv4.32.gtf.gz #基因組注釋文件路徑
-T 1 -p --countReadPairs -g gene_id -t exon 
-o 03featurecountsQuatification/PR1 02hisat2Mapping/PR1.srt.bam

以PR1為例旺入，將多余的信息除去，只保留基因ID和reads數(shù)凯力，得到count文件茵瘾，然后同樣步驟處理其他樣本，最后將所有的reads定量數(shù)據(jù)合并為一個(gè)countstable文件咐鹤。

cat PR1 | cut -f 1,7 > PR1.count
paste PR1.count PR2.count PR3.count SR1.count SR2.count SR3.count > countstable           
#合并不同樣本的定量信息到一個(gè)文件

awk '{$3="";$5="";$7="";$9="";$11="";print $0}' countstable > RNA-Seq_Practice_countstable     
#提取指定列到新文件

將生成的RNA-Seq_Practice_countstable保存到本地拗秘，然后計(jì)算FPKM和TPM值，在R語言中進(jìn)行相關(guān)計(jì)算慷暂。

FPKM（Fragments Per Kilobase of exon model per Million mapped fragments）表示每千個(gè)堿基的轉(zhuǎn)錄每百萬映射讀取的fragments聘殖，該方法是利用每個(gè)樣本的總fragments數(shù)進(jìn)行校正。優(yōu)點(diǎn)是消除樣本間和基因間的差異行瑞，但是該方法主要是計(jì)算的結(jié)果可能與真實(shí)表達(dá)水平有偏差奸腺。

TPM (transcript per million) 表示每百萬轉(zhuǎn)錄本中來自于某基因的轉(zhuǎn)錄本數(shù)目，該方法先對(duì)每個(gè)基因的reads數(shù)用基因的長(zhǎng)度進(jìn)行校正血久，然后再利用校正后的reads數(shù)與校正后該樣本所有的reads數(shù)求商突照。優(yōu)點(diǎn)是消除了外顯子長(zhǎng)度造成的差異和樣本間測(cè)序總reads counts不同造成的差異，缺點(diǎn)為該法不是采用比對(duì)到基因組上的總reads counts氧吐，所以有時(shí)候不夠精確讹蘑。

具體計(jì)算方法如下所述：

> df <- read.table("RNA-Seq_Practice_countstable") #讀入數(shù)據(jù)文件
> rownames(df) <- df$V1   #修改數(shù)據(jù)表的行名為基因ID（第一列V1）
> df <- df[,c(2,3,4,5,6,7)] #現(xiàn)在行名已經(jīng)為基因ID末盔，因此刪除數(shù)據(jù)的第一列
> colnames(df) <- c("PR1","PR2","PR3","SR1","SR2","SR3") #修改行名
> dim(df); names(df)   #查看df的數(shù)據(jù)結(jié)構(gòu)和變量名稱是否正確
[1] 41768     6   
[1] "PR1"    "PR2"    "PR3"    "SR1"    "SR2"    "SR3"   
> head(df[,1:6],2)      #輸出前兩行數(shù)據(jù)進(jìn)行預(yù)覽
             PR1 PR2 PR3 SR1 SR2  #PR和SR共六列數(shù)據(jù)
Zm00001d027   0   0   0   0   0
Zm00001d021   0   0   0   0   0

featureCounts_meta <- df[,1:6]  #提取基因信息
df <- df[rowSums(df)>0,]     #去除表達(dá)量均為0的行
prefix <-"maize_PR_SR"   #設(shè)置輸出文件前綴名

# ----- TPM計(jì)算 ------
kb <- nrow(featureCounts_meta)/ 1000
rpk <- df / kb
tpm <- t(t(rpk)/colSums(rpk) * 1000000)
write.table(tpm, paste0(prefix,"_tpm.xls"), quote=F, sep="\t", row.names=T, col.names=T )

# ----- FPKM計(jì)算 ------
fpkm <- t(t(rpk)/colSums(df) * 10^6)
write.table(fpkm,file= paste0(prefix, "_fpkm.xls"), quote=F, sep="\t", row.names=T, col.names=T )

6. 基因差異表達(dá)分析

利用DESeq2軟件基于樣本的原始reads數(shù)目計(jì)算差異表達(dá)基因，利用run-featurecounts.R腳本對(duì)每個(gè)樣本的reads數(shù)進(jìn)行定量座慰，然后利用abundance_estimates_to_matrix.pl腳本合并所有的定量結(jié)果陨舱，最后利用DESeq2進(jìn)行差異表達(dá)基因的分析。

首先版仔，將所需的腳本文件全部拷貝到工作目錄下游盲，同時(shí)將樣本信息相關(guān)文件也拷貝保存。

cp -r xxx/script_featurecounts .
cp -r xxx/script-mergecounts .
cp xxx/RNA-Seq_Practice/contrast.txt .
cp xxx/RNA-Seq_Practice/sampleinfo.txt .

新建一個(gè)文件夾蛮粮，運(yùn)行R語言腳本益缎，對(duì)樣本進(jìn)行表達(dá)定量分析，以PR1為例然想，輸入文件為02文件夾中生成的bam文件和基因組注釋文件莺奔，生成輸出文件PR1，其他的幾個(gè)樣本按照同樣的方法進(jìn)行處理变泄。

mkdir R-featurecounts
Rscript script_featurecounts/run-featurecounts.R 
-b ../02_hisat2Mapping/PR1.srt.bam  #輸入文件
-g ../Ref/Ref.gtf.gz  #基因組注釋文件
-o R-featurecounts/PR1  #輸出文件

運(yùn)行結(jié)束后令哟，在輸出文件夾內(nèi)將會(huì)生成count文件和log文件，利用perl語言腳本abundance_estimates_to_matrix.pl合并所有樣本的定量結(jié)果杖刷。

perl script-mergecounts/abundance_estimates_to_matrix.pl 
     --est_method featureCounts R-featurecounts/*.count 

利用trinity 軟件中的DESeq2分析差異表達(dá)基因励饵，輸入文件有counts矩陣、所選方法DESeq2滑燃、樣本信息表sampleinfo.txt役听、分組信息contrasts.txt。

perl xxx/run_DE_analysis.pl 
--matrix matrix.counts.matrix   
--method DESeq2 --samples_file sampleinfo.txt  #樣本信息表文件
--contrasts contrasts.txt   #這里需要注意文件名稱是否正確

上述流程輸出的結(jié)果存放在隨機(jī)文件夾中表窘，進(jìn)入該文件夾后典予，可以看到生成的差異表達(dá)基因結(jié)果matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results文件壁熄，然后利用awk命令提取出滿足條件（同時(shí)滿足：|log2FoldChange| > 1实抡，且padj < 0.05）的差異表達(dá)基因,生成新的輸出文件件PRvsSR_DEGs炒刁。

drwxr-xr-x 307 Nov 29 09:54 DESeq2.197264.dir  #隨機(jī)新目錄
$ ls -f  #查看新目錄下文件列表
matrix.counts.matrix.PR_vs_SR.DESeq2.Rscript
matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results  #差異表達(dá)基因結(jié)果文件
matrix.counts.matrix.PR_vs_SR.DESeq2.count_matrix
matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results.MA_n_Volcano.pdf
awk 'sqrt($7*$7)>1 && $11 <0.05 {print $0}'  
./matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results > PRvsSR.DESeq2_DEGs  
#利用awk命令提取滿足條件的指定數(shù)據(jù)

將最終得出的差異表達(dá)基因數(shù)據(jù)文件PRvsSR.DESeq2_DEGs保存到電腦盖喷，進(jìn)行GO富集和KEGG富集分析。繪制火山圖時(shí)使用的是matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results

7. GO富集分析

利用在線網(wǎng)站(http://shiny.host.ugeneyun.cn/GOmaize/)進(jìn)行GO富集分析店煞，輸入文件為差異表達(dá)基因的ID信息表轰胁，此處不進(jìn)行顯著性過濾艰毒，原因是前期已經(jīng)進(jìn)行篩選了既琴。提交運(yùn)行后生成GOenrich結(jié)果占婉，共有270條Term，下圖隨機(jī)展示5條甫恩，可以得出GO編號(hào)逆济、功能描述、顯著性和分類信息。然后奖慌，利用網(wǎng)站的繪圖功能制作GO富集氣泡圖抛虫，保存到本地GO.pdf。

8. KEGG分析

在R語言（4.2.2）中使用clusterProfiler简僧、KEGG.db建椰、ggplot2三個(gè)R包進(jìn)行KEGG富集分析，先構(gòu)建一個(gè)OrgDb本地索引庫涎劈，然后利用enrichKEGG函數(shù)進(jìn)行富集分析广凸，利用在線網(wǎng)站(https://www.maizegdb.org/gene_center/gene)將基因ID轉(zhuǎn)換為GenBank Gene類型阅茶，保存為文件gene.txt蛛枚，輸出文件為pathway數(shù)據(jù)表gene.pathway.csv 和KEGG富集結(jié)果gene.pathway .pdf

remotes::install_github("YuLab-SMU/createKEGGdb",force = TRUE)
createKEGGdb::create_kegg_db('zma')  
#首先使用該軟件制作玉米KEGG富集包
library(clusterProfiler) 
library(ggplot2)   # 繪圖需用R包
library(KEGG.db)   # 該步驟需要手動(dòng)在Rstudio中安裝KEGG.db包

file="gene.txt"  #輸入文件，包含18條差異表達(dá)的基因信息
gene <- read.delim("gene.txt",header = F,row.names = 1)
> head(gene,2)
               GenBank.Gene
Zm00001d0218    1036647
Zm00001d0234    1002408

kk <- enrichKEGG(gene= gene$GenBank.Gene,organism= 'zma',
                 pvalueCutoff = 1,
                 qvalueCutoff = 1,
                 use_internal_data =T)  #進(jìn)行KEGG富集分析
write.csv(kk,file = paste0("gene.pathway.csv"),row.names = F)
p=dotplot(kk);p  #繪制點(diǎn)圖脸哀，并輸出plot查看
ggsave(p,filename = "gene.pathway.pdf"),width = 8,height = 7) #保存圖片

9. 差異表達(dá)基因火山圖

通過ggplot2軟件繪制基因差異表達(dá)火山圖蹦浦，使用DESeq2軟件生成的matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results為輸入文件，該文件內(nèi)geneid撞蜂、sampleA盲镶、sampleB、baseMeanA蝌诡、baseMeanB溉贿、baseMean、log2FoldChange浦旱、lfcSE宇色、stat、pvalue颁湖、padj等信息宣蠕，將其導(dǎo)入R語言中，提取1,6,7,8,9,10,11列甥捺，然后設(shè)置數(shù)據(jù)的行列名稱

library(ggplot2)  # 載入R包
res <- read.table("matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results")
res <- res[,c(1,6,7,8,9,10,11)]           
#篩選數(shù)據(jù)抢蚀，提取指定列
rownames(res) <- res$V1                    
#命名數(shù)據(jù)框的行名為基因ID
colnames(res) <- c("geneid","baseMean",   
"log2FoldChange","lfcSE","stat","pvalue","padj")  # 修改列名
deseq_res <- as.data.frame(res[order(res$padj), ]) # 按照p值大小排序
write.table(deseq_res, 'DESeq2-test.txt', row.names = FALSE, 
sep = '\t', quote = FALSE)             
#保存差異表達(dá)基因數(shù)據(jù)

設(shè)置差異表達(dá)分類sig，最后利用ggplot函數(shù)繪制火山圖镰禾，保存為'volcano.svg'

deseq_res <- read.delim('DESeq2-test.txt', sep = '\t')
deseq_res[which(deseq_res$log2FoldChange >= 1 & deseq_res$padj <
 0.05),'sig'] <- 'up (p.adj < 0.05, log2FC >= 1)'        
#上調(diào)基因
deseq_res[which(deseq_res$log2FoldChange <= -1 & deseq_res$padj <
     0.05),'sig'] <- 'down (p.adj < 0.05, log2FC <= -1)'     
#下調(diào)基因
deseq_res[which(abs(deseq_res$log2FoldChange) < 1 | deseq_res$padj >= 0.05),'sig'] <- 'no diff'                 
#無顯著差異
p <- ggplot(deseq_res, aes(log2FoldChange, -log(padj, 10))) +
  geom_point(aes(color = sig), alpha = 0.6, size = 1) +
  scale_color_manual(values = c('blue2', 'gray30', 'red2')) +
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(color = 'black') 
        legend.position = c(0.26, 0.92)) +
  theme(legend.title = element_blank(), 
        legend.key = element_rect(fill = 'transparent'), 
        legend.background = element_rect(fill = 'transparent')) +
  geom_vline(xintercept = c(-1, 1), color = 'gray', size = 0.25) +
  geom_hline(yintercept = -log(0.05, 10), color = 'gray') +
  labs(x = 'log2 Fold Change', y = '-log10 p-value', color = NA) +
  xlim(-5, 5)
ggsave('volcano.svg', p, width = 4, height = 3)     
# 保存火山圖

參考資料：
[1] Brown J, Pirrung M, McCue L A. FQC Dashboard: integrates FastQC results into a web-based, interactive, and extensible FASTQ quality control tool[J]. Bioinformatics, 2017, 33(19): 3137-3139.
[2] Bolger A M, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data[J]. Bioinformatics, 2014, 30(15): 2114-2120.
[3] Kim D, Paggi J M, Park C, et al. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype[J]. Nature biotechnology, 2019, 37(8): 907-915.
[4] Liao Y, Smyth G K, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features[J]. Bioinformatics, 2014, 30(7): 923-930.
[5] Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[J]. Genome biology, 2014, 15(12): 1-21.
[6] Yu G, Wang L G, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters[J]. Omics: a journal of integrative biology, 2012, 16(5): 284-287.
[7] Wickham H. Data analysis[M]//ggplot2. Springer, Cham, 2016: 189-201.

本文由mdnice多平臺(tái)發(fā)布