转录组分析学习笔记（持续补充）

时间 2019-11-24

标签转录分析学习笔记持续补充繁體版

原文原文链接

转录组分析流程(有参和无参de novo)

得到测序数据，Fastq格式，称之为Raw data。

质量检测

比对Mapping

Quantification|Quantitation

差别表达分析

补充：开始项目以前，先确立合理的文件目录结构。php

【1】Raw Data 处理

理论知识

高通量测序之因此可以可以达到如此高的通量的缘由就是他把原来几十M，几百M，甚至几个G的基因组经过物理或化学的方式打算成几百bp的短序列，而后同时测序。在测序过程当中，机器会对每次读取的结果赋予一个值，用于代表它有多大把握结果是对的。从理论上都是前面质量好，后面质量差。而且在某些GC比例高的区域，测序质量会大幅度下降。所以，咱们在正式的数据分析以前须要对分析结果进行质控。html

Fastq 文件

测序给的“原始数据”，称之为Raw Data。java

FASTQ是基于文本的，保存生物序列（一般是核酸序列）和其测序质量信息的标准格式。其序列以及质量信息都是使用一个ASCII字符标示，最初由Sanger开发，目的是将FASTA序列与质量数据放到一块儿，目前已经成为高通量测序结果的事实标准。
python

FASTQ文件中以四行最为一个基本单元，并对应一条序列的测序信息，各行记录信息以下：ios

第一行记录序列标识以及相关的描述信息，以‘@’开头，为了保证后续分析软件可以区分每条序列，单个序列的标识必须具备惟一性；
第二行为碱基序列；
第三行以‘+’开头，后面是序列标示符、描述信息，或者什么也不加；
第四行，是质量信息，长度和第二行的序列相对应，每个序列都有一个质量评分，根据评分体系的不一样，每一个字符的含义表示的数字也不相同。
碱基质量得分与错误率的换算关系： Q = -10log10p（p表示测序的错误率，Q表示碱基质量分数）
ASCII值与碱基质量得分之间的关系：
Phred64 Q=ASCII转换后的数值-64
Phred33 Q=ASCII转换后的数值-33
目前illumina使用的碱基质量格式为phred+33, 和Sanger的质量基本一致（老数据建议查看清楚再进行后续处理）。git

SRA 文件

Sequence Read Archive (SRA) makes biological sequence data available to the research community to enhance reproducibility and allow for new discoveries by comparing data sets. The SRA stores raw sequencing data and alignment information from high-throughput sequencing platforms, including Roche 454 GS System®, Illumina Genome Analyzer®, Applied Biosystems SOLiD System®, Helicos Heliscope®, Complete Genomics®, and Pacific Biosciences SMRT®.github

具体参考NCBI关于SRA的介绍数据库

SRA文件转换成fastq文件

咱们从NCBI等数据库下载的原始数据不少为SRA格式，须要转换成fastq文件。经常使用工具为：NCBI SRA Toolkit
简单介绍使用方法（具体的坑踩过才能记住！！！）express

单个文件转换(PE)sass

fastq-dump --gzip --split-3 -O outputdir -A file1.sra
多个文件批量转换

for I in seq 56 62
do
fastq-dump --gzip –split-3 -O ./fastq/ -A SRR35899${I}.sra
done
--split-3：若是是双端测序数据，则输出两个文件，若是不是则只输出一个文件
--gzip：输出格式为gzip的压缩文件（fastqc软件能够直接识别gzip压缩的文件）
-A：accession序列号，输入的文件
-O：outdir输出文件夹，指定输出路径

FastQC（测序质量分析）：多个文件批量进行

$ fastqc -q -t 4 -o ./fastqc_result/ .fastq.gz &
-t：调用核心数目
-q：安静运行，运行过程当中不会生成报告，在结束时将报告生成一个文件
-o ../path/to/file ：文件输出位置
. fq.gz：输入文件，当前目录下全部名字中有“ .fq.gz ”的文件

查看QC结果

一、单个查看：鼠标双击打开html文件查看

二、批量查看：使用 moltiqc软件: moltiqc *fastqc.zip

Fastqc结果报告关注重点：
1）.basic statistics
2）.per base sequence quality
3）.per base sequcence content
4）.adaptor content
5）.sequence duplication levels
最主要的几个指标是GC含量，Q20和Q30的比例以及是否存在接头（adaptor）、index以及其余物种序列的污染等。

Per base sequence quality，各位置碱基质量，每一个read各位置碱基的测序质量。横轴碱基的位置，纵轴是质量分数，Quality score=-10log10p（p表明错误率），因此当质量分数为40的时候，p就是0.0001，质量算高了。红色线表明中位数，蓝色表明平均数，黄色是25%-75%区间，触须是10%-90%区间。若任一位置的下四分位数低于10或者中位数低于25，出现“警告”；若任一位置的下四分位数低于5或者中位数低于20，出现“失败，Fail”。

Per base sequence content，碱基分布，对全部reads的每个位置，统计ATCG四种碱基的分布，横轴为位置，纵轴为碱基含量，正常状况下每一个位置每种碱基出现的几率是相近的，四条线应该平行且相近。当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，每每提示咱们有overrepresented sequence的污染。本结果前10个位置，每种碱基频率有明显的差异，说明有污染。当任一位置的A/T比例与G/C比例相差超过10%，报'WARN'；当任一位置的A/T比例与G/C比例相差超过20%，报'FAIL'。

Adapter content，接头含量：在此处可查看数据中使用接头信息
具体接头查询地点：github-fastqc 或者：Download common Illumina adapters

TruSeq Universal Adapter：
AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
Illumina Small RNA 3p Adapter：
1 ATCTCGTATGCCGTCTTCTGCTTG

另外，通常而言RNA-Seq数据在sequence deplication levels 未经过是比较正常的。毕竟一个基因会大量表达，会测到不少遍

Trimmomatic处理fastq数据

根据FastQC质控报告，利用Trimmomatic软件处理数据。trimmomatic，是java软件，因此直接Google找到其官网，而后下载二进制版本解压便可使用！这个软件设计就是为了illumina的测序数据的，由于它自带的adaptor文件有限！通常只去除TruSeq Universal Adapter 这个接头，运行的时候，不报错才算是成功的！
官网有例子，很简单的：http://www.usadellab.org/cms/?page=trimmomatic

Paired End:
java -jar trimmomatic-0.35.jar PE -phred33 input_forward.fq.gz input_reverse.fq.gz output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 ## 因此只须要把参数放对位置便可！

This will perform the following:

Remove adapters (ILLUMINACLIP:TruSeq3-PE.fa:2:30:10)
Remove leading low quality or N bases (below quality 3) (LEADING:3)
Remove trailing low quality or N bases (below quality 3) (TRAILING:3)
Scan the read with a 4-base wide sliding window, cutting when the average quality per base drops below 15 (SLIDINGWINDOW:4:15)
Drop reads below the 36 bases long (MINLEN:36)

！！！记住指定接头文件必定要用全路径哦！！！
固然，也能够用cutadapt这个python软件来去除接头的，可是它有一个弊端，须要本身指定接头文件。

【2】下载参考基因组及基因注释

在 UCSC 下载 hg19 参考基因组；

从 gencode 数据库下载基因注释文件，而且用 IGV 去查看感兴趣的基因的结构，好比TP53,KRAS,EGFR 等等。

截图几个基因的 IGV 可视化结构

下载 ENSEMBL，NCBI 的 gtf，也导入 IGV 看看，截图基因结构

了解 IGV 常识
来源于生信技能树：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
Y大宽连接：https://www.jianshu.com/p/02a92e4ead4b

【3】比对Mapping

把整理好的数据和参考基因组序列进行比对，寻找每一个reads的最佳匹配位置。可使用HISAT2，tophat2，STAR等软件。

具体参照Y叔介绍序列比对：Hisat2 pipeline（全流程，不只仅是Mapping）
Hisat2(比对,sam) ==> SAMtools (bam, sort, index) ==> htseq-count(reads 计数, 表达矩阵) ==> R(合并多个矩阵) ==> DEseq2(筛选差别表达基因并注释)

【4】表达量分析（软件：HTSeq，RSEM等）

RSEM (RNA-Seq by Expectation-Maximization)

RSEM^1,2 is an RNA-Seq transcript quantification program developed in 2009. You need a server with Linux/Mac OS. To run RSEM, your server should have C++, Perl and R installed. In addition, you need at least one aligner to align RNA-Seq reads for you. RSEM can call Bowtie, Bowtie 2 or STAR for you if you have them installed. Last but not least, you need to install the latest version of RSEM.
github连接

HTseq

htseq-count 是一款用于reads计数的软件，他能对位于基因组上的一些单位的reads数进行统计，这里所说的单位主要是指染色体上的一组位置区间（咱们常见的就是gene exon）。
HTSeq follows install conventions of many Python packages. In the best case, it should install from PyPI like this:

pip install HTSeq

If this does not work, please open an issue on Github .

基本用法：
htseq-count [options]

htseq-counts跟bedtools的区别

bedtools无论三七二十一，只要你的reads比对到基因组的坐标跟目的基因坐标有交叉，就算你一个reads，不须要管你是否是multiple mapping的。可是htseq就谨慎不少，并且还能够挑选model，通常来讲，它会把multiple mapping的reads归类到 not unique aligned里面。

最后htseq-counts使用的时候有一些参数尤为须要注意：

软件官网说明书： https://htseq.readthedocs.io/en/release_0.11.1/

参考gtf文件能够是gencode或者是ensembl数据库的，可是尤为要注释chr的问题，并且版本问题，gtf/gff格式无所谓。比对后的文件必定要进行sort，推荐必定要sort -n，根据reads的name来sort
-f sam/bam: 若是对bam文件进行counts，必须保证服务器的python安装了正确的pysam模块
-r name/pos: 通常状况下咱们的bam都是按照参考基因组的pos来sort的，可是这个软件默认倒是reads的name，很坑，通常建议从新把bam文件sort一下，而不是选择 -r pos，由于-r pos实在是太消耗内存了。
-s yes/no/reverse, 这也是巨坑的参数，默认是yes，通常人拿到的数据都是no，因此千万要注意！！！
-t 选择gff/gtf文件的第3列，通常是exon，也能够是gene，transcript ，这个不多调整的。
-i 这个须要修改，否则默认是ensembl的基因ID，通常人看不懂，能够改成gene_name，前提是你的gff文件里面有gene_name这个属性。

【5】差别表达分析（edgeR, DEseq2, EBseq等）

DEseq2

安装（坑本身踩！！！）

if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
BiocManager::install("DEseq2")

须要准备2个table：一个是countData，一个是colData

DESeq流程

> dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)
> dds <- DESeq(dds)
> res = results(dds, contrast=c("condition", "control", "treat"))
> 或下面命令
> res= results(dds)
> res = res[order(res$pvalue),]

> head(res)
> summary(res)
> 全部结果先进行输出
> write.csv(res,file="All_results.csv")
> table(res$padj<0.05)

提取差别基因

> diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
    或
> diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
> dim(diff_gene_deseq2)
> head(diff_gene_deseq2)
> write.csv(diff_gene_deseq2,file= "DEG_treat_vs_control.csv")

EBSeq

Bioconductor软件包安装：

if (!requireNamespace("BiocManager", quietly = TRUE))
             install.packages("BiocManager")
     BiocManager::install("EBSeq", version = "3.8")

TPR relatively independent of sample size and presence of outliers.

Poor FDR control in most situations, relatively unaffected by outliers.

Medium computational time requirement, increases slightly with sample size.

标准化方法不一样+检验不一样=多种组合/软件，用以前须要结合本身的样本量来考虑，多参考有类似实验设计的文献，经常使用的方法都跑一下，本身评估下结果差别，再作定夺。（研究原本就是充满了不肯定性，一切都只能用“可能性”来定义，因此，采用一样参数仍然没法彻底重复出文献中的结果也是常见。即便你心有芥蒂...）by fatlady：https://www.jianshu.com/p/364650e8bd3e

其余软件包

DEGseq(http://www.bioconductor.org/packages/release/bioc/html/DEGseq.html)
NOISeq(http://www.bioconductor.org/packages/release/bioc/html/NOISeq.html)
Ballgown (https://github.com/alyssafrazee/ballgown)是R语言中基因差别表达分析的工具，能利用RNA-Seq实验的数据(StringTie, RSEM, Cufflinks)的结果预测基因、转录本的差别表达。然而Ballgown并无不能很好地检测差别外显子，而 DEXseq、rMATS和MISO能够很好解决该问题。

【6】基因注释

正在纠结ID转换.....

【7】基因富集分析(gene set enrichment analysis, GSEA)

用clusterProfiler作富集分析

（1）Bioconductor安装软件包：http://www.bioconductor.org/install/

【7】数据可视化处理

根据不一样的PIPELINE选择合适的方式(R)进行可视化。
比较好的教程：

https://www.jianshu.com/p/72c871663e62

【8】转录组分析流程举例

有参转录组分析

1.首推Y叔分享的教程：RNA-seq分析简洁版

文章很精辟，在此不赘述。

2. Hisat2-Stringtie-Ballgown

案例来源：Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., & Salzberg, S. L. (2016). Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nature protocols, 11(9), 1650-67.doi: 10.1038/nprot.2016.095
数据和软件准备：
```
$ sudo yum install axel        
 $ axel ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.gz
 $ tar –zxvf chrX_data.tar.gz
 #下载和安装miniconda
 $ wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
 #下载完成后在终端中安装
 $ bash Miniconda-latest-Linux-x86_64.sh
 #按照提示安装，完成后
 $ source ~/.bashrc   #使以上的安装当即生效
 #输入如下命令检验miniconda是否安装成功
 $ conda list
```
利用conda install 软件名+版本号安装软件便可，咱们须要安装hisat二、stringtie、samtools三个软件，安装的命令为：
```
$ conda install hisat2
 $ conda install stringtie 
 $ conda install samtools
```

具体分析流程：
首先使用hisat2进行比对：

$ for i in {188044,188104,188234,188245,188257,188273,188337,188383,188401,188428,188454,204916};do hisat2 -p 4 -x chrX_data/indexes/chrX_tran -1 chrX_data/samples/ERR${i}_chrX_1.fastq.gz -2 chrX_data/samples/ERR${i}_chrX_2.fastq.gz -S ERR${i}_chrX.sam
    done

脚本执行完可获得12个sam文件。
经过samtools将sam文件转换为bam文件，做为stringtie的输入文件，具体脚本为：

$ for i in {188044,188104,188234,188245,188257,188273,188337,188383,188401,188428,188454,204916};do samtools sort -@ 4 -o ERR${i}_chrX.bam ERR${i}_chrX.sam done

此处sort默认输出的bam文件是按基因组位置排序，一样tophat的输出的bam文件也是按此顺序排序的，而sort -n 则是按reads的ID排序。

接下来利用stringtie对转录组进行组装，会针对每一个bam文件生成一个gtf文件：

$ for i in {188044,188104,188234,188245,188257,188273,188337,188383,188401,188428,188454,204916};do
    stringtie -p 8 -G ./genes/chrX.gtf -o ERR${i}_chrX.gtf -l ERR${i} ERR${i}_chrX.bam
    done

将输出的12个GTF文件的文件名录入到mergelist.txt文件中，而后利用软件stringtie将12个含有转录本信息的gtf文件合并成一个gtf。

$ stringtie --merge -p 8 -G ./genes/chrX.gtf -o stringtie_merged.gtf ./mergelist.txt

参数--merge 为转录本合并模式。在合并模式下，stringtie将全部样品的GTF/GFF文件列表做为输入，并将这些转录本合并/组装成非冗余的转录本集合。这种模式被用于新的差别分析流程中，用以生成一个跨多个RNA-Seq样品的全局的、统一的转录本。
接下来，从新组装转录本并估算基因表达丰度，并为ballgown建立读入文件。利用组装好的非冗余的转录本文件即stringtie_merged.gtf 和12个bam文件，执行下面的脚本:

$ for i in {188044,188104,188234,188245,188257,188273,188337,188383,188401,188428,188454,204916};do
    stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR${i}/ERR${i}_chrX.gtf ERR${i}_chrX.bam
  done

接下来要用到R语言分析（除dplyr，都为Bioconductor包）：

> library(RSkittleBrewer)
> library(ballgown)
> library(genefilter)
> library(dplyr)
> library(devtools)
#设置R语言的工做路径
> setwd("F:/data/R") 
#读取表型数据如右图所示
> read.csv("geuvadis_phenodata.csv")
> pheno_data<-read.csv("F:/data/R/geuvadis_phenodata.csv ")
#dataDir告知数据路径， samplePattern则依据样本的名字来，pheno_data    则指明了样本数据的关系，这个里面第一列样本名须要和ballgown下面的文件夹的样本名同样，否则会报错
> bg_chrX = ballgown(dataDir = “F:/data/R/ballgown", samplePattern = "ERR", pData=pheno_data) 
#滤掉低丰度的基因，这里选择过滤掉样本间差别少于一个转录本的数据
> bg_chrX_filt=subset(bg_chrX,"rowVars(texpr(bg_chrX))>1",genomesubset=TRUE)
#确认组间有差别的转录本，在这里咱们比较male和famle之间的基因差别，指定的分析参数为“transcripts”，主变量是“sex”，修正变量是“population”，getFC能够指定输出结果显示组间表达量的foldchange。
> result_transcripts=stattest(bg_chrX_filt,feature = "transcript",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
#查看有差别转录本的输出结果，以下图
 > result_transcripts

肯定各组间有差别的基因

>result_genes=stattest(bg_chrX_filt,feature = "gene",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
#为组间有差别的转录本添加基因名，以下图
> result_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt), geneIDs = ballgown::geneIDs(bg_chrX_filt), result_transcripts)
#按照p-value从小到大排序
> result_transcripts=arrange(result_transcripts,pval) 
> result_genes=arrange(result_genes,pval)
#把两个结果写入到csv文件中
> write.csv(result_transcripts, "chrX_transcript_results.csv",row.names=FALSE)
> write.csv(result_genes, "chrX_gene_results.csv",row.names=FALSE)
#从以上的输出中筛选出q值小于0.05的transcripts和genes，即性别间表达有差别的基因
> subset(result_transcripts,result_transcripts$qval<0.05)
> subset(result_genes,result_genes$qval<0.05)
#赋予调色板五个指定颜色
> tropical= c('darkorange', 'dodgerblue','hotpink', 'limegreen', 'yellow') 
> palette(tropical)
#以FPKM为参考值做图，以性别做为区分条件
> fpkm = texpr(bg_chrX,meas="FPKM")
#方便做图将其log转换，+1是为了不出现log2(0)的状况
> fpkm = log2(fpkm+1) 
> boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')

#查看单个转录本在样品中的分布，如图，并绘制箱线图
> ballgown::transcriptNames(bg_chrX)[12] 
> ballgown::geneNames(bg_chrX)[12]
>plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2), main=paste(ballgown::geneNames(bg_chrX)[12],' : ', ballgown::transcriptNames(bg_chrX)[12]),pch=19, xlab="Sex", ylab='log2(FPKM+1)')
>points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)), col=as.numeric(pheno_data$sex))

# plotTranscripts函数能够根据指定基因的id画出在特定区段的转录本，能够经过sample函数指定看在某个样本中的表达状况，这里选用id=1750, sample=ERR188234
> plotTranscripts(ballgown::geneIDs(bg_chrX)[1750], bg_chrX, main=c('Gene MSTRG.575 in sample ERR188234'), sample=c('ERR188234'))

3. 用tophat和cufflinks分析RNAseq数据

这个流程比较老，参考文章：

Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R., Pimentel, H., Salzberg, S. L., Rinn, J. L., … Pachter, L. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nature protocols, 7(3), 562-78. doi:10.1038/nprot.2012.016
https://www.nature.com/articles/nprot.2012.016

无参转录组分析

de novo组装

1. Trinity

Trinity是由 the Broad Institute 开发的转录组de novo组装软件，由三个独立的软件模块组成：Inchworm,Chrysalis和Butterfly。三个软件依次来处理大规模的RNA-seq的reads数据。Trinity的简要工做流程为：

Inchworm: 将RNA-seq的原始reads数据组装成Unique序列；
Chrysalis: 将上一步生成的contigs聚类，而后对每一个类构建Bruijn图；
Butterfly: 处理这些Bruijn图，依据图中reads和成对的reads来寻找路径，从而获得具备可变剪接的全长转录子，同时将旁系同源基因的转录子分开。

目前最经常使用 Illumina RNA-Seq data de novo组装软件。案例：

http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-554
https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-016-2633-2
http://www.nature.com/articles/srep08259
http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0128659

Trinity download

Build Trinity by typing 'make' in the base installation directory.
Assemble RNA-Seq data like so:

Trinity --seqType fq --left reads_1.fq --right reads_2.fq --CPU 6 --max_memory 20G

Find assembled transcripts as: 'trinity_out_dir/Trinity.fasta'

Trinity --seqType fq --max_memory 100G --CPU 50 --min_kmer_cov 3 --left   FCHK2FVCCXY_L3_WHDAVllgEAAARAAPEI-96_1.fq.gz,FCHK2FVCCXY_L3_WHDAVllgEAABRAAPEI-97_1.fq.gz,FCHK2FVCCXY_L3_WHDAVllgEAABRAAPEI-97_1.fq.gz  --right FCHK2FVCCXY_L3_WHDAVllgEAAARAAPEI-96_2.fq.gz,FCHK2FVCCXY_L3_WHDAVllgEAABRAAPEI-97_2.fq.gz,FCHK2FVCCXY_L3_WHDAVllgEAACRAAPEI-98_2.fq.gz --output gongtong_trinity_out  --group_pairs_distance 230 --no_version_check  --verbose --min_contig_length 250 --min_glue 3 --no_distributed_trinity_exec ~/bio/trinityrnaseq-Trinity-v2.4.0/trinity-plugins/parafly/bin/ParaFly -c recursive_trinity.cmds -CPU 50 -v

对于Trinity获得的转录本序列，Trinity官网推荐取每条基因中最长的转录本做为Unigene，并用这些Unigene去进行功能注释，可是在计算表达量的时候咱们依然会用到全部的转录本。
后续流程参考：无参转录组GO、KEGG富集分析——diamond+idmapping+GOstats

注意事项

Trinity分步运行
当数据量比较大的时候，trinity运行的时间会很长，同时，内存不够等状况出现的时候有可能程序运行崩溃。最好是分步运行。下一步会接着前一步进行下去。

Stage 1: generate the kmer-catalog and run Inchworm: –no_run_chrysalis

Stage 2: Chrysalis clustering of inchworm contigs and mapping reads: –no_run_quantifygraph

Stage 3: Chrysalis deBruijn graph construction: –no_run_butterfly

Stage 4: Run butterfly, generate final Trinity.fasta file. (exclude –no_ options)

计算资源
Ideally, you will have access to a large-memory server, ideally having ~1G of RAM per 1M reads to be assembled (but often, much less memory may be required).

The assembly from start to finish can take anywhere from ~1/2 hour to 1 hour per million reads (your mileage may vary). 我的记录了一次，使用dell服务器，64GB RAM，24 threads : 53M 的reads，运行了16.5h（平均3.2M/h），内存使用峰值为43G.

2. est2Assembly 和gsassembler 软件

案例： 2014 巴西橡胶树的研究，是一个综合多组织样本的RNA库，ployT建库，454测序，用的是est2Assembly 和gsassembler 软件作组装，用 NCBI RefSeq, Plant Protein Database 作注释，由于没有分组，因此没必要作差别分析，只须要找SNV和SSR标记便可，最后也是作GO/KEGG注释。

参考资料

https://www.cnblogs.com/chenpeng1024/p/9166988.html
HOPTOP转录组入门（三）：你懂质量控制吗？-转录组-生信技能树
转录组入门3-测序数据质量检查 | 分享自为知笔记
PANDA姐的转录组入门(3):了解fastq测序数据
(3)转录组之数据质控-转录组-生信技能树
转录组（3）：了解fastq测序数据 - 简书
RNA-seq分析简洁版
htseq-count的使用
RNA-seq(7): DEseq2筛选差别表达基因并注释(bioMart)
Count-Based Differential Expression Analysis of RNA-seq Data
http://guangchuangyu.github.io/2015/05/use-clusterprofiler-as-an-universal-enrichment-analysis-tool/
http://guangchuangyu.github.io/2016/01/go-analysis-using-clusterprofiler/