RNA-seq转录组数据分析丨一套完整的案例流程

2022-12-11 本文已影响0人生信分析笔记

今天分享的学习笔记是一套转录组分析简单流程，适用于初学者入门阅读，从原始测序数据开始，经过质控、序列比对、定量表达、差异表达、功能富集等一系列分析步骤，最终获得基因表达信息。本文所有数据都经过特殊修改，仅供学习参考使用。

转录组是在特定时空条件下细胞中基因转录表达产物，广义的转录组包括信使RNA，核糖体RNA，转运RNA及非编码RNA，狭义上是指所有mRNA的集合，转录组分析能够获得不同基因的表达情况。

1. 数据来源

假设有两个不同组织（PR和SR），每个组织各区三个样本，一共六个样本，利用illumina平台进行转录组测序，得到双端测序数据。数据原始格式为.fq,共有12条测序数据文件（每个样本产生两条）

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

创建工作文件夹，并新建子步骤目录：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   
#批量创建工作文件夹

2. 测序数据质量评估

利用fastQC软件对获得的fastq序列文件进行质量分析，生成html格式的结果报告，其中含有各项指标，以下用PR1样本为例。

fastqc *.fq        #批量对fq后缀文件运行fastqc程序
输出结果：PR1_1_fastqc.html  
Filename    PR1_1.fq  
File type   Conventional base calls
Encoding    Illumina 1.5
Total Sequences 105300       #总序列数
Sequences flagged as poor quality   0
Sequence length 90              #序列长度
%GC 52                #GC碱基含量

质量评估报告，使用浏览器打开：

3. 过滤低质量序列

利用Trimmomatic软件除去序列文件中的接头（adapter），并对碱基进行合适的修改，然后对碱基进行修剪，对低质量的序列进行过滤。

time java -jar trimmomatic-0.39.jar PE  #trimmomatic是依赖java环境运行
-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

运行程序对序列文件中低质量的序列进行过滤，将输出结果存储到01_trimmomaticFiltering，每一个样本序列会生成如下输出文件，summary文件包含过滤的总结信息。

-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

打开其中一个PR1.summary文件进行查看，其中Input Read Pairs表示过滤之前的reads数，Both Surviving Reads表示过滤之后的reads数。

$ cat PR1.summary  #查看总结文件
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. 比对到参考基因组

利用hisat2软件，将fasta序列比对到参考基因组。首先需要构建索引文件，下载或者拷贝参考基因组序列文件Ref和基因组注释文件，通过hisat2-build命令构建索引文件。

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比对到参考基因组上。输入文件为两个fasta序列文件，将运行过程中的输出提示重定向到log文件。

-S设置输出文件为sam
-x设置参考基因组
-p设置运行的线程数量

time hisat2 -p 1 #线程数为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  
#输出日志

得到的日志文件中包含比对成功的reads数量和比对率等信息：

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%   #总比对率8.76%

上述步骤完成后，会在当前目录下生成多个sam格式文件，SAM（sequence Alignment/mapping)格式是高通量测序中存放比对数据的标准格式。另外，bam是sam的二进制格式，可以减小sam文件的大小，因此利用samtools对sam进一步处理，得到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. 基因表达定量分析

利用featuresCounts软件，对reads进行精确的计数，最后将所有样本的reads数合并为一个文件，将RNA-Seq_Practice_countstable文件导出，根据FPKM和TPM的计算公式定量每个基因在每个样本中的表达量。利用网站（https://www.ncbi.nlm.nih.gov/gene）将基因ID转化为Gene description，从而了解其功能和相关信息。

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数，得到count文件，然后同样步骤处理其他样本，最后将所有的reads定量数据合并为一个countstable文件。

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

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

将生成的RNA-Seq_Practice_countstable保存到本地，然后计算FPKM和TPM值，在R语言中进行相关计算。

FPKM（Fragments Per Kilobase of exon model per Million mapped fragments）表示每千个碱基的转录每百万映射读取的fragments，该方法是利用每个样本的总fragments数进行校正。优点是消除样本间和基因间的差异，但是该方法主要是计算的结果可能与真实表达水平有偏差。

TPM (transcript per million) 表示每百万转录本中来自于某基因的转录本数目，该方法先对每个基因的reads数用基因的长度进行校正，然后再利用校正后的reads数与校正后该样本所有的reads数求商。优点是消除了外显子长度造成的差异和样本间测序总reads counts不同造成的差异，缺点为该法不是采用比对到基因组上的总reads counts，所以有时候不够精确。

具体计算方法如下所述：

> df <- read.table("RNA-Seq_Practice_countstable") #读入数据文件
> rownames(df) <- df$V1   #修改数据表的行名为基因ID（第一列V1）
> df <- df[,c(2,3,4,5,6,7)] #现在行名已经为基因ID，因此删除数据的第一列
> colnames(df) <- c("PR1","PR2","PR3","SR1","SR2","SR3") #修改行名
> dim(df); names(df)   #查看df的数据结构和变量名称是否正确
[1] 41768     6   
[1] "PR1"    "PR2"    "PR3"    "SR1"    "SR2"    "SR3"   
> head(df[,1:6],2)      #输出前两行数据进行预览
             PR1 PR2 PR3 SR1 SR2  #PR和SR共六列数据
Zm00001d027   0   0   0   0   0
Zm00001d021   0   0   0   0   0

featureCounts_meta <- df[,1:6]  #提取基因信息
df <- df[rowSums(df)>0,]     #去除表达量均为0的行
prefix <-"maize_PR_SR"   #设置输出文件前缀名

# ----- TPM计算 ------
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计算 ------
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. 基因差异表达分析

利用DESeq2软件基于样本的原始reads数目计算差异表达基因，利用run-featurecounts.R脚本对每个样本的reads数进行定量，然后利用abundance_estimates_to_matrix.pl脚本合并所有的定量结果，最后利用DESeq2进行差异表达基因的分析。

首先，将所需的脚本文件全部拷贝到工作目录下，同时将样本信息相关文件也拷贝保存。

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

新建一个文件夹，运行R语言脚本，对样本进行表达定量分析，以PR1为例，输入文件为02文件夹中生成的bam文件和基因组注释文件，生成输出文件PR1，其他的几个样本按照同样的方法进行处理。

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

运行结束后，在输出文件夹内将会生成count文件和log文件，利用perl语言脚本abundance_estimates_to_matrix.pl合并所有样本的定量结果。

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

利用trinity 软件中的DESeq2分析差异表达基因，输入文件有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   #这里需要注意文件名称是否正确

上述流程输出的结果存放在随机文件夹中，进入该文件夹后，可以看到生成的差异表达基因结果matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results文件，然后利用awk命令提取出满足条件（同时满足：|log2FoldChange| > 1，且padj < 0.05）的差异表达基因,生成新的输出文件件PRvsSR_DEGs。

drwxr-xr-x 307 Nov 29 09:54 DESeq2.197264.dir  #随机新目录
$ ls -f  #查看新目录下文件列表
matrix.counts.matrix.PR_vs_SR.DESeq2.Rscript
matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results  #差异表达基因结果文件
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命令提取满足条件的指定数据

将最终得出的差异表达基因数据文件PRvsSR.DESeq2_DEGs保存到电脑，进行GO富集和KEGG富集分析。绘制火山图时使用的是matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results

7. GO富集分析

利用在线网站(http://shiny.host.ugeneyun.cn/GOmaize/)进行GO富集分析，输入文件为差异表达基因的ID信息表，此处不进行显著性过滤，原因是前期已经进行筛选了。提交运行后生成GOenrich结果，共有270条Term，下图随机展示5条，可以得出GO编号、功能描述、显著性和分类信息。然后，利用网站的绘图功能制作GO富集气泡图，保存到本地GO.pdf。

8. KEGG分析

在R语言（4.2.2）中使用clusterProfiler、KEGG.db、ggplot2三个R包进行KEGG富集分析，先构建一个OrgDb本地索引库，然后利用enrichKEGG函数进行富集分析，利用在线网站(https://www.maizegdb.org/gene_center/gene)将基因ID转换为GenBank Gene类型，保存为文件gene.txt，输出文件为pathway数据表gene.pathway.csv 和KEGG富集结果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)   # 该步骤需要手动在Rstudio中安装KEGG.db包

file="gene.txt"  #输入文件，包含18条差异表达的基因信息
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)  #进行KEGG富集分析
write.csv(kk,file = paste0("gene.pathway.csv"),row.names = F)
p=dotplot(kk);p  #绘制点图，并输出plot查看
ggsave(p,filename = "gene.pathway.pdf"),width = 8,height = 7) #保存图片

9. 差异表达基因火山图

通过ggplot2软件绘制基因差异表达火山图，使用DESeq2软件生成的matrix.counts.matrix.PR_vs_SR.DESeq2.DE_results为输入文件，该文件内geneid、sampleA、sampleB、baseMeanA、baseMeanB、baseMean、log2FoldChange、lfcSE、stat、pvalue、padj等信息，将其导入R语言中，提取1,6,7,8,9,10,11列，然后设置数据的行列名称

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)]           
#筛选数据，提取指定列
rownames(res) <- res$V1                    
#命名数据框的行名为基因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)             
#保存差异表达基因数据

设置差异表达分类sig，最后利用ggplot函数绘制火山图，保存为'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)'        
#上调基因
deseq_res[which(deseq_res$log2FoldChange <= -1 & deseq_res$padj <
     0.05),'sig'] <- 'down (p.adj < 0.05, log2FC <= -1)'     
#下调基因
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多平台发布