轉錄組分析學習筆記（持續補充）

時間 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註釋。