转录组分析（三）Hisat2+StringTie+Ballgow

一、分析流程

该分析流程主要根据2016年发表在Nature Protocols上的一篇名为Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown的文章撰写的，主要用到以下三个软件：

1. HISAT (http://ccb.jhu.edu/software/hisat/index.shtml)利用大量FM索引，以覆盖整个基因组，能够将RNA-Seq的读取与基因组进行快速比对，相较于STAR、Tophat，该软件比对速度快，占用内存少。
2. StringTie(http://ccb.jhu.edu/software/stringtie/)能够应用流神经网络算法和可选的de novo组装进行转录本组装并预计表达水平。与Cufflinks等程序相比，StringTie实现了更完整、更准确的基因重建，并更好地预测了表达水平。
3. Ballgown (https://github.com/alyssafrazee/ballgown)是R语言中基因差异表达分析的工具，能利用RNA-Seq实验的数据(StringTie, RSEM, Cufflinks)的结果预测基因、转录本的差异表达。然而Ballgown并没有不能很好地检测差异外显子，而 DEXseq、rMATS和MISO可以很好解决该问题。

1、使用HISAT将读段匹配到参考基因组上,使用者可以提供注释文件，但HISAT依旧会检测注释文件没有列出来的剪切位点。
2、比对上的reads将会被呈递给StringTie进行转录本组装，StringTie单独的对每个样本进行组装，在组装的过程中顺带估算每个基因及isoform的表达水平。
3、所有的转录本都被呈递给StringTie的merge函数进行merge，这一步是必须的，因为有些样本的转录本可能仅仅被部分reads覆盖，无法被第二步的StringTie组装出来。merge步骤可以创建出所有样本里面都有的转录本，方便下一步的对比。
4、merge的数据再一次被呈递给StringTie，StringTie可以利用merge的数据重新估算转录本的丰度，还能额外的提供转录本reads数量的数据给下一步的ballgown。
5、Ballgown从上一步获得所有转录本及其丰度，根据实验条件进行分类统计。

1. hisat2比对：

# 建立index， 必须选项是基因组所在文件路径和输出的前缀，其实hisat2-buld在运行的时候也会自己寻找exons和splice_sites
hisat2-build Oryza_sativa.IRGSP-1.0.genome 
hisat2 --dta  -p 6 --max-intronlen 5000000 -x Oryza_sativa.IRGSP-1.0.genome -1 C1-1_good_1.fq -2 C1-1_good_2.fq -S C1-1.HISAT_aln.sam  >hisat2_running.log 2>&1

2. samtools

最终我们获得了12个sam文件（此处只列举其中一个样本）然后通过samtools将sam文件转换为bam文件，作为stringtie的输入文件，具体脚本为：

samtools view -F 4 -Su C1-1.HISAT_aln.sam  | samtools sort -T C1-1.accepted_hits -o C1-1.accepted_hits.bam  &&   samtools index C1-1.accepted_hits.bam.bai

3. StringTie组装、预测新基因

3.1 组装
接下来利用stringtie对转录组进行组装，会针对每个bam文件生成一个gtf文件，它主要记录了转录本的组装信息

stringtie C1-1.accepted_hits.bam -G Oryza_sativa.IRGSP-1.0.gtf -l C1-1 -o C1-1.transcripts.gtf

3.2 合并
然后利用软件stringtie将12个含有转录本信息的gtf文件合并成一个gtf，此时需要预先将12个GTF文件的文件名录入到mergelist.txt文件中，下载的数据中已经给出该文件，执行完会多出一个GTF文件，即tringtie_merged.gtf：

stringtie --merge -G Oryza_sativa.IRGSP-1.0.gtf -F 0.1 -T 0.1 -i -o StringTie_merged.gtf mergelist.txt

3.3 gffcompare—与已知的转录本注释文件进行比较
比较不同样本的转录本定量信息需要先将转录本信息储存为相同的格式，一般组装软件的输出结果都是gtf或gff。由于在组装的过程中产生了大量的新的转录本信息，而我们仅通过肉眼观察其唯一的注释信息----染色体上的起始位置，很显然无法阐明其中蕴含的生物学意义，因此我们需要将它们与已知的转录本注释文件---annotation.gtf进行比较，将新得到的转录本与注释好的转录本之间建立联系，这样可以让我们更好地发现新的转录本。

gffcompare -r Oryza_sativa.IRGSP-1.0.gtf StringTie_merged.gtf

输出文件六个，前四个文件可以指定保存位置，后两个文件是跟输入的gtf文件保存在一个位置，并且都是以-o提供的前缀开头的

gffcmp.annotated.gtf：包含了class code信息，该文件一般用于下文继续stringtie
gffcmp.stats：包含了feature的统计信息，也包含了找到新的外显子、内含子的数目，其中有两个统计量sensitivity和precision，定义为 Sensitivity is defned as the proportion of genes from the annotation that are correctly reconstructed,whereas precision (also known as positive predictive value) captures the proportion of the output that overlaps the annotation
gffcompare.loci：见说明书
gffcompare.tracking：见说明书
gffcompare_result.refmap：这个文件包含四列信息，第一列ref_gene_id是gene symbol ，无symbol的给出的是ensemble的gene id; 第二列ref_id是指ensemble的transcript id; 第三列class_code 是“=”和“c”;第四列是cuff_id_list。这个文件指组装后与参考基因组几乎完全匹配的转录本
gffcompare_result.tmap:包含了转录本的定量信息，如cov，FPKM等，可用于定量或筛选新转录本

3.4 筛选新基因
那我们需要筛选出新的转录本，那该如何筛呢？这个可以从GTF文件的class codes着手，该信息记录了每个转录本相对于已知转录本的位置信息。通过这个class_code 我们一般选在3种类型的转录本，分别是:
i : 内含子区的转录本
u: 基因间区的新转录本
x: 已知外显子的反义链转录本
具体代码实现：

#perl get_track.pl -gtf  gffcmp.annotated.gtf -index rice  -out rice.newGene.track.list.info
$gtf = abs_path($gtf);
my %newGene_track;
open IN,"$gtf" || die $!;
while (<IN>) {
    chomp;
    next if (/^$/ || /^\#/);
    my @tmp = split/\t+/,$_;
    my @info = split/;/,$tmp[8];
    my ($gene_track,$geneID,$iso_track,$isoID,$class_code);
    if ($tmp[8]=~/class_code \"u\";/) {
        $tmp[8]=~/transcript_id\s\"([^\"]+)\";\sgene_id\s\"([^\"]+)\";\sxloc/;
        $gene_track=$2;
        $iso_track=$1;
        $newGene_track{$gene_track}{$iso_track}=1;
        next;
    }
}
close IN;
open OUT,">$out" || die $!;
my $i=1;
foreach my $gene_track (sort keys %newGene_track) {
    my $gene_name = "$index"."_newGene"."_$i";
    my $j = 1;
    foreach my $iso_track (sort keys %{$newGene_track{$gene_track}}) {
        my $iso_name = "$gene_name".".$j";
        print OUT "$gene_track\t$gene_name\t$iso_track\t$iso_name\n";
        $j++;
    }
    $i++;
}
close OUT;

3.5 转换格式
通过读取gffcmp.annotated.gtf文件对应生成新基因的gff或gtf格式文件rice.newGene_final.gff
并且读取参考序列将gff文件转为对应的转录本序列rice.newGene.transcript.fa
3.6 TransDecoder预测CDS
识别转录本序列中的潜在的编码区域，也就是预测CDS。具体教程参考

https://www.bioinfo-scrounger.com/archives/106/

#提取FASTA序列中最长的ORF
 TransDecoder.LongOrfs -t rice.newGene.transcript.fa -O longest_orfs.pep

3.7 进行过滤
通过判断蛋白序列中ORF长度，设定50的阈值过滤掉那些较短的编码蛋白新基因，降低假阳性。（随着ORF长度的减少，其假阳性的概率也会随之上升）
3.8 merge known and new gff

cat rice.newGene_final.filtered.gff Oryza_sativa.IRGSP-1.0.gff >final.gff 

4. Ballgown差异表达分析

4.1 Ballgown文件准备

stringtie C1-1.accepted_hits.bam -eB -G final.gff   -o ./tringTie_asm.gtf -A ./gene_abundence.tab -C ./known.cov_refs.gtf

参数说明：

-e  限制reads比对的处理，仅估计和输出与用-G选项给出的参考转录本匹配的组装转录本。使用该选项，则会跳过处理与参考转录本不匹配的组装转录本，这将大大的提升了处理速度。
-B  应用该选项，则会输出Ballgown输入表文件（* .ctab），其中包含用-G选项给出的参考转录本的覆盖率数据。（有关这些文件的说明，请参阅Ballgown文档。）
-o [<path/>]<out.gtf> 设置StringTie组装转录本的输出GTF文件的路径和文件名。
-G <ref_ann.gff>    使用参考注释基因文件指导组装过程，格式GTF/GFF3。输出文件中既包含已知表达的转录本，也包含新的转录本。选项-B，-b，-e，-C需要此选项（详情如下）
-A <gene_abund.tab> 输出基因丰度的文件（Coverage, FPKM, TPM）
-C <cov_refs.gtf>   输出所有转录本对应的reads覆盖度的文件，此处的转录本是指参考注释基因文件中提供的转录本。(需要参数 -G).

由Stringtie生成的Ballgown可读的表达文件如下：

e_data.ctab: 外显子水平表达值
i_data.ctab:内含子水平表达值
t_data.ctab:转录组水平表达值
e2t.ctab:表中有两列，e_id和t_id，表示哪些外显子属于哪些转录本。这些id与e_data和t_data表中的id匹配。
i2t.ctab:表中有两列，i_id和t_id，表示哪些内含子属于哪些转录本。这些id与i_data和t_data表中的id匹配。

4.2 差异分析

参考Ballgown官网说明
https://github.com/alyssafrazee/ballgown