edger和deseq2_转录组分析（二）Hisat2+DESeq2/EdgeR

一、序列比对

在2016年的一篇综述A survey of best practices for RNA-seq data analysis，提到目前有三种RNA数据分析的策略。那个时候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已经被它的作者推荐改用HISAT进行替代。

1. Hisat2教程

1.1 下载安装

#conda直接安装

conda install hisat2

#源码下载安装

wget wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-source.zip

unzip hisat2-2.1.0-source.zip

make

1.2 构建index

直接下载现有的insex或通过Hisat2的方法进行创建

# 其实hisat2-buld在运行的时候也会自己寻找exons和splice_sites，但是先做的目的是为了提高运行效率

extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &

extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf &

# 建立index，必须选项是基因组所在文件路径和输出的前缀

hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19

1.3正式比对

hisat2基本用法就是hisat2 [options]* -x {-1 -2 | -U } [-S ]，基本就是提供index的位置，PE数据或者是SE数据存放位置。然而其他可选参数却是进阶的一大名堂。新手就用默认参数呗。

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

1.4 Hisat2输出结果

比对之后会输出如下结果，解读一下就是全部数据都是100%的，2.88%的配对数据一次都没有比对，94.20%的数据比是唯一比对，2.92%是多个比对。然后如果不按照顺序来，有4.96%的比对。之后把剩下的部分用单端数据进行比对的话，65.57%数据没比对上，33.23%的数据比对一次，1.20%比对超过一次。零零总总的加起来是98.20%的比对。

20182824 reads; of these:

20182824 (100.00%) were paired; of these:

581893 (2.88%) aligned concordantly 0 times

19011569 (94.20%) aligned concordantly exactly 1 time

589362 (2.92%) aligned concordantly >1 times

----

581893 pairs aligned concordantly 0 times; of these:

28886 (4.96%) aligned discordantly 1 time

----

553007 pairs aligned 0 times concordantly or discordantly; of these:

1106014 mates make up the pairs; of these:

725197 (65.57%) aligned 0 times

367552 (33.23%) aligned exactly 1 time

13265 (1.20%) aligned >1 times

98.20% overall alignment rate

2. SAMtools三板斧

SAM(sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式，当然他可以用于存放未比对的数据。而目前处理SAM格式的工具主要是SAMTools

view: BAM-SAM/SAM-BAM 转换和提取部分比对

sort: 比对排序

merge: 聚合多个排序比对

index: 索引排序比对

faidx: 建立FASTA索引，提取部分序列

tview: 文本格式查看序列

#最常用的三板斧就是格式转换，排序，索引。

samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam

samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam

samtools index SRR35899${i}_sorted.bam

3. BAM/SAM文件格式

SAM文件主要由两个部分构成:

header：标记了该SAM文件的一些基本信息，比如版本、按照什么方式排序的、Reference信息等等。

本体：每行为一个reads，不同列记录了不同的信息，列与列之间通过tab分隔。

每列的含义：

MAPQ值：

表示为mapping的质量值，mapping Quality, It equals -10log10Pr{mapping position is wrong}, rounded to the nearest integer, A value 255 indicates that the mapping quality is not available. 该值的计算方法是mapping的错误率的-10log10值，之后四舍五入得到的整数，如果值为255表示mapping值是不可用的，如果是unmapped read则MAPQ为0，一般在使用bwa mem或bwa aln(bwa 0.7.12-r1039版本)生成的sam文件，第五列为60表示mapping率最高，一般结果是这一列的数值是从0到60，且0和60这两个数字出现次数最多

想把小于2的都丢弃：

samtools view -bSq 2 file.sam > filtered.bam

flag的含义：

1 ：代表这个序列采用的是PE双端测序

2：代表这个序列和参考序列完全匹配，没有插入缺失

4：代表这个序列没有mapping到参考序列上

8：代表这个序列的另一端序列没有比对到参考序列上，比如这条序列是R1,它对应的R2端序列没有比对到参考序列上

16：代表这个序列比对到参考序列的负链上

32 ：代表这个序列对应的另一端序列比对到参考序列的负链上

64 ：代表这个序列是R1端序列， read1;

128 : 代表这个序列是R2端序列，read2；

256：代表这个序列不是主要的比对，一条序列可能比对到参考序列的多个位置，只有一个是首要的比对位置，其他都是次要的

512：代表这个序列在QC时失败了，被过滤不掉了(# 这个标签不常用)

1024: 代表这个序列是PCR重复序列(#这个标签不常用)

2048: 代表这个序列是补充的比对(#这个标签具体什么意思，没搞清楚，但是不常用)

cigar的含义:

cigar中会包含数字，代表了特定match持续了多少nt；以及不同的字符，代表了不同的match情况。

30S512M216N12S (30nt soft clip -> 512nt exact match -> 216nt skipped region -> 12nt soft clip)

30S (30nt soft clip)

40M (40nt exact match)

其中不同的字符及其含义如下：

参考：

https://www.jianshu.com/p/a584d31418f3

https://www.jianshu.com/p/9c87bba244d8

二、htseq-count的使用

HTSeq作为一款可以处理高通量数据的python包，由Simon Anders, Paul Theodor Pyl, Wolfgang Huber等人携手推出HTSeq — A Python framework to work with high-throughput sequencing data。自发布以来就备受广大分析人员青睐，其提供了许多功能给那些熟悉python的大佬们去自信修改使用，同时也兼顾着给小白们提供了两个可以拿来可用的可执行文件 htseq-count(计数) 和 htseq-qa(质量分析)。

具体参考：

https://www.cnblogs.com/triple-y/p/9338890.html

https://blog.csdn.net/herokoking/article/details/78257714

三、基因差异表达分析

1. DESeq2(DESeq2不支持无生物学重复的数据)

library("DESeq2")

#directory