RNA-seq【分析记录】

本文详细介绍了RNA-seq数据分析的完整流程,包括原始数据质量控制、转录本比对、表达定量以及差异表达分析。重点讨论了Trimmomatic的过滤策略、STAR和HiSAT2的比对工具、RSEM和HTSeq的表达量计算,以及DESeq2的差异基因分析。此外,还提到了链特异性文库的构建和无参转录本组装工具Trinity的应用。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

视频学习笔记

在这里插入图片描述

Isoforms at the RNA level are readily characterized by cDNA transcript studies.(secondary structures in cellular RNAs)

Alternative RNA splicing
In different human tissues, there are four distinctive isoforms of TH for both TH mRNA and protein. TH mRNA isoforms are produced by differential splicing of a single-copy gene using two splicing donor sites in the first exon and inclusion/exclusion of the second exon. The alternative RNA splicing produces mRNA variants from which proteins are produced that differ in amino acid length, potential phosphorylation, and subtle kinetic properties.


链特异性文库

由于转录组文库构建过程中要将mRNA反转录为双链的cDNA,因此,在测序时,即会同时得到cDNA双链的序列信息,这就消除了mRNA单链的特性,无法区分cDNA中的正义链和反义链。

因此提出了另一种文库构建方式,链特异性文库,其在文库制备和测序中保留mRNA 的方向信息,使得有效数据增多,基因表达定量更准确;基因注释更准确,新基因识别等分析更可靠;同该方法能够鉴定与开放阅读框反义的ncRNA、发现反义转录本。

其建库原理与普通转录组类似,不同之处在于合成第二链cDNA时,用dUTP代替dTTP,此时第二链cDNA上布满了含dUTP的位点。

在接头连接后用一种特定的酶,其可在尿嘧啶位置产生一个单核苷酸缺口,从而消化掉第二链cDNA,只保留第一链cDNA,随后进行PCR扩增纯化,得到只含第一链cDNA信息的文库。


概述

  • 预处理
    • FASTQC数据质控
    • 数据过滤Trimmomatic
  • 有参转录组序列比对
    • 不同比对软件的比较
    • 常用序列比对算法
    • 基因组比对
      • STAR
      • hisat2
    • 转录本比对 RSEM
  • 无参转录组
    • 转录本从头拼接原理
    • 拼接方法 Trinity
  • 表达定量
    • RNA-seq常用统计定量单位
    • 基因组比对
      • Htseq-Count
      • FeatureCount
    • 转录本比对 RSEM
  • 数据分析
    • 思路?
    • 差异表达
      + Deseq标准化原理
    • 富集分析
      + GO
      + 通路富集分析
    • 数据可视化展示
      • IGV,基因浏览器

前期准备

  • 构建项目目录
	mkdir 00ref 参考序列fasta,注释文件gtf/gff
	mkdir 01raw_data 
	mkdir 02clean_data
	mkdir 03align_out
	mkdir 04quantification
	mkdir 05Deqseq ...
  • 参考序列/注释信息下载
    • Ensembl数据库,JGI数据库,模式植物拟南芥(tair,ARAPORT)
  • 原始数据上传,检查完整性md5值
    • 给自己的文件生成md5值:md5sum *gz > md5.txt
    • 比对已有的md5值(check):md5sum -c md5.txt

数据预处理

质量控制
  • FASTQC
fastqc sample*.gz
for i in `ls *gz` ;do fastqc $i ;done
ls *gz |xargs -I [] echo 'nohup fastqc [] -o ./ &' > fastqc.sh
bash fastqc.sh
multiqc ./
  • MultiQC
质量过滤
  • Trimmomatic特点:多线程,快。主要用于去除illumina平台接头。根据碱基质量对fastq进行筛选,支持SE/PE,支持gzip/bzip2。

过滤依据:
ILLUMINACLIP:过滤reads中illumina接头
LEADING:reads开头切除质量值低于阈值的碱基
TRAILING:reads末尾切除质量低于阈值碱基
SLIDINGWINDOW:从reads的5‘端滑窗过滤,切除碱基质量平均值低于阈值的滑窗
MINLEN:丢弃剪切后长度低于阈值reads
TOPHRED33:碱基质量体系phred33
TOPHRED64:碱基质量体系phred64

  • 接头序列的选择
    • “Illumina Single End” / “Illumina Paired End” (看fastqc信息): “TruSeq2-SE.fa” and “TruSeq2-PE.fa”
    • “TruSeq Universal Adapter” / “TruSeq Adapter, Index …” (看fastqc信息): “TruSeq3-SE.fa” and “TruSeq3-PE.fa”
    • 去接头参数的选择:true(选择),false
trimmomatic PE -threads 4 \
s1_R1.fastq.gz s1_R2.fastq.gz \
cleandata/s1_paired_clean_R1.fastq.gz \
cleandata/s1_unpair_clean_R1.fastq.gz \
cleandata/s1_paired_clean_R2.fastq.gz \
cleandata/s1_unpair_clean_R2.fastq.gz \
ILLUMINACLIP:path/miniconda3/share/trimmomatic-0.36-5/adapters/TruSeq3-PE.fa:2:30:10:1:true \
LEADING:3 TRAILING:3 \
SLIDINGWINDOW:4:20 MINLEN:50 TOPHERD33

序列比对

无参转录组
  • 使用拼接工具组装转录本trinity
基于基因组比对(以染色体为单位)
  • STAR
    • 第一个通过算法优化将比对时间大幅降低的比对工具
    • 提供完善的输出内容
    • 消耗相对大的内存
  • Hisat2
    • tophat继任者
    • STAR启发下的后起之秀,耗时少内存低
    • 缺点:输出结果仅为比对文件,结果单一
基于转录本比对(以转录本为单位)
  • RSEM
    • 需要提前借助基因组和注释信息准备相关文件
    • 结合bowtie2和STAR进行比对和定量分析
STAR实例
  • 建立索引
#拟南芥
#build index
STAR --runThreadN 6 --runMode genomeGenerate \#索引模式mode
--genomeDir arab_STAR_genome \
--genomeFastaFiles 00ref/TAIR10.Chr.all.fasta \
--sjdbGTFfile 00ref/Araport11_GFF3_genes_transposons.201606.gtf \
--sjdbOverhang 149#reads长度-1:150-1,用于可变剪切预测?默认100

请添加图片描述

  • 比对
#STAR align 简洁版
STAR --runThreadN 5 --genomeDir arab_STAR_genome \
--readFilesCommand zcat \#文件是gz需要解压缩
--readFilesIn 02clean_data/s1_paired_clean_R1.fastq.gz \
02clean_data/s1_paired_clean_R2.fastq.gz \
--outFileNamePrefix 03align_out/s1_

请添加图片描述

#STAR align 复杂版
STAR --runThreadN 5 --genomeDir arab_STAR_genome \
--readFilesCommand zcat \#文件是gz需要解压缩
--readFilesIn 02clean_data/s2_paired_clean_R1.fastq.gz \
02clean_data/s2_paired_clean_R2.fastq.gz \
--outFileNamePrefix 03align_out/s2_ \
--outSAMtype BAM SortedByCoordinate \
--outBAMsortingThreadN 5 \
--quantMode TranscriptomeSAM GeneCounts

#--quantMode TranscriptomeSAM will output alignments translated into transcript coordinates, 为了使用RSEM进行定量分析做准备。

–quantMode TranscriptomeSAM
GeneCounts:这个参数需要特别注意,如果你后续要使用RSEM进行定量,就需要添加TranscriptomeSAM,不然生成的BAM文件RSEM无法使用的,GeneCounts生成用于后续featureCounts分析的输入文件

"FATAL ERROR: could not create output file ... Solution: check that you have permission to write this file"
#先查看你的系统打开文件限制数目
$ ulimit -n
1024
#设置一个更大的数目
$ ulimit -n 10000

请添加图片描述

  • 查看比对文件

表达定量

处理原始比对文件
  • picard/samtools
  • sam转bam
  • bam排序
  • 去除比对得分较低序列
  • 如果需要去除重复reads
三个流程
  • STAR+RSEM(先比对再定量)
    • 输出结果可以选择转录本定量或基因定量
    • 定量单位包括feature count,FPKM,TPM
    • 操作相对复杂
  • STAR+HTseq(先比对再定量)
    • 输出结果为原始read count
    • 结果可用于差异表达分析
    • 操作相对简单
  • Kallisto (free-alignment)
    • 速度快省内存
    • 基于转录本定量
    • 不产生bam不方便后续分析
STAR+RSEM实例
  • 准备定量分析所需文件
#rsem prepare reference构建索引文件
rsem-prepare-reference --gtf 00ref/***.gtf \
00ref/***.fasta \
arab_RSEM/arab_rsem

请添加图片描述

  • 利用STAR结果进行定量分析
rsem-calculate-expression --paired-end --no-bam-output \
--alignment -p 5 \
-q 03align_out/s2_Aligned.toTranscriptome.out.bam \
arab_RSEM/arab_rsem \
04rsem_out/s2_rsem

请添加图片描述
基于基因定量结果,基于转录本定量结果,统计信息目录

Kallisto实例
  • 利用转录本参考序列文件构建索引
kallisto index -I arab_kallisto ../arab_RSEM/arab_rsem.transcripts.fa

请添加图片描述

  • 进行无比对定量分析
kallisto quant -I arab_kallisto/arab_kallisto -o 05kallisto_out/s2 \
02clean_data/s2_paired_clean_R1.fastq.gz \
02clean_data/s2_paired_clean_R2.fastq.gz 

请添加图片描述
请添加图片描述

Htseq-count
htseq-count -r pos -m union -f bam -s no \
-q 03align_out/s2Aligned.sortedByCoord.out.bam >
05htseq_out/s2.htseq.out
STAR + featureCounts = STAR + HTSeq 升级版
  • 安装: conda install subread 优点:速度很快
featureCounts -p -a path/***.gtf -o out.counts.txt -T 6 -t exon -g gene_id *.Aligned.sortedByCoord.out.bam

差异分析

表达定量结果转表达矩阵
  • RSEM自带脚本
#构建表达矩阵
rsem-generate-data-matrix *_rsem.genes.results > output.matrix
  • 去除所有样本表达量为0的基因
#删除未检测到gene
awk 'BEGIN{printf "geneid\ta1\ta2\tb1\tb2\n"}{if($2+$3+$4+$5>0)print $0}' output.matrix > deseq2_input.txt
edgeR
DESeq2
mkdir 06deseq_out

deseq.r
请添加图片描述

###查看文件的行数
wc -l output.matrix 
wc -l deseq2_input.txt 

mv deseq2_input.txt ../06deseq_out/
mv deseq2.r ../06deseq_out/

abundance_estimates_to_matrix.pl
run_DE_analysis.pl
R脚本

### 设置当前目录;
打开一个文件夹中的R脚本,session
或 输入路径;setwd("~/00")



# read.table ,文件有header,第一行为row.name;
input_data <- read.table("deseq2_input.txt", header=TRUE, row.names=1)

# 取整,函数round
input_data <-round(input_data,digits = 0)

# 准备文件
# as.matrix 将输入文件转换为表达矩阵;
input_data <- as.matrix(input_data)

# 控制条件:因子(对照2,实验2)
condition <- factor(c(rep("ct1", 2), rep("exp", 2)))

# input_data根据控制条件构建data.frame
coldata <- data.frame(row.names=colnames(input_data), condition)

library(DESeq2)
# build deseq input matrix 构建输入矩阵
#countData作为矩阵的input_data;colData Data.frame格式;控制条件design;
dds <- DESeqDataSetFromMatrix(countData=input_data,colData=coldata, design=~condition)

# DESeq2 进行差异分析
dds <- DESeq(dds)

#实际包含3步,请自行学习;
# 提取结果
#dds dataset格式转换为DESeq2 中result数据格式,矫正值默认0.1
res <- results(dds,alpha=0.05)
#查看res(DESeqresults格式),可以看到上下调基因
summary(res)

#res(resultset)按照校正后的P值排序
res <- res[order(res$padj), ]
res

# 将进过矫正后的表达量结果加进去;
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)
names(resdata)[1] <- "Gene"

#查看(resdata)
# output result 输出结果

write.table(resdata,file="diffexpr-results.txt",sep = "\t",quote = F, row.names = F)

#可视化展示
plotMA(res)
maplot <- function (res, thresh=0.05, labelsig=TRUE,...){
with(res,plot(baseMean, log2FoldChange, pch=20, cex=.5, log="x", ...))
with(subset(res, padj<thresh), points(baseMean, log2FoldChange, col="red", pch=20,cex=1.5))

}
png("diffexpr-malot.png",1500,1000,pointsize=20)
maplot(resdata, main="MA Plot")
dev.off()

#画火山图
install.packages("ggrepel")
library(ggplot2)
library(ggrepel)

resdata$significant <- as.factor(resdata$padj<0.05 & abs(resdata$log2FoldChange) > 1)
ggplot(data=resdata, aes(x=log2FoldChange, y =-log10(padj),color =significant)) +
geom_point() +
ylim(0, 8)+
scale_color_manual(values =c("black","red"))+
geom_hline(yintercept = log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
geom_vline(xintercept = c(1,-1),lty=4,lwd=0.6,alpha=0.8)+
theme_bw()+
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"))+
labs(title="Volcanoplot_biotrainee (by sunyi)", x="log2 (fold change)",y="-log10 (padj)")+

theme(plot.title = element_text(hjust = 0.5))+
geom_text_repel(data=subset(resdata, -log10(padj) > 6), aes(label=Gene),col="black",alpha =0.8)


差异基因的提取
#查看符合条件的差异基因
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |head
#查看差异基因有多少行
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |wc -l
#提取某几列
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 |head
#上调基因的提取;
awk '{if($3>1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > up.gene.txt
#下调基因的提取;
awk '{if($3<-1 && $7<0.05)print $0}' diffexpr-results.txt |cut -f 1,2,4,7 > down.gene.txt

根据第1列是Geneid,第7,8,9,10列是counts数,用awk提取出geneID和counts。
awk -F '\t' '{print $1,$7,$8,$9,$10}' OFS='\t' out_counts.txt >subread_matrix.out

<think>好的,我现在需要帮助用户了解RNA-seq数据分析的完整流程,包括步骤和工具推荐。首先,我得回顾自己掌握的知识,并确保符合用户提到的系统级指令,比如正确使用LaTeX格式、生成相关问题等。 用户的问题是关于RNA-seq分析的详细流程,步骤和工具推荐。我需要先理清整个流程的大致步骤,从数据准备到下游分析,每一步需要什么工具,以及注意事项。同时,用户提供的三个引用资料可能涉及具体步骤和工具,比如引用2提到的STAR和FeatureCounts,Fastqc,引用3中的文件夹结构和一些分析步骤。需要结合这些信息来组织回答。 首先,我应该按照典型的RNA-seq流程来分步骤说明。通常流程包括:数据准备、质量控制、数据过滤、比对、定量表达、差异分析、功能注释等。每个步骤需要对应的工具,比如FastQC用于质量控制,Trim Galore或Trimmomatic进行过滤,STAR或Hisat2进行比对,FeatureCounts或HTSeq-count进行定量,DESeq2或edgeR用于差异分析,然后进行GO/KEGG富集分析等。 需要检查用户提供的引用内容,确保提到的工具和步骤被涵盖。例如,引用2提到STAR+FeatureCounts的组合,引用3中提到了Salmon和Hisat2,还有文件夹结构的准备。这些都需要整合到回答中。 同时,要注意回答的结构,分步骤说明,每个步骤的要点和工具推荐。可能还需要注意事项或常见问题,比如质量控制的重要性,比对工具的选择依据,差异分析工具的比较等。 然后,生成相关问题部分需要考虑用户可能进一步想知道的内容,比如各步骤的详细参数设置、工具之间的比较、数据可视化方法、处理大数据时的优化策略等。需要确保问题覆盖用户可能感兴趣的扩展方向。 另外,要确保符合系统级指令中的格式要求,比如行内数学公式用$...$,独立公式用$$...$$,但RNA-seq流程可能较少涉及数学公式,不过如果有统计方法如差异表达中的p值计算,可能需要提到,比如DESeq2使用的负二项分布模型,此时可能需要使用LaTeX公式,例如:差异表达分析可能涉及计算$p$值,并通过多重检验校正得到$q$值[^1]。 还要检查是否有需要单独成段的公式,如果有的话使用$$...$$格式。不过RNA-seq流程可能更多是步骤和工具,数学公式可能不多,但如果有必要的话要正确应用。 最后,确保回答使用中文,结构清晰,引用标识自然添加在引用的内容末尾。例如,当提到数据质控步骤时,可以引用引用2中的Fastqc步骤,或者在提到文件夹结构时引用引用3中的内容。 总结起来,回答的结构应该是: 1. 准备工作(文件夹结构,引用3) 2. 质量控制(Fastqc,引用2) 3. 数据过滤(Trim Galore/fastp,引用3) 4. 比对(STAR/Hisat2/Subjunc,引用2和3) 5. 定量(FeatureCounts/Salmon,引用2和3) 6. 差异分析(DESeq2/edgeR/limma) 7. 功能分析(clusterProfiler) 8. 注意事项(数据备份、版本控制、参数优化) 9. 工具推荐表格 10. 相关问题 每个步骤中需要提到工具、作用、可能的注意事项,并适当添加引用。例如,在数据准备阶段,引用3提到建立文件夹结构的重要性,可以引用;在定量阶段,引用2提到FeatureCounts,可以加上引用标识。 现在需要组织这些内容,确保流程逻辑正确,工具推荐合理,并自然插入引用。同时,生成的相关问题要覆盖流程的优化、工具比较、扩展分析等方向。</think>### RNA-seq数据分析完整流程及工具推荐 RNA-seq分析流程可分为**上游分析**(数据质控与预处理)和**下游分析**(表达量与生物学解释)两大部分,以下为详细步骤及工具推荐: --- #### 一、上游分析:数据准备与预处理 1. **准备工作(文件组织)** 建立标准化文件夹结构,便于管理原始数据、中间文件和结果,例如: ```bash mkdir -p data/rawdata data/cleandata/trim_galore mkdir -p Expression/featureCounts Mapping/Hisat2 ``` 强调版本控制和备份原始数据[^3]。 2. **质量控制(Quality Control)** - **工具**:`FastQC`(可视化质控报告)、`MultiQC`(整合多样本报告) - **目的**:检查测序数据质量(如碱基质量分布、GC含量、重复序列比例)[^2]。 3. **数据过滤与修剪** - **工具**:`Trim Galore`(基于Cutadapt的自动修剪)、`fastp`(高通量过滤) - **操作**:去除低质量碱基(如Phred评分<$20$)、接头序列、N碱基过多片段。 4. **序列比对(Alignment)** - **工具选择**: - 基因组比对:`STAR`(快速长读长)、`Hisat2`(内存高效) - 转录组比对:`Subjunc`(支持可变剪切) - **输出**:BAM/SAM文件,需用`samtools`排序和去重。 --- #### 二、下游分析:定量与生物学解释 5. **表达定量(Quantification)** - **基于比对的方法**: - `FeatureCounts`(计数基因区域 reads) - `HTSeq-count`(灵活处理重叠区域) - **无比对的方法**:`Salmon`(快速转录本丰度估计)。 6. **差异表达分析** - **工具**: - `DESeq2`(基于负二项分布,适用于复杂设计) - `edgeR`(高效处理小样本) - `limma-voom`(适用于log转化后的数据) - **关键输出**:差异基因列表($|\log_2(\text{FC})| >1$且$q\text{-value}<0.05$)[^1]。 7. **功能富集分析** - **工具**:`clusterProfiler`(GO/KEGG富集)、`GSEA`(基因集富集) - **可视化**:`ggplot2`绘制富集条形图或气泡图。 --- #### 三、注意事项 - **数据备份**:原始数据不可修改,中间文件建议压缩存储。 - **版本控制**:记录软件版本(如`STAR --version`)和分析参数。 - **计算资源**:比对步骤需较高内存(如STAR建议≥32GB)。 --- #### 工具推荐对比表 | 步骤 | 推荐工具 | 特点 | |---------------------|-----------------------------------|--------------------------------| | 质量控制 | FastQC + MultiQC | 可视化直观,支持多样本整合 | | 序列比对 | STAR(基因组) / Salmon(无比对)| 平衡速度与准确性 | | 差异分析 | DESeq2 | 统计严谨,支持复杂实验设计 | | 功能富集 | clusterProfiler | 支持多种数据库,集成可视化 | ---
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值