转录组分析是生信分析的一个基础和常见的分析流程，一般从下机数据开始，经过一系列的质控，基因组比对，差异分析等过程得到实验组与对照组（或者肿瘤中的转移组和原癌组）中差异表达的基因集，然后在基于该基因集结合我们的研究方向进行进一步的功能分析，信号通路分析，细胞间通讯，实验验证等得到一个可解释的实验结果。这篇文章主要就涉及到的基本流程和相应软件进行一个简单的介绍，后续更深入的知识点还需要小伙伴们自己去挖掘哦。
基本分析流程：
- 原始数据 -- Rawdata
- 质控 -- QC -- Trimmomatic，FastQC
- 基因组比对 -- Aligment -- STAR
- 比对结果量化 -- FeatureCount -- featurecount
- 两组间的基因表达差异分析 -- DEGs -- limma
- 差异基因的功能富集分析 -- 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
横轴为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比对到参考基因组上，一般分为有参考基因组比对和无参考基因组比对。我们常用的就是有参考基因组的比对。
- 如果数据库中存在该物种的物种基因组（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
比对率统计
总比对率和唯一比对率是我们关注的，尤其是总比对率。
- 总比对率低，一般可分为2种情况，总比对率低于10% ，可能是弄错了物种来源，比如鼠源细胞当成了人源的；总比对率低于60%，说明有外源基因组的污染，最常见的就是细胞培养的时候受到支原体污染。
比对区域分布统计
比对reads在每条染色体上的密度分布情况
可以通过IGV软件可视化查看BAM文件信息
具体以后再详细展开……
四. 比对结果的外显子（基因）表达计数 FeatureCount -- featurecount
在高通量测序分析中用于下游分析的关键信息是比对到每个genomic feature（外显子的表达量、featurecounts是一款使用于RNA-seq和DNA-seq的read summarization工具，应用了高效率的染色体哈希算法和feature区块技术因等）中的read数目。
featurecounts 速度快，需要的内存空间少，同时可以用于单端和双端的数据
通过STAR等软件比对得到的BAM文件，利用计数软件（如 featurecount, HTseq, stringtie），获取各基因的 reads count。
假设前提：认为比对到基因的reads 数与该基因的表达量成正相关。
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。
五. 不同分组间的基因差异表达分析 DEGs -- limma
根据样本的实验设计，利用limma包对不同样本组之间筛选差异表达的已知转录本。满足pvalue <= 0.05和大于等于2倍差异表达范围筛选两组之间的差异转录本
差异基因可视化 -- 火山图
针对组间筛选的不同差异基因集， 采用venn方法进行共同和特有的比较分析
针对样本组间筛选的差异表达基因， 采用对基因和样本进行双向层级聚类并且用热图显示
六. 差异基因的功能富集分析 -- 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中的差异基因都显著富集到了哪些通路上面。
还可以更详细地查看通路富集情况
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")
富集通路的条形图，可以显示差异基因富集到的通路有哪些。
