前提已經(jīng)獲得了trinity生成的.fasta文件(trinity_out_dir_all.Trinity.fasta)
在此，我們采用的是Trinity自帶的腳本進(jìn)行定量計(jì)算。
一：首先生成如下的一個(gè)**samples.txt**文檔，包含了所有的樣本：
cond_J??? cond_J1??? J1_1.fq.gz??? J1_2.fq.gz
cond_J??? cond_J2??? J2_1.fq.gz??? J2_2.fq.gz
cond_J??? cond_J3??? J3_1.fq.gz??? J3_2.fq.gz
cond_SF??? cond_SF1??? SF1_1.fq.gz??? SF1_2.fq.gz
cond_SF??? cond_SF2??? SF2_1.fq.gz??? SF2_2.fq.gz
cond_SF??? cond_SF3??? SF3_1.fq.gz??? SF3_2.fq.gz
cond_SM??? cond_SM1??? SM1_1.fq.gz??? SM1_2.fq.gz
cond_SM??? cond_SM2??? SM2_1.fq.gz??? SM2_2.fq.gz
cond_SM??? cond_SM3??? SM3_1.fq.gz??? SM3_2.fq.gz
二：用trinity自帶的腳本align_and_estimate_abundance.pl進(jìn)行轉(zhuǎn)錄本表達(dá)定量。
$align_and_estimate_abundance.pl? #指明調(diào)用的腳本
--transcripts? /data1/spider/ytbiosoft/data/trinity.all/trinity_out_dir_all.Trinity.fasta? #讀取之前Trinity拼接的轉(zhuǎn)錄本
--seqType fq --samples_file /data/spider_data/tamu/sample.txt? #讀取第一步創(chuàng)建的樣品信息
--est_method? RSEM? #進(jìn)行表達(dá)量計(jì)算的軟件是RSEM
--aln_method bowtie2? #由于RSEM是通過比對(duì)進(jìn)行的表達(dá)量計(jì)算，因此會(huì)采用的bowtie2進(jìn)行比對(duì)
--trinity_mode? #這個(gè)加上會(huì)采用Trinitymode以調(diào)用前期assembly過程中的一個(gè)gene_trans_map文件
--prep_reference? #會(huì)根據(jù)拼接的fasta文件構(gòu)建index?
--output_dir /data1/spider/ytbiosoft/data/RSEMResult? #輸出文件夾，如果不指定路徑的話，由于會(huì)采用讀取樣品信息，因此會(huì)輸入到樣品信息的文件夾
--thread_count 6? #線程數(shù)
三：最后會(huì)生成6個(gè)文件夾（我有6個(gè)樣），每個(gè)文件夾內(nèi)容如下，最主要的是RSEM.genes.results與RSEM.isoforms.results。
├── bowtie2.bam #bowtie2 生成的 bam文件
├── bowtie2.bam.for_rsem.bam #用于RSEM計(jì)算的 bam文件
├── bowtie2.bam.ok
├── RSEM.genes.results #基于基因的EM Reads Count
├── RSEM.isoforms.results #基于轉(zhuǎn)錄本的 EM Reads Count
├── RSEM.isoforms.results.ok
└── RSEM.stat
? ? ? ? ├── RSEM.cnt
? ? ? ? ├── RSEM.model
? ? ? ?└── RSEM.theta
四：由于我們?cè)诿恳粋€(gè)文件夾中的Reads count 沒有經(jīng)樣本間的均一化，因此需要做一個(gè)樣本均一化，構(gòu)建轉(zhuǎn)錄本-基因表達(dá)矩陣并得到不同樣本中的均一化表達(dá)數(shù)據(jù)TMM是后期要做的一個(gè)工作。我們采用以下的矩陣的到了三個(gè)結(jié)果
* 首先手動(dòng)整理到一個(gè)文件夾中:
$find * -name '*.isoforms.results'> quant.file? #這個(gè)地方，我們采用了find命令將子文件夾中的isoforms基因表達(dá)量結(jié)果全部查找出來(lái)然后路徑放到一個(gè)文件中（quant.file）后期要使用
可以$cat quant.file文件看一下內(nèi)容
* 再調(diào)用trinity自帶腳本abundance_estimates_to_matrix.pl構(gòu)建轉(zhuǎn)錄本-基因表達(dá)矩陣
abundance_estimates_to_matrix.pl? #采用的腳本
--est_method RSEM? #由于是對(duì)RSEM的結(jié)果進(jìn)行矩陣構(gòu)建，因此需要說明這個(gè)
--gene_trans_map trinity_out_dir_all.Trinity.fasta.gene_trans_map? #通過這個(gè)map構(gòu)建基因的表達(dá)量矩陣
--name_sample_by_basedir? #**這個(gè)必須要選**，不然會(huì)導(dǎo)致程序沒辦法合并之前的結(jié)果進(jìn)行計(jì)算
--quant_files quant.file? #這個(gè)指明需要的上游文件的位置
* 經(jīng)過計(jì)算后得到的結(jié)果如下：
同理做*.genes.results的表達(dá)矩陣。
最后獲得*.genes.results以及*.isoforms.result的表達(dá)矩陣：
* 以上結(jié)果中分為基因的表達(dá)矩陣和轉(zhuǎn)錄本的表達(dá)矩陣。
其一： '.counts.matrix' 文件用于后期的差異表達(dá)分析
其二： '.TMM.EXPR.matrix'文件可以用于其他基因表達(dá)的分析
到此，就把相關(guān)的基因表達(dá)均一化這邊的工作做完了，應(yīng)該說經(jīng)過這個(gè)工作后所有的基因的表達(dá)量已經(jīng)完成計(jì)算，后期就是通過差異表達(dá)的軟件進(jìn)行一個(gè)差異分析就可以進(jìn)入可視化階段了。
歡迎聯(lián)系互相學(xué)習(xí)：909474045@qq.com
