mRNA-seq学习（五）：基于比对的转录本组装

Cufflinks基于reads与参考基因组的比对来进行全长转录本的组装，不需要额外提供gff/gtf这种描述gene/transcript的feature文件，软件会生成gff文件。

1. 比对

第一步仍然是建索引比对，这里用到tophat2（以bowtie2为搜索引擎）。

bowtie2-build -f Arabidopsis_thaliana.TAIR10.dna.toplevel.fa TAIR10.bybowtie
mv Arabidopsis_thaliana.TAIR10.dna.toplevel.fa TAIR10.bybowtie.fa #和一些比对软件（比如（二）中介绍的bowtie2）不同的是，tophat比对需要参考基因组跟索引有相同的前缀；而像bowtie2只需要利用索引即可，所以没有这个要求
$ ls
TAIR10.bybowtie.1.bt2  TAIR10.bybowtie.2.bt2  TAIR10.bybowtie.3.bt2  TAIR10.bybowtie.4.bt2  TAIR10.bybowtie.fa  TAIR10.bybowtie.rev.1.bt2  TAIR10.bybowtie.rev.2.bt2

for i in `seq 2 7`
do
tophat2 -r 230 -p 8 -o ./SRR328680${i} ~/learn_rnaseq/mRNA-seq/bowtie2/ref/TAIR10.bybowtie \ #r表示insert size
 ~/learn_rnaseq/mRNA-seq/data/SRR328680${i}_1.fastq.gz \
 ~/learn_rnaseq/mRNA-seq/data/SRR328680${i}_2.fastq.gz
done

#输出文件，其中accepted_hits.bam后面会用到
$ ls -t
unmapped.bam  logs/  accepted_hits.bam  junctions.bed  deletions.bed  insertions.bed  align_summary.txt  prep_reads.info

2. 组装

2.1 官方流程

Cufflinks的流程图，可以看出当有多个处理（或重复）时，流程建议分开跑tophat和cufflinks，然后再把分别组装得到的转录本合并。

Cufflinks在组装的时候可以提供额外的gff文件用于指导组装过程。

for i in {2..7}
do
cd ~/learn_rnaseq/mRNA-seq/map/map_tophat2/SRR328680${i}
cufflinks -p 2 -o ./creat/ ./accepted_hits.bam &
done

之后每个creat文件夹下面会生成4个文件

genes.fpkm_tracking  isoforms.fpkm_tracking  skipped.gtf  transcripts.gtf

接下来将这6个样本的transcripts.gtf合并

$ cat list
SRR3286802_transcripts.gtf
SRR3286803_transcripts.gtf
SRR3286804_transcripts.gtf
SRR3286805_transcripts.gtf
SRR3286806_transcripts.gtf
SRR3286807_transcripts.gtf

nohup cuffmerge -p 8 list -s ~/learn_rnaseq/mRNA-seq/bowtie2/ref/TAIR10.bybowtie.fa -g ~/learn_rnaseq/mRNA-seq/ref/new.GCF_000001735.4_TAIR10.1.gtf &

#之后生成一个merged_asm文件夹，merged.gtf为合并后的转录本gtf文件
$ ls merged_asm/
logs/  merged.gtf

之后做差异表达，结果相当丰富

nohup cuffdiff ./merged_asm/merged.gtf \
 SRR3286802_accepted_hits.bam,SRR3286803_accepted_hits.bam,SRR3286804_accepted_hits.bam \
 SRR3286805_accepted_hits.bam,SRR3286806_accepted_hits.bam,SRR3286807_accepted_hits.bam &