RNA-seq分析RNA seq好知识

RNAseq常规流程

2021-05-21  本文已影响0人  QXPLUS

转录组分析是生信分析的一个基础和常见的分析流程，一般从下机数据开始，经过一系列的质控，基因组比对，差异分析等过程得到实验组与对照组（或者肿瘤中的转移组和原癌组）中差异表达的基因集，然后在基于该基因集结合我们的研究方向进行进一步的功能分析，信号通路分析，细胞间通讯，实验验证等得到一个可解释的实验结果。这篇文章主要就涉及到的基本流程和相应软件进行一个简单的介绍，后续更深入的知识点还需要小伙伴们自己去挖掘哦。

基本分析流程：

  1. 原始数据 -- Rawdata
  2. 质控 -- QC -- Trimmomatic，FastQC
  3. 基因组比对 -- Aligment -- STAR
  4. 比对结果量化 -- FeatureCount -- featurecount
  5. 两组间的基因表达差异分析 -- DEGs -- limma
  6. 差异基因的功能富集分析 -- GO,KEGG -- clusterProfiler

一. 测序数据下机后处理

首先，我们拿到的原始序列文件就是fq.gz 文件

1.下载fastq.gz rawdata

wget url

2.每个样本创建一个以样本名命名的文件夹（eg RNA-1）, 并将相应的._1.fq.gz/._2.fq.gz 放入文件夹（eg RNA-1），用作后续QC处理。

mkdir sample1 sample2 sample3 control1 control2 control3
for samp in {sample1 sample2 sample3 control1 control2 control3}
  do
  mv "../rawdata/*"$samp"*_1.fq.gz" $samp 
  mv "../rawdata/*"$samp"*_2.fq.gz" $samp 
  done

二. 低质量数据过滤和质量评估QC -- Trimmomatic/FastQC

1. Trimmomatic

Trimmomatic: A flexible read trimming tool for Illumina NGS data

### 切除接头序列
java -jar /home/Software/Trimmomatic-0.39/trimmomatic-0.39.jar  \
PE  \    #  rawdata为双端测序
-phred33 \   #  Fastq文件的质量值格式
-threads 6 \   #  6 线程
-trimlog Trimmomatic.log \   # 日志文件保存位置
$read1 $read2 \     # *_1.fq.gz, *_2.fq.gz  需要过滤的Fastq文件
$trim_read1 \        # *_1.clean.fq.gz  过滤后的Fastq文件
output_unpaired_1.fq.gz \    #  read1的5‘->3'和read2的3‘->5'未能匹配上的Fastq文件
$trim_read2 \        
output_unpaired_2.fq.gz \
ILLUMINACLIP:/home/Software/Trimmomatic-0.39/adapters/TruSeq3-PE-2.fa:2:30:10:1:TRUE  \   # 去除illumina测序平台下的TruSeq3接头序列，不知道可以问测序公司。
# 2:30:10即表示，在比对接头序列时允许有两个位置的碱基发生错配，双端测序的两条reads与接头序列匹配率超过30%的话，就会被切除掉，单条reads如果与接头序列的匹配率超过10%，也会被切除掉。
LEADING:0 \    # 切除reads 5’端质量值低于0的碱基, 即不对5’端剪切
TRAILING:25 \    # 切除reads 3’端质量值低于25的碱基
SLIDINGWINDOW:50:25 \  #以50bp为窗口进行滑窗统计，切除碱基平均质量低于25的窗口及之后的序列。
MINLEN:50 \   #去除过滤后长度低于50的reads
HEADCROP:0 \   # 切除reads开头的0个碱基
2>trim.log 1>trim.err   # 0 表示stdin标准输入；1 表示stdout标准输出；2 表示stderr标准错误

使用参数

2.FastQC (测序数据质控)

结果解读
我们对Trim过滤之后的序列进行质量评估，评估结果如下：
FastQC 会生成一个fastqc_report.html的结果报告，下面是软件对质控结果进行判断：绿色代表PASS；黄色代表WARN；红色代表FAIL

image.png

简单说一下报告中常见的几张图吧

横轴为read长度，纵轴为质量得分。说明trim后的Reads质量较高，可以用于后续分析。
一般要求所有位置的10%分位数大于20，即大于最多允许该位置10%的序列低于Q20。绿色区域说明数据质量较高。一般但随着测序长度的加深，碱基的错赔率会增加。

横轴表示Q值，纵轴表示每个值对应的read数目，当测序结果主要集中在高分中，证明测序质量良好

GC含量-GC含量对应的read数。蓝色为经验分布给出的理论值，红色是测序数据的真实值，但两条线接近重合时，说明数据质量较高；当红色出现双峰可能是样品被污染或不纯。

其他更多信息可以去官网自行查看。

三. Clean reads的基因组比对 Aligment -- STAR

针对每个样本，利用 STAR软件 将预处理序列与测序物种的参考基因组序列进行序列比对。

1. 比对原理

将QC后的reads比对到参考基因组上，一般分为有参考基因组比对和无参考基因组比对。我们常用的就是有参考基因组的比对。


aligment.png

在参考基因组比对的时候，会优先考虑大片段比对，如果比对不上，再考虑将片段切割成不同区段进行比对。

RNAseq 常用的比对工具有：Hisat2, STAR 等

2. STAR 比对

STAR 
--runThreadN 20 \      # 线程数
--genomeDir /home/database/ref \   # 参考基因组
--readFilesCommand gunzip -c \  # *.gz 文件解压缩
--readFilesIn *_1.fq.gz *_2.fq.gz \  #输入文件，空格隔开
--outSAMtype BAM SortedByCoordinate \  #  输出bam文件（默认输出sam文件）
--outFileNamePrefix Aligment \  # 输出文件的位置和前缀

结果解读

BAM/SAM文件
bam是sam的压缩格式，在内容上是一样的，linux下可使用samtools, bamtools 进行查看与转换。
bam文件分为2个部分，header和record.

比对率统计


mapcount.png

总比对率和唯一比对率是我们关注的，尤其是总比对率。

比对区域分布统计


image.png

比对reads在每条染色体上的密度分布情况


image.png

可以通过IGV软件可视化查看BAM文件信息
具体以后再详细展开……

四. 比对结果的外显子（基因）表达计数 FeatureCount -- featurecount

在高通量测序分析中用于下游分析的关键信息是比对到每个genomic feature（外显子的表达量、featurecounts是一款使用于RNA-seq和DNA-seq的read summarization工具，应用了高效率的染色体哈希算法和feature区块技术因等）中的read数目。
featurecounts 速度快，需要的内存空间少，同时可以用于单端和双端的数据

通过STAR等软件比对得到的BAM文件，利用计数软件（如 featurecount, HTseq, stringtie），获取各基因的 reads count。
假设前提：认为比对到基因的reads 数与该基因的表达量成正相关。


readscount.png

featurecount的调用方法

/home/Software/featureCounts \
-T 4 \    # 线程数目，1~32
-a $gtf \  # 参考gtf文件
-o read.count \   # 输出文件的名字，输出文件的内容为read 的统计数目
-p \     # 只能用在paired-end的情况中，会统计fragment而不统计read  (-p -B -C 同时使用、不使用)
-B \    # 在-p选择的条件下，只有两端read都比对上的fragment才会被统计
-C \    # 如果-C被设置，那融合的fragment就不会被计数，这个只有在-p被设置的条件下使用
-t exon \   # 设置feature-type，-t指定的必须是gtf中有的feature，同时read只有落到这些feature上才会被统计到，默认是“exon”
-g gene_name \  # 默认为gene_id，注意！选择gtf中提供的id identifier！！！
-s 2  \   # 2条链比对
Align.sorted.bam   # 输入的bam/sam文件，支持多个文件输入

得到read.count数据之后，还需要进行进一步处理才能够进行分析。

由于实验中每个样本的上样量不同，或者基因的长度,测序深度不同等导致基因count本身就存在着较大差异，所以，我们是不能够直接只用count 数据进行后续的基因表达差异分析的。因此我们需要对 count 数据进行标准化。

countToFpkm <- function(counts, effLen) {N <- sum(counts)
     exp( log(counts) + log(1e9) - log(effLen) - log(N) )}

转录本表达量计算采用FPKM计算度量指标(FPKM- Fragments Per Kilobase of transcript per Million fragments mapped)。FPKM含义是以每百万比配成对序列每1Kbp长度做转录本表达量指标，其中转录本长度和总比配成对read数目用于归一化表达量数值。用cuffnorm程序分别计算各个样本FPKM表达量取log2。

genecount.png

五. 不同分组间的基因差异表达分析 DEGs -- limma

根据样本的实验设计，利用limma包对不同样本组之间筛选差异表达的已知转录本。满足pvalue <= 0.05和大于等于2倍差异表达范围筛选两组之间的差异转录本


image.png

差异基因可视化 -- 火山图


volcano.png wenn.png

针对组间筛选的不同差异基因集， 采用venn方法进行共同和特有的比较分析

heatmap.png

针对样本组间筛选的差异表达基因， 采用对基因和样本进行双向层级聚类并且用热图显示

六. 差异基因的功能富集分析 -- GO,KEGG

拿到差异基因之后，需要将这些基因从微观映射到宏观，从分子表达水平上升到功能分析，这样也具有更好的解释性。GO和KEGG就是基于不同的分类思想而储存的基因相关功能的数据库。

1. GO功能分析

GO功能分析是针对全基因/转录本和差异基因/转录本进行功能注释和归类。GO数据库，全称是Gene Ontology(基因本体)，他们把基因的功能分成了三个部分分别是：细胞组分（cellular component, CC）、分子功能（molecular function, MF）、生物过程（biological process, BP）。利用GO数据库，我们就可以得到我们的目标基因在CC, MF和BP三个层面上，主要和什么有关。

library(clusterProfiler)
ego<-enrichGO(gene=genelist,OrgDb="org.Hs.eg.db",minGSSize = 2,keyType = "ENTREZID",ont ="ALL",pvalueCutoff = 0.1, readable=TRUE)
write.table(ego@result,"GO_enrichment.xls",row.names=F,quote=F,sep="\t")

条形图可以显示，DEGs中的差异基因都显著富集到了哪些通路上面。


GO.png

还可以更详细地查看通路富集情况


GO-tree.png

2.KEGG富集通路分析

KEGG数据库：除了对基因本身功能的注释，我们也知道基因会参与人体的各个通路，基于人体通路而形成的数据库就是通路相关的数据库。而KEGG就是通路相关的数据库的一种。其实通路数据库有很多，类似于wikipathway，reactome都是相关的通路数据库。只是因为KEGG比较被人熟知，所以基本上都做这个分析的.

library(clusterProfiler)
kk<-enrichKEGG(gene=genelist,organism='hsa',keyType="ncbi-geneid",pvalueCutoff = 0.05,pAdjustMethod = "BH",  minGSSize = 2, maxGSSize = 500,qvalueCutoff = 0.2, use_internal_data = FALSE)
res_pathview<-kk@result
kk<-setReadable(kk,OrgDb="org.Hs.eg.db",keyType="ENTREZID") ###mapping gene id to gene symbol
res<-kk@result
colnames(res)[1]<-"pathwayID"
write.table(res,"KEGG_enrichment.xls",quote=F,row.names=F,sep="\t")

富集通路的条形图，可以显示差异基因富集到的通路有哪些。


keggpathway.png
上一篇 下一篇

猜你喜欢

热点阅读