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