RNAseq常规流程

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

基本分析流程:

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

一. 测序数据下机后处理

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

  • 一般研究人员会把自己的样本(实验组/对照组)送到测序公司进行测序,然后测序公司会返回给我们一个fq.gz文件,如下. 每个样本都有两条链的测序数据 *_1.fq.gz 和 *_2.fq.gz , 我们要把他们单独放在一个文件夹下作为一个样本进行后续分析。


    fq.gz.png

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

  • 因为公司给的数据里面,是没有去掉接头的,所以需要跑之前先把接头去掉,过滤质量后再用。
    使用Trimmomatic(或者Trim Galore)切除数据中的接头序列和低质量序列,专门清洗illumina测序数据的常用工具
    优点:
    1、可直接对.fq.gz压缩文件操作
    2、结果默认生成在同文件的目录下,不生成文件夹,直接生成四个*.fq.gz文件。若需要有文件夹,需要额外建立文件夹的指令。参数-o,为文件夹(文件)
    trim_result.png

    Trimmomatic的执行文件是一个Java文件
### 切除接头序列
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标准错误

使用参数

  • PE/SE
    rawdata为 Paired-End,还是 Single-End的reads,其输入和输出参数稍有不同。
  • -threads
    设置多线程运行数
  • -phred33
    设置碱基的质量格式,可选pred64
  • ILLUMINACLIP:TruSeq3-PE.fa:2:30:10
    切除adapter序列。参数后面分别接adapter序列的fasta文件:允许的最大mismatch数:palindrome模式下匹配碱基数阈值:simple模式下的匹配碱基数阈值。
  • LEADING:3
    切除首端碱基质量小于3的碱基
  • TRAILING:3
    切除尾端碱基质量小于3的碱基
  • SLIDINGWINDOW:4:15
    Perform a sliding window trimming。Windows的size是4个碱基,其平均碱基质量小于15,则切除。
  • MINLEN:50
    最小的reads长度
  • CROP:
    保留reads到指定的长度
  • HEADCROP:
    在reads的首端切除指定的长度
  • TOPHRED33
    将碱基质量转换为pred33格式
  • TOPHRED64
    将碱基质量转换为pred64格式

2.FastQC (测序数据质控)

  • FastQC 是高通量测序数据的高级质控工具,输入FastQ,SAM,BAM文件,输出对测序数据评估的网页报告。
    一般在Trimmomatic之后,用FastQC软件对预处理数据进行质量控制分析。
    优点:
    1、可直接对.fq.gz压缩文件操作
    2、结果默认路径生成在同文件的目录下,生成的文件夹为XX_fastqc,里面有fastqc_data.txt , fastqc_report.html等重要文件
    fastqc $trim_read1 $trim_read2
    fqc.trim.png

    fqc_detail.trim.png

    3、可生成在指定的文件夹(应该说是路径吧),需要额外建立文件夹,用-o参数,为文件夹(文件)
    fastqc -t 12 -o out_path sample1_1.fq sample1_2.fq
    使用说明
  • -t --threads:线程数
  • -o --outdir:输出路径
  • --extract:结果文件解压缩
  • --noextract:结果文件压缩
  • -f --format:输入文件格式.支持bam,sam,fastq文件格式

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

image.png

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

  • 每个碱基位置的质量分数分布图


    image.png

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

  • 平均质量分数分布图


    image.png

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

  • 每个碱基位置的碱基N百分比


    image.png

    N占比越低,数据质量越高

  • 序列的核酸组成偏好性


    image.png

    正常情况下,我们期望的是A=T=C=G,GC含量较高时,经常会出现A=T,C=G的情况。

  • 每条reads的GC占比


    image.png

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

  • 重复reads统计


    image.png

    横轴表示重复的次数,纵轴表示重复的reads的数目

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

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

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

1. 比对原理

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


aligment.png
  • 如果数据库中存在该物种的物种基因组(fasta)及其基因组注释文件(gff, gtf),那么,我们就可以用该基因组信息作为参考基因组进行clean reads序列比对,便可以知道基因组的哪些位置转录了RNA序列,通过对reads进行计数,而得到基因表达的定量数据。
  • 如果没有可用的参考基因组及其注释信息,则需要考虑无参考基因组的比对方法进行比对。首先先将reads 进行组装拼接,得到从头组装的组装转录本,然后再将clean reads比对到组装转录本上得到比对结果。

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

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 \  # 输出文件的位置和前缀

结果解读

  • 比对结果:Aligned.out.sam或者Aligned.out.bam
  • 比对完以后的统计信息: Log.final.out
  • 剪切的信息:SJ.out.tab

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

  • header,以@开头,包含:
    @HD 版本信息;
    @SQ 比对参考序列信息;
    @RG 测序平台信息;
    @PG 序列比对使用的软件信息
  • record为基因组比对信息,每一行信息为:
    reads ID
    比对到基因组上的信息
    比对到的染色体
    比对上的位置信息
    比对质量值
    碱基组成
    碱基质量值等


    bam.png

比对率统计


mapcount.png

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

  • 总比对率低,一般可分为2种情况,总比对率低于10% ,可能是弄错了物种来源,比如鼠源细胞当成了人源的;总比对率低于60%,说明有外源基因组的污染,最常见的就是细胞培养的时候受到支原体污染

比对区域分布统计


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就是基于不同的分类思想而储存的基因相关功能的数据库。

  • GO和KEGG是两个数据库,里面有每个基因相关的功能信息,而富集分析就是一个把这些功能进行进行整合计算的算法。
  • GO和KEGG富集的R包clusterProfiler

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
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,817评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,329评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,354评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,498评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,600评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,829评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,979评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,722评论 0 266
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,189评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,519评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,654评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,329评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,940评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,762评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,993评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,382评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,543评论 2 349

推荐阅读更多精彩内容