RNA检测流程

风风是超人

已于 2022-08-09 14:12:44 修改

阅读量1.8k

点赞数 1

分类专栏： pipline 文章标签：经验分享

于 2022-04-19 14:44:36 首次发布

本文链接：https://blog.csdn.net/qq_42962326/article/details/124268133

版权

pipline 专栏收录该内容

1 篇文章 1 订阅

订阅专栏

1 比对

RNA数据的比对软件，和DNA数据的比对软件有何区别？
由于转录过程中存在可变剪接，测得的reads有很大一部分跨越不同外显子，所有RNA数据的比对软件要支持reads的splice比对。

Hist2 下载
https://github.com/DaehwanKimLab/hisat2

./hisat2/hisat2-2.2.1/hisat2 \
 -x /hg19index/ucsc.hg19 \
 -1 fq1.gz \
 -2 fq2.gz \
 -S out.sam \
 --dta --rna-strandness

-x ：对hg19参考基因组做的索引

cd  hg19index
hisat2/hisat2-2.2.1/hisat2-build  -p  4  ucsc.hg19.fasta   ucsc.hg19

2 sam转bam

samtools view -bS out.sam > out.bam
samtools sort out.bam out.sorted

3 定量分析

比对完成后，根据基因在参考基因组上的位置信息，便可统计每个基因从起始到终止范围内覆盖的reads数目。采用subread软件中的featureCounts工具对基因进行表达水平的定量，featureCounts主要使用-Q 10 -B -C参数，分别过滤掉比对质量值低于10的reads，非成对比对上的reads，比对到基因组多个区域的reads。

featureCounts利用高效的染色体拆分和特征模块的技术，先对染色体拆分成bin，在bin中拆分blocks，在blocks中拆分feature。它有两种重要的特征值：features和meta-features。其中features通常是对应于基因注释文件(GTF)的exon，meta-features通常是对应于注释文件的gene。在每个染色体上的features根据起始位置进行排序，将染色体分成大小为128kb的没有重叠的bins，根据feaatures的起始位置分配到对应的bin中，其中图中的实线表示的features。在每一个bin中连续的features中数量相等的分类成相同的模块(Block)，其中在每一个bin中Block的数量是与在这个bin中features的数目的平方根相等。图中搜索的reads首先与基因组的bins进行比较，然后再与bin中的blocks比较，最后与blocks中的features进行比较，最终的read的定位结果是定位到了第i个bin中的第一个blocks中，且read与blocks中的两个features有重叠，即junction read。
在这里插入图片描述

3.1 featureCounts 进行基因表达定量(read count)

3.1.1 下载 subread

https://sourceforge.net/projects/subread/

3.1.2 下载GTF注释文件

ensembl官网
http://grch37.ensembl.org/index.html
在这里插入图片描述

在这里插入图片描述

wget http://ftp.ensembl.org/pub/grch37/current/gtf/homo_sapiens/Homo_sapiens.GRCh37.87.gtf.gz

gtf：
在这里插入图片描述

3.1.3 分析 featureCounts

featureCounts/subread-1.5.0-p3-Linux-x86_64/bin/featureCounts \
 -a Homo_sapiens.GRCh37.87.gtf \
 -o sample.featureCounts.txt \
 -Q 10 -B -C \
 input.sorted.bam

在featureCounts 软件中，有两个核心概念：

feature
metafeature

feature指的是基因组区间的最小单位，比如exon;
而metafeature可以看做是许多的feature构成的区间，比如属于同一个gene的外显子的组合。

-a参数指定的区间注释文件，默认是gtf格式;-T参数指定线程数，默认是1;-t参数指定想要统计的feature的名称,取值范围是gtf 文件中的第3列的值，默认是exon;-g参数指定想要统计的meta-feature的名称，取值范围参考gtf第9列注释信息，gtf的第9列为key=value的格式，-g参数可能的取值就是所有的key, 默认值是gene_id。
在这里插入图片描述

4. stringTie 进行FPKM定量

安装stringTie
https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual#input

./stringTie/stringtie-2.2.1.Linux_x86_64/stringtie  \
 -G Homo_sapiens.GRCh37.87.gtf \
 -e -b ballgown_out_dir\
 -o out.gtf \
 input.bam

在输出的 out.gtf 和 ballgown_out_dir 目录中有FPKM结果

或者使用 rsem-calculate-expression 软件，可以同时进行readcount计算，FPKM,TPM 计算。

5. star-fusion 检测融合

融合基因是指两个基因的全部或部分序列融合而成的嵌合基因，一般由染色体易位、缺失等原因所致。融合基因首次发现于血液系统的恶性肿瘤中，其中以慢性粒细胞白血病中BCR-ABL的基因融合最为经典，治疗慢性粒细胞白血病的药物伊马替尼/格列卫，其作用靶点就是该融合基因。高通量RNA测序技术因其通量高、成本低、检测精度高和检测范围广等优点大大加快了融合基因的研究，诺禾致源使用STAR-Fusion软件进行融合基因的检测，STAR-Fusion是利用STAR比对的融合输出结果来检测融合转录本的软件包，STAR-Fusion分析流程如下图所示。分为SATR比对，STAR-Fusion.predict，STAR-Fusion.filter，具体如下。

STAR比对：是STAR比对的标准比对过程中的一个延伸，第一步，利用MMP算法（连续最大可比对种子搜索）将种子准确的定位到参考基因组上，下一步，选择基因组比对窗口将固定的种子进行聚类。

STAR-Fusion.predict：将Split reads又称JunctionReads(含有两个基因融合断点的reads，流程图中S=3即为Split reads的数目)和Discordant pair又称SpanningFrags(不一致比对的reads，即reads的两端比对到不同的两个基因，图中J=2即为Discordant pair的数目）比对到参考基因组的注释文件。STAR-Fusion根据支持融合断点的Split reads和Discordant pair的数目应用最少reads支持准则对融合基因进行筛选。

STAR-Fusion.filter：这一步是融合基因检测的最后一步，主要是筛选出最可靠的融合基因，过滤到最初预测的不可靠的候选的融合基因。过滤的过程主要是先按照基因对和断点的距离进行分类，然后按照在推断的断点reads的支持和比对的范围以及推断的融合基因的序列的相似性，过滤掉配对混乱的融合基因对。

FusionInspector 是一个癌症转录分析工具包，软件会对STAR-Fusion的预测结果通过再检测，再次打分，进一步对融合基因的分析结果进行校正分析。具体方法：输入的文件是我们用STAR-Fusion得到的融合基因的列表，参考基因组文件，样品的clean data。FusionInspector会根据STAR-Fusion得到的融合基因，对每对融合基因提取在基因组上的区域并构建融合基因的最小融合的contig序列。将reads比对到这些候选的融合的contig上。最后对每个融合基因支持的JunctionReads与SpanningFrags进行识别。得到最终的融合基因列表。

STAR-Fusion官网推荐对STAR-Fusion的预测结果再次用校验工具FusionInspector进行进一步的矫正分析。通过两次的矫正保证了融合基因结果的准确性。我们融合基因分析的过程就是利用FusionInspector对STAR-Fusion预测结果进行校正。
在这里插入图片描述

5.1 下载 STAR

网址：https://github.com/alexdobin/STAR/

wget https://github.com/alexdobin/STAR/archive/refs/tags/2.7.9a.tar.gz

5.2 下载 star-fusion

网址：https://github.com/STAR-Fusion/STAR-Fusion
下载带 FULL 的压缩包

wget https://github.com/STAR-Fusion/STAR-Fusion/releases/download/v1.10.0/STAR-Fusion-v1.10.0.FULL.tar.gz

5.3 下载star-fusion依赖的资源包

wget -c https://data.broadinstitute.org/Trinity/CTAT_RESOURCE_LIB/GRCh37_gencode_v19_CTAT_lib_Mar012021.plug-n-play.tar.gz

5.4 分析fusion

从fastq开始，通过star-fusion 内嵌的STAR进行比对，然后用star-fusion检测融合

/STAR-Fusion-v1.10.0/STAR-Fusion \
 --genome_lib_dir 下载的资源包路径/ctat_genome_lib_build_dir \
 --STAR_PATH 下载的STAR/STAR \
 --FusionInspector inspect \
 --left_fq fq1.gz \
 --right_fq fq2.gz \
 --output_dir outdir

只利用 STAR的结果分析

/STARfusion/STAR-Fusion-v1.10.0/STAR-Fusion \
 --genome_lib_dir ./ctat_genome_lib_build_dir \
 --STAR_PATH ./path to/STAR \
 -J Chimeric.out.junction \  # STAR的输出结果
 --output_dir outdir

最终输出的 star-fusion.fusion_predictions.tsv 为检测结果

6 差异基因分析（DESeq2或edgeR）

DESeq2 进行分析
转载：https://www.sohu.com/a/400667791_120736615

整理read count数据，只能为正整数，不能为小数或者负数。分为对照组和实验组。
在这里插入图片描述

.libPaths("path_to/software/conda_py38/conda3/envs/r3.6/lib/R/library")
library(DESeq2)

dat <- read.delim('deseq2TestDATA.xls', row.names = 1, sep = '\t', check.names = FALSE)
coldata <- data.frame(condition = factor(rep(c('control', 'treat'), each = 3), levels = c('control', 'treat')))

##DESeq2 默认流程
#第一步，构建 DESeqDataSet 对象

dds <- DESeqDataSetFromMatrix(countData = dat, colData = coldata, design= ~condition)

#第二步，计算差异倍数并获得 p 值

#备注：parallel = TRUE 可以多线程运行，在数据量较大时建议开启
dds1 <- DESeq(dds, fitType = 'mean', minReplicatesForReplace = 7, parallel = FALSE)

#注意，需将 treat 在前，control 在后，意为 treat 相较于 control 中哪些基因上调/下调

res <- results(dds1, contrast = c('condition', 'treat', 'control'))

#输出表格至本地
#差异分析结果保存在“res1”中，包含了基因id、标准化后的基因表达值、log2转化后的差异倍数（Fold Change）值、显著性p值以及校正后p值padj(adjusted p values)（默认FDR校正）等主要信息
res1 <- data.frame(res, stringsAsFactors = FALSE, check.names = FALSE)

write.table(res1, 'control_treat.DESeq2.txt', col.names = NA, sep = '\t', quote = FALSE)


##筛选差异表达基因

#首先对表格排个序，按 padj 值升序排序，相同 padj 值下继续按 log2FC 降序排序

res1 <- res1[order(res1$padj, res1$log2FoldChange, decreasing = c(FALSE, TRUE)), ]

#log2FC≥1 & padj<0.01 标识 up，代表显著上调的基因

#log2FC≤-1 & padj<0.01 标识 down，代表显著下调的基因

#其余标识 none，代表非差异的基因

res1[which(res1$log2FoldChange >= 1 & res1$padj < 0.01),'sig'] <- 'up'

res1[which(res1$log2FoldChange <= -1 & res1$padj < 0.01),'sig'] <- 'down'

res1[which(abs(res1$log2FoldChange) <= 1 | res1$padj >= 0.01),'sig'] <- 'none'

#输出选择的差异基因总表

res1_select <- subset(res1, sig %in% c('up', 'down'))

write.table(res1_select, file = 'control_treat.DESeq2.select.txt', sep = '\t', col.names = NA, quote = FALSE)

#根据 up 和 down 分开输出

res1_up <- subset(res1, sig == 'up')

res1_down <- subset(res1, sig == 'down')

write.table(res1_up, file = 'control_treat.DESeq2.up.txt', sep = '\t', col.names = NA, quote = FALSE)

write.table(res1_down, file = 'control_treat.DESeq2.down.txt', sep = '\t', col.names = NA, quote = FALSE)


##ggplot2 差异火山图

library(ggplot2)

#默认情况下，横轴展示 log2FoldChange，纵轴展示 -log10 转化后的 padj

p <- ggplot(data = res1, aes(x = log2FoldChange, y = -log10(padj), color = sig)) +
geom_point(size = 1) + #绘制散点图
scale_color_manual(values = c('red', 'gray', 'green'), limits = c('up', 'none', 'down')) + #自定义点的颜色
labs(x = 'log2 Fold Change', y = '-log10 adjust p-value', title = 'treat vs control', color = '') + #坐标轴标题
theme(plot.title = element_text(hjust = 0.5, size = 14), panel.grid = element_blank(), #背景色、网格线、图例等主题修改
panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = c(-1, 1), lty = 3, color = 'black') + #添加阈值线
geom_hline(yintercept = 2, lty = 3, color = 'black') +
xlim(-12, 12) + ylim(0, 35) #定义刻度边界

png(file="huoshan.png", bg="transparent")
p
dev.off()

在这里插入图片描述

6.2 使用edgeR进行无重复差异表达分析

转载：http://events.jianshu.io/p/c1036d39f2b9

我们需要安装两个R包，一个是edgeR, 一个是airway. 其中airway是一个数据集包, 功能就是提供一个用于分析的数据

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

if (!requireNamespace("edgeR", quietly = TRUE))
    BiocManager::install("edgeR")

if (!requireNamespace("airway", quietly = TRUE))
    BiocManager::install("airway")

加载R包

library(edgeR)
library(airway)

在这里插入图片描述

我个人是不太推荐没有重复的差异表达分析，毕竟统计学上的p值是为了证明两个样本的差异是真实存在而不是抽样误差导致, 但是你单个样本如何计算变异呢?

因此每当别人提问的时候, 我个人的建议就是定性看看倍数变化吧. 但是如果真的强行要算p值, 其实也不是不行, edgeR就是一种选择.

差异表达分析
不同差异表达分析工具的目标就是预测出dispersion(离散值), 有了离散值就能够计算p值. 那么dispersion怎么计算呢? edgeR给了几个方法

方法一:

根据经验给定一个值(BCV, square-root-dispersion). edgeR给的建议是, 如果你是人类数据, 且实验做的很好(无过多的其他因素影响), 设置为0.4, 如果是遗传上相似的模式物种, 设置为0.1, 如果是技术重复, 那么设置为0.01

这里用的数据是人类, 此处设置为 0.4

y_bcv <- y
bcv <- 0.4
et <- exactTest(y_bcv, dispersion = bcv ^ 2)

我们用decideTestsDGE看下有多少基因上调, 多少基因下调. 设置p.value=0.05

gene1 <- decideTestsDGE(et, p.value = 0.05, lfc = 0)
summary(gene1)

         2-1
Down       4
NotSig 14733
Up        19

差异基因少的可怜, 只有4+19个。我们可以尝试调整下BCV

y_bcv <- y
bcv <- 0.2
et2 <- exactTest(y_bcv, dispersion = bcv ^ 2)
gene2 <- decideTestsDGE(et2, p.value = 0.05, lfc = 0)
summary(gene2)

         2-1
Down     159
NotSig 14380
Up       217

这个时候的差异基因上升到了159 + 217个.

和有重复表达进行比较

由于我们的数据集原来是有重复的，那么我们就可以比较下无重复和有重复之间会相差多少基因

group <- meta_info$dex
y_rep <- DGEList(counts=expr_matrix, group = group)
keep <- rowSums(cpm(y_rep)>1) >= 5
y_rep <- y_rep[keep, , keep.lib.sizes=FALSE]
y_rep <- calcNormFactors(y_rep)
design <- model.matrix(~group)
y_rep <- estimateDisp(y_rep, design = design)
fit <- glmQLFit(y_rep, design)
qlf.2vs1 <- glmQLFTest(fit, coef=2)

我们看下差异基因的数目, 是1050 + 869, 明显多于之前.

res <- decideTestsDGE(qlf.2vs1, p.value = 0.05, lfc = 0)
summary(res)

7 富集分析

https://www.jianshu.com/p/e133ab3169fa

https://zhuanlan.zhihu.com/p/347148653
https://cloud.tencent.com/developer/article/1838918

.libPaths("path_to/conda_py38/conda3/envs/r4.1/lib/R/library")
library(org.Hs.eg.db)
library(clusterProfiler)
library(enrichplot)

data<-read.table("control_treat.DESeq2.txt",header = T,sep="\t")
data$fold_change <- data$log2FoldChange
data <- na.omit(data)
gene_List<-data$fold_change
names(gene_List)<-data$X
gene_List<-sort(gene_List,decreasing = T)

gene_map<-bitr(names(gene_List), fromType = "SYMBOL", toType = "ENTREZID", org.Hs.eg.db)
gene_map<-dplyr::distinct(gene_map, SYMBOL, .keep_all=T)

gene_List.entrez<-gene_List
names(gene_List.entrez)<-gene_map[,2]
gene_List.entrez<-gene_List.entrez[!is.na(names(gene_List.entrez))]

kk2 <- gseKEGG(geneList     = gene_List.entrez,
               organism     = 'hsa',
               nPerm        = 1000,
               minGSSize    = 120,
               pvalueCutoff = 0.05,
               verbose      = FALSE)

write.table(kk2, 'GSEA.kegg.txt', col.names = TRUE, row.names=FALSE, sep = '\t', quote = FALSE)


#gseagene <- gseNCG(gene_List.entrez, nPerm=10000)
#p2 <- dotplot(gseagene, showCategory=30) + ggtitle("dotplot for GSEA")
#dotplot(kk2,title="Enrichment KEGG_dot")

png(file="gsea.png", bg="transparent")
dotplot(kk2,title="Enrichment KEGG_dot")
dev.off()

pdf(file="EnrichmentScore.pdf")#, bg="transparent", height=500,width=500)
gseaplot2(kk2,1:3,subplots = 1:2, base_size = 20,pvalue_table = F)
dev.off()

#png(file="barplot.png", bg="transparent")
#barplot(kk2, title="Enrichment KEGG_bar")
#dev.off()
#png(file="emapplot.png", bg="transparent")
#emapplot(kk2)
#dev.off()

#png(file="heatplot.png", bg="transparent")
#heatplot(kk2)
#dev.off()

风风是超人

关注

1
点赞
踩
6

收藏

觉得还不错? 一键收藏
打赏
0
评论
RNA检测流程

1 比对Hist2 下载https://github.com/DaehwanKimLab/hisat2./hisat2/hisat2-2.2.1/hisat2 \ -x /hg19index/ucsc.hg19 \ -1 fq1.gz \ -2 fq2.gz \ -S out.sam \ --dta --rna-strandness-x ：对hg19参考基因组做的索引cd hg19indexhisat2/hisat2-2.2.1/hisat2-build -p 4 ucsc.
复制链接

扫一扫