1. 数据

手头有两个转录组数据，分别是某一真菌物种(该物种没有基因组数据)的双核菌丝阶段和子实体阶段的转录组数据。测序平台是Illumina Hiseq2000；插入片段长度200bp，测序的reads长度90bp。得到的数据文件为：

/home/user/RNA-seq/mycelium_reads1.fastq
/home/user/RNA-seq/mycelium_reads2.fastq
/home/user/RNA-seq/fruitingbody_reads1.fastq
/home/user/RNA-seq/fruitingbody_reads2.fastq

2. 数据的预处理

2.1. 去除N的比例大于5%的reads;去除低质量reads（质量值Q≤20的碱基数占整个read的50％以上）；
2.2. 根据FastQC对上述过滤后的reads的质量检测，去除reads首尾各10bp的碱基，得到的预处理数据为：

/home/user/RNA-seq/clean_reads/mycelium_reads1.fastq
/home/user/RNA-seq/clean_reads/mycelium_reads2.fastq
/home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq
/home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq

3. 使用Trinity和TGICL进行转录组的组装

3.1. 对mycelium数据和fruitingbody数据分别进行转录组的组装

$ pwd
/home/user/RNA-seq/

$ mkdir assembly
$ cd assembly

$ Trinity.pl --jaccard_clip --seqType fq --JM 50G --SS_lib_type FR --CPU 24 --inchworm_cpu 24 --bflyCPU 24 --group_pairs_distance 500 -min_contig_length 200 --output mycelium_contig --left /home/user/RNA-seq/clean_reads/mycelium_reads1.fastq --right /home/user/RNA-seq/clean_reads/mycelium_reads2.fastq
$ Trinity.pl --jaccard_clip --seqType fq --JM 50G --SS_lib_type FR --CPU 24 --inchworm_cpu 24 --bflyCPU 24 --group_pairs_distance 500 -min_contig_length 200 --output fruitingbody_contig --left /home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq --right /home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq

当然，需要将生成的两个转录组的Trinity.fasta序列按长度进行排序；统一更改的fasta
头名称；更改fasta文件名

$ cp mycelium_contig/Trinity.fasta mycelium_contigs.fasta
$ cp fruitingbody_contig/Trinity.fasta fruitingbody_contigs.fasta

3.2. 使用TGICL将两个转录组序列合并

$ pwd
/home/user/RNA-seq/assembly

$ mkdir all_tissue_contig
$ cd all_tissue_contig

$ cat ../mycelium_contigs.fasta ../fruitingbody_contigs.fasta > all.contigs.fasta
$ tgicl -F all.contigs.fasta

tgicl生成了两个有用的文件： asm_1/contigs  和 all.contigs.fasta.single
tons。其中前者是聚类后的contigs结果；后者是没有聚类的单独的contigs的序列名，需
要分别到../mycelium_contigs.fasta 和 ../fruitingbody_contigs.fasta文
件中提取出相应的序列: all.contigs.fasta.singletons.mycelium.fasta 和 
all.contigs.fasta.singletons.fruitingbody.fasta

$ cat asm_1/contigs all.contigs.fasta.singletons.mycelium.fasta all.contigs.fasta.singletons.fruitingbody.fasta > all_contigs.fasta

在对all_contigs.fasta进行fasta头的重命名，并将序列按长度排序；同时要得到all_
contigs.fasta中的序列和mycelium_contigs.fasta，fruitingbody_contigs.
fasta中序列的对应关系。

$ cp all_contigs.fasta ../

3.3. 至此得出转录组的组装结果：

/home/user/RNA-seq/assembly/all_contigs.fasta

4. 使用cufflinks来分析差异表达基因

4.1 使用tophat来将预处理后的reads比对到转录组序列上

$ pwd
/home/user/RNA-seq/

$ mkdir cufflinks
$ cd cufflinks

$ bowtie2-build ../all_contigs.fasta all_contigs

$ tophat -o tophat_out_mycelium -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean_reads/mycelium_reads1.fastq /home/user/RNA-seq/clean_reads/mycelium_reads2.fastq
$ tophat -o tophat_out_fruitingbody -r 60 --mate-std-dev 80 -p 24 --library-type fr-unstranded all.contigs /home/user/RNA-seq/clean_reads/fruitingbody_reads1.fastq /home/user/RNA-seq/clean_reads/fruitingbody_reads2.fastq

4.2 自制一个transcripts.gtf文件
由于cuffdiff运行需要一个参考序列的transcripts.gtf文件: 该文件有9列，使用tab分隔；使exon的范围为整个contig。其格式如：

AA.auricula_all_contig_1     chenlianfu  exon  1  9401  .  .  .  gene_id "A.auricula_all_contig_1"; transcript_id "A.auricula_all_contig_1";
A.auricula_all_contig_10     chenlianfu  exon  1  4464  .  .  .  gene_id "A.auricula_all_contig_10"; transcript_id "A.auricula_all_contig_10";
A.auricula_all_contig_100    chenlianfu  exon  1  3090  .  .  .  gene_id "A.auricula_all_contig_100"; transcript_id "A.auricula_all_contig_100";
A.auricula_all_contig_1000   chenlianfu  exon  1  1768  .  .  .  gene_id "A.auricula_all_contig_1000"; transcript_id "A.auricula_all_contig_1000";
A.auricula_all_contig_10000  chenlianfu  exon  1  586   .  .  .  gene_id "A.auricula_all_contig_10000"; transcript_id "A.auricula_all_contig_10000"; 
A.auricula_all_contig_10001  chenlianfu  exon  1  586   .  .  .  gene_id "A.auricula_all_contig_10001"; transcript_id "A.auricula_all_contig_10001";

4.3 使用cuffdiff来分析转录子的表达量和差异表达基因

$ pwd
/home/user/RNA-seq/cufflinks

$ cuffdiff -L mycelium,fruitingbody --library-type fr-unstranded -p 8 -o cuffdiff ./transcriptome.gtf ./tophat_out_mycelium/accepted_hits.bam ./tophat_out_fruitingbody/accepted_hits.bam

至此得出转录子的表达量数据和差异表达分析

/home/user/RNA-seq/cufflinks/cuffdiff/gene_exp.diff