目的：
对illumina数据进行处理，利用 RNA-Seq 发现新的 RNA 变体和剪接位点，或量化 mRNA 以进行基因表达分析等。对两组或多组样本的转录组数据，通过差异表达分析和对所发现的差异表达基因集合进行功能富集分析以推断生物学功能。
数据准备：
数据下载：
- Human genome(GRCh38/hg3)：Index of /goldenPath/hg38/chromosomes (ucsc.edu)
- Mouse genome (GRCm39/mm39)：Index of /goldenPath/mm39/bigZips (ucsc.edu)
- D. melanogaster genome(BDGP Release 6 + ISO1 MT/dm6)：Index of /goldenPath/dm6/bigZips (ucsc.edu)
2）注释文件（.gtf/gff）下载
- Human genome Annotation: Index of /pub/release-104/gtf/homo_sapiens/ (ensembl.org)
数据质控：
可采用Trim Galore或fastqc进行质控：
Fastqc软件参数：
#示例：fastqc -o qc/ *.fastq.gz
# 基本格式
# fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
# 主要是包括前面的各种选项和最后面的可以加入N个文件
# -o --outdir FastQC生成的报告文件的储存路径，生成的报告的文件名是根据输入来定的
# --extract 生成的报告默认会打包成1个压缩文件，使用这个参数是让程序不打包
# -t --threads 选择程序运行的线程数，每个线程会占用250MB内存，越多越快咯
# -c --contaminants 污染物选项，输入的是一个文件，格式是Name [Tab] Sequence，里面是可能的污染序列，如果有这个选项，FastQC会在计算时候评估污染的情况，并在统计的时候进行分析，一般用不到
# -a --adapters 也是输入一个文件，文件的格式Name [Tab] Sequence，储存的是测序的adpater序列信息，如果不输入，目前版本的FastQC就按照通用引物来评估序列时候有adapter的残留
# -q --quiet 安静运行模式，一般不选这个选项的时候，程序会实时报告运行的状况。
Trim Galore软件参数：
trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 \
--paired $dir/cmp/01raw_data/$fq1 $dir/cmp/01raw_data/$fq2 \
--gzip -o $input_data
--hardtrim5：用于从序列的3’端切除碱基，示意如下
before: 5’-CCTAAGGAAACAAGTACACTCCACACATGCATA-3’
--hardtrim5 20: CCTAAGGAAACAAGTACACT
通过hardtrim5参数可以将序列截取成固定长度。
--hardtrim3参数：从序列的5’端切除碱基。
--quality：设定Phred quality score阈值，默认为20。
--phred33：选择-phred33或者-phred64，表示测序平台使用的Phred quality score。
--adapter：输入adapter序列。也可以不输入，Trim Galore!会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个，也直接显式输入这三种平台，即--illumina、--nextera和--small_rna。
--stringency：设定可以忍受的前后adapter重叠的碱基数，默认为1。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到。
--length：设定输出reads长度阈值，小于设定值会被抛弃。
--paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准。
--retain_unpaired：对于双端测序结果，一对reads中，如果一个read达到标准，但是对应的另一个要被抛弃，达到标准的read会被单独保存为一个文件。
--gzip和--dont_gzip：清洗后的数据zip打包或者不打包。
--output_dir：输入目录。需要提前建立目录，否则运行会报错。
-- trim-n : 移除read一端的reads。
一般步骤：
- 去除read 3‘端低质量碱基
- 去除adapter(接头)序列
有三种常见接头：
Illumina: AGATCGGAAGAGC
Small RNA: TGGAATTCTCGG
Nextera: CTGTCTCTTATA
- 去除很短的序列
- 其他
基本步骤：
Reads align
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1
chrX_data/samples/ERR188044_chrX_1.fastq.gz -2
chrX_data/samples/ERR188044_chrX_2.fastq.gz -S
ERR188044_chrX.sam
Sort sam to bam
samtools sort -@ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam
Assemble transcripts
$ stringtie -p 8 -G chrX_data/genes/chrX.gtf -o
ERR188044_chrX.gtf –l ERR188044 ERR188044_chrX.bam
Merge transcripts
stringtie --merge -p 8 -G chrX_data/genes/chrX.gtf -o
stringtie_merged.gtf chrX_data/mergelist.txt
Examine the transcripts compare to reference annotation
$ gffcompare –r chrX_data/genes/chrX.gtf –G –o merged
stringtie_merged.gtf
Estimate transcript abundances ,create table counts for Ballgown:
$ stringtie –e –B -p 8 -G stringtie_merged.gtf -o
ballgown/ERR188044/ERR188044_chrX.gtf ERR188044_chrX.bam
differential expression analysis
$ R
R version 3.2.2 (2015-08–14) -- "Fire Safety"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.4.0 (64-bit)
library(ballgown)
library(RSkittleBrewer)
library(genefilter)
library(dplyr)
library(devtools)
#Load the phenotype data for the samples
pheno_data = read.csv("geuvadis_phenodata.csv")
#Read in the expression data that was calculated by StringTie
bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR",
pData=pheno_data)
#Filter to remove low abundance genes
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >
1",genomesubset=TRUE)
#Identify transcripts that show statistically significant differences between
groups.
results_transcripts = stattest(bg_chrX_filt,
feature="transcript", covariate="sex", adjustvars =
c("population"), getFC=TRUE, meas="FPKM")
#Identify genes that show statistically significant differences between groups.
results_genes = stattest(bg_chrX_filt, feature="gene",
covariate="sex", adjustvars = c("population"), getFC=TRUE,
meas="FPKM")
#Add gene names and gene IDs to the results_transcripts data frame:
results_transcripts =
data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),
geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)
#Sort the results from smallest p-value to largest:
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)
#Write the results to a CSV (comma-separated values) file that can be shared
and distributed:
write.csv(results_transcripts, "chrX_transcript_results.csv",
row.names=FALSE)
write.csv(results_genes, "chrX_gene_results.csv",
row.names=FALSE)
#Identify transcripts and genes with a q-value of less than 0.05:
subset(results_transcripts,results_transcripts$qval<0.05)
subset(results_genes,results_genes$qval<0.05)
#Make the plots pretty.
#Show the distribution of gene abundances
And so on.
基本分析：
有参考基因组，需要组装新的转录本：
HISAT2 -> StringTie -> Ballgown(或)htseq-count + DESeq2
注：样本量如果大于10，推荐使用TACO，代替cuffmerge或stringtie
HISAT2的使用：
示例：
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran -1
chrX_data/samples/ERR188044_chrX_1.fastq.gz -2
chrX_data/samples/ERR188044_chrX_2.fastq.gz -S
ERR188044_chrX.sam
主要参数：
-x <hisat2-idx>
参考基因组索引文件的前缀（目录及文件前缀）。
-1 <m1>
双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
-2 <m2>
双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
-U <r>
单端数据文件。若有多组数据，使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。
–sra-acc <SRA accession number>
输入SRA登录号，比如SRR353653，SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型，进行比对。
-S <hit>
指定输出的SAM文件（目录及名称）。
-t/--time print wall-clock time taken by search phases
-p/--threads <int> number of alignment threads to launch 指定进程数目（8/16）
Samtools的使用
示例：
#version 1.3 or newer
samtools sort -@ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam
#version 1.2 or older
samtools view -bS ERR188044_chrX.sam > ERR188044_chrX_unsorted.sam
samtools sort -@ 8 ERR188044_chrX_unsorted.bam ERR188044_chrX
主要参数：
参考：Samtools使用大全
Stringtie的使用
主要参数：
|
参数 |
描述 |
|
-G <ref_ann.gff> |
使用注释好的gtf文件辅助组装，在-e未设置的条件下，输出中包括注释文件中的转录本和预测的新型转录本 |
|
-o [<path/>]<out.gtf> |
输出文件的名字，最好是全路径，默认输出为标准输出 |
|
-l <label> |
为输出的转录本设置前缀名，默认为STRG |
|
-p <int> |
线程数，默认为1 |
|
-A <gene_abund.tab> |
对输出的gtf统计基因表达量，并以一个tab分割的文件输出，这里需要提交输出的文件名 |
|
-C <cov_refs.gtf> |
对输出的gtf中属于-G提交的参考gtf的转录本统一输出到该文件，这里需要提交一个文件名 |
|
-B |
是否需要输出Ballgown可以识别的文件,在-b设置的情况下，使用-o的路径输出 |
|
-b <path> |
对Ballgown输出的文件指定一个路径保存 |
|
-e |
我认为是最需要注意的参数！！只统计可以匹配-G提交的参考gtf中的转录本，不再对新的转录本做预测，这可以加快程序的运行速度 |
|
-m <int> |
对预测的转录本设置最小长度，默认为200 |
Ballgown的使用
参考：Ballgown.pdf
有参考基因组，不需组装新的转录本：
参考：https://github.com/hnnd/NGS/blob/master/pip_b.md
无参考基因组：
trinity > rsem > DEseq2
参考：https://github.com/hnnd/NGS/blob/master/pip_c.md
其他分析：
- 可变剪切分析
- 融合基因分析
- GO富集分析
- Pathway富集分析
基本分析流程[1]：
Concept
Abundance estimate：表现为FPKM值，即每百万映射读取的转录片段每千碱基。
Total：clean data双端总reads数目Duplicate：重复的reads数目。
Mapped：比对到参考基因组上的总reads数目（比例）。
Properly mapped：比对到参考基因组且方向正确的reads数目（比例）。
PE mapped：双端reads比对到参考基因组上的reads数目（比例）。
SE mapped：仅单端read比对到参考基因组上的reads数目（比例）。
with mate mapped to a different chr：比对到不同染色体的reads数目。
with mate mapped to a different chr (mapQ>=5)：比对到不同染色体且比对质量不低于5的reads数目。
Average_sequencing_depth:比对到参考基因组的平均测序深度（测序数据量/基因组大小）。
Coverage：比对数据对全基因组区域的覆盖度（碱基覆盖长度占全基因组碱基总长的比例）。
Coverage_at_least_4X:全基因组区域中碱基覆盖深度不低于4X的比例。
Coverage_at_least_10X：全基因组区域中碱基覆盖深度不低于10X的比例。
Coverage_at_least_20X：全基因组区域中碱基覆盖深度不低于20X的比例。
RNA-seq可行的分析流程与方向
全转录组是指特定细胞在特定状态下所能转录出来的所有RNA的总和，包括mRNA和非编码RNA（non-coding RNA），针对非编码RNA的研究主要集中在具有调控作用的miRNA，lncRNA和circRNA。
基于高通量测序的RNA测序（RNA-Seq）是一种高度灵敏和准确的检测全转录组表达的工具，研究人员可以在单次分析中同时分析同一样本中的mRNA，miRNA，lncRNA，circRNA，检测已知的和发现新的特征，使转录本亚型、基因融合、单核苷酸位点变异和其他特征不受先验知识的限制。研究人员能够检测疾病状态、治疗反应、各种环境条件下以及一系列其他研究设计中以前无法检测的变化。全面快速地获得RNA序列和表达丰度信息，通过关联分析，还可深入挖掘生命现象背后的转录调控问题。
靶向RNA测序则是对感兴趣的转录本表达进行高度准确且特异测定的方法，同时提供定量和定性信息，能实现表达分析、融合基因鉴定。通过对特定基因的表达谱进行分析，可以评估疾病相关的变异和表观遗传学改变。基因融合和基因表达变化分析为癌症功能相关的改变提供了特定的视角。
Small RNA主要包括miRNA、piRNA、tsRNA（tRF&tiRNA）、snRNA和snoRNA等，是一类不具有蛋白编码能力的RNA分子，能调控基因表达，在细胞生长、发育和代谢等基础生物学过程中都扮演着重要的角色，甚至在癌症等相关疾病形成过程中也起着关键的作用。高通量测序技术省去了烦琐的Small RNA克隆文库构建过程，可以一次性生成上百万条Small RNA序列，无需预先假设，能够快速鉴定特定条件下表达的已知Small RNA并发现新的Small RNA，以单碱基分辨率鉴定变异，同时还可以研究不同条件下Small RNA的表达差异，并可以联合转录组测序表达谱数据进行关联分析，在转录和转录后水平研究基因调控。
基因结构水平分析：与参考基因组比对、新转录本预测、可变剪切分析、SNP/InDel注释分析、融合基因分析、癌基因注释等；基因表达水平分析：RNA-seq整体质量评估、差异基因分析、GO/KEGG富集分析、蛋白互作网络分析、转录因子注释等。
1、原始数据质控
以原始数据为研究对象，对于低质量序列，未检测序列，接头序列进行过滤，并对于过滤前后数据的碱基质量、GC含量、长度分布、接头留存和Duplication比率等指标进行分析。
2、RNA基因组比对（RNA Mapping）
采用Hisat2/Mapsplice/Star/Tophat2等算法进行基因组比对，得到基因组比对的bam文件，并基于bam文件进行信息统计，得到基因组比对率、reads在基因结构和染色体上的分布结果。
3、表达量统计（Expression）
采用HTSeq以及基因组注释的gff3文件，根据单端或双端测序类型，选择RPKM或FPKM的标化方式对基因表达量进行统计。基于统计结果，分析得到样本间相关性、 RPKM/FPKM密度和丰度等分析结果，反映单个样本基因表达水平分布和离散程度，以及不同样本整体基因表达水平的差异。
4、差异基因筛选（Dif Gene Analysis）
采用DESeq2/DESeq/EBSeq/EdgeR/Limma等算法进行差异筛选，得到满足差异倍数（Fold Change）以及FDR阈值的差异基因，并基于差异筛选结果以及样本的FPKM或RPKM，进行火山图分析（Volcano Plot）以及聚类图分析（Heatmap）。
5、功能分析（GO Analysis）
采用NCBI/UNIPROT/SWISSPROT/AMIGO等GO数据库，对于差异基因进行功能分析，从而得到差异基因所显著性富集的功能条目（GO Term）。
6、信号通路分析（Pathway Analysis）
整合KEGG、NCBI、EMBL等数据库，深入优化所需算法，对差异基因进行信号通路分析，得到差异基因所显著性富集的信号通路条目。
7、GO-Tree分析
采用GO数据库中GO-term的上下级层级从属关系，进行GO-Tree绘制，得到显著性差异功能的功能簇以及层级从属关系。
8、Path-Act-Network分析
采用KEGG数据库记载的信号通路上下游调控关系，进行Path-Act-Network绘制，得到宏观上的显著性信号通路的上下游调控关系。
9、共表达网络分析（Co-Exp-Network Analysis）高级分析
基于GO Analysis和Pathway Analysis得到的显著性条目，以及研究者感兴趣条目，以这些条目中基因的表达值为研究目标，进行共表达网络和K-Core分析，从而得到基因间的相关性和基因的核心度，再以Co-Expression.txt和K-Core为研究对象，采用Cytoscape进行图形化展示，得到Co-Expression-Network。
10、基因间相互作用关系网络(Gene-Act-Network Analysis) 高级分析
基于GO Analysis和Pathway Analysis得到的显著性条目，以研究者感兴趣的相关表型基因为研究对象，采用KEGG数据库基因间关系注释，帮助研究者绘制Gene-Act-Network，快速定位“核心”基因。
11、韦恩分析
以各分组间的基因为研究对象，采用韦恩作图分析的方法，可找出各分组间共有或者特有的差异表达基因并进行深入分析。
12、趋势分析
以各差异分组间的韦恩基因的FPKM值为研究对象，采用STEM算法，进行趋势分析，得到按照样本逻辑顺序所在趋势。
13、加权基因共表达网络分析（WGCNA）分析
基于加权的表达相关性，进行层级聚类分析，并根据设定标准切分聚类结果，获得不同的基因模块，采用聚类树的分枝和不同颜色来鉴定高度协同变化的基因集。如果有表型信息，还可以计算基因模块与表型相关性，鉴定性状相关的模块，并根据基因集的内连性和基因集与表型之间的关联鉴定候补生物标记基因或治疗靶点。
关于GO富集时的pathway选择
不要仅仅基于Pathway富集分析的结果解读数据，人为的解读和挑选是必不可少的。因为生物数据的解读，在现阶段更多是生物学问题，而不是数学问题。原因大体如下：
（1）基因调控是个系统，不要仅仅看成1个1个孤立的pathway。
在今年参加的第二届植物抗逆会议上，1个教授说了一句话，我认为很有道理。“在植物体内其实根本就不存在pathway，什么脱落酸通路，水杨酸通路，其实这些调控因子都是相互联通，相互影响的，是个整体。只是我们人类为了研究方便，人为将这些系统拆分各个子集。 ” 所以，如果你真的将pathway看成1个个破碎的途径，以为某种处理只会影响某个pathway，富集分析必须在数学上或统计学上得到1个指向性很强的结论，那是不大可能的。
具体说了，说基因调控是个系统，可以从两个层面进行解读：
a）1个基因的改变可以造成整个系统的改变；
举几个例子：
把1个生命活动必须的蛋白敲除后，整个细胞会发生紊乱。而植物抗病应激，也往往是1个受体蛋白识别了病原的外源蛋白，然后导致整个细胞系统的变化。
b）1个基因往往有多个功能，但执行具体的功能往往是不同蛋白复合物共同作用的结果。
例如。基因X理论上在不同情况下，有可能参与A、B、C通路。在某个生物处理下，或许基因X 只在A通路里起作用。但如果进行基因注释的话，X同样也会被注释到B、C。所以，富集分析的结果总是会涉及特别多的通路。例如，研究人的项目，无论什么研究背景，常常会富集到 帕金森综合症通路。不是你的材料真的得了帕金森综合症，只是那些与你实验处理相关的基因，在一定条件下也可以参与到帕金森综合症的过程，所以被注释到了这个通路里。
小结：所以，我们也看到了。无论什么实验处理，总有可能导致整个系统的变化。同时，基因的通路注释也有欺骗性。那么，从这一堆冗余信息中，想得到与我们研究相关的结论，离不开人为的筛选也解读。从那个复杂的整体中，筛选出核心的局部片段，这是个技术活。“这样的话是否存在一个问题就是在结果的解释上比较主观，也会因自身背景知识的不足而漏掉一些新颖的结果”。那当然，同样的结果1个外行可能什么都没有看见。但1个资深的学者可能会把握到很精彩的内容。好像任何领域都是如此，除了提高内功好像没有其他捷径。
（2）pathway富集分析的统计假设，并非在任何情况下都适用
pathway富集分析，在生物学上的假设是：1个pathway 上游基因的改变，会导致下游相关基因改变，从而改变通路中大量基因的表达，达到统计学上富集的效果。但很多pathway中，基因A、B、C并不是相互调控的关系，而是共同参与某个过程的不同部分。
例如，代谢物X的合成修饰。基因A、B、C分步骤参与合成的3个步骤。基因A 给X前体加了 羟基，然后传递到下游；基因B又给X前体加了苯环，再传递到下游；基因C又给X的前体 加了个乙酰基 ，完成X的合成。那么，基因A、B、C是参与了的相同的通路。如果基因A发生表达量变化，会直接调控影响B、C的表达量变化吗？ 看来很有可能不会，所以从RNA-seq差异分析的富集分析结果中，这个通路是不显著的。那么基因A的表达变化是否有生物学意义？ 当然有，因为代谢物X的合成的确受影响了。
类似的例子，理论上DNA差异甲基化的结果，就不能看pathway富集分析的结果。1个pathway 1个基因的DNA甲基化变化，就足以改变这个通路的基因表达，而不需要整个通路的甲基化都发生变化。DNA甲基化、组蛋白CHIP-seq的结果，其实只看功能注释、或通路注释就足够了，不需要考虑富集。
所以，我们还是要观察、理解某个核心pathway中基因的相互作用，才能判断其中的基因变化是否有生物学意义， 而不仅仅看富集分析的p value或Q value。
（3）目前的pathway是不完整的。
目前KEGG等数据库，收录的是已有的研究结果。但这些pathway的信息，远没有到达完善的水准。大部分通路只是了解1个大概的调控途径，而中间有什么转录因子参与、是否还有其他代谢物的生成，都是不知道的。这些通路的完整性，也会影响pathway富集分析结果。例如，基因A发生变化了，看起来下游基因没有变化。也许是还有其他的调控在起作用，只是这些调控作用现在还不知道而已。
总结：pathway 和 GO富集分析结果的解读，应该从生物学意义的角度出发，P value 和 Q value只是个参考而已，那些不显著的通路也值得解读（从功能注释的角度解读，而不是从富集分析的角度解读）。只要结果可以解释，有意义，不用太迷信P value。
If you want these code to practice it ，please contact me：zengcj2035@163.com