Trimmomatic，Hisat2，SAMtools，Stringtie/Htseq，Deseq
首先先下个conda，便于下载某些软件：
1.先将这个链接中的命令跟着走一遍（安装的时候如果失败/出错，重新试一次就可以了，它经常出错）：用conda安装RNA-seq所需要的工具
2.安装软件
conda search xxx
conda install xxx
3.使用：刚安装完的时候直接输入软件名就可以直接打开，但重启终端之后需要先启动环境
source activate
conda activate python2
一、数据质控：Trimmomatic
1.下载安装：
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip
unzip Trimmomatic-0.38.zip
2.输入文件：sample_R1.fastq、sample_R2.fastq
（java -jar <path to trimmomatic.jar> ）PE [-threads <threads] [-phred33 | -phred64] [-trimlog<logFile>] >] [-basein <inputBase> | <input 1> <input 2>] [-baseout <outputBase> |<paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> <step 2> ...
-threads 设置多线程运行数 ，-phred可不选
产生4个结果文件
参考：高通量测序数据质控神器Trimmomatic、Trimmomatic下载安装、Trimmomatic基本功能的使用
二、数据比对：HISAT2
需要用到的数据：基因组序列，转录组序列；结果得到SAM文件
基因组比对软件常用bwa，转录组比对软件常用bowtie2、hisat2等，其中有参考基因组的常用hisat2，没有参考基因组的常用bowtie2。
此处是要将转录组比对到基因组上去。
1.可以使用conda下载：
conda install hisat2
PATH=/home/freshman/lbq/hisat2:$PATH
PATH=/home/freshman/lbq/hisat2/*:$PATH
2.建立基因组、转录组索引：生成的文件会默认建立在你现在所在的目录下
（hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19）
hisat2-build –p 4 ~/data/cdna.fasta genome
hisat2-build –p 4 ~/data/dna.fasta gdna
3.对比命令及结果如下：
hisat2 -f -p16 -x ~/data/gdna -U cdna.fasta -S gly.sam
三、将SAM转换为BAM文件：SAMtools
需要HISAT产生的SAM文件，转成BAM文件后对BAM文件进行分析
1.使用conda下载samtools
conda install samtools
2.利用将sam文件转为bam文件（samtools同样可以用conda下载）
samtools view -S gly.sam -b > gly.bam
3.查看bam文件
cd data
samtools flagstat gly.bam
参考：SAMTools
四、重复序列注释：StringTie
需要sort过的BAM文件与GTF文件
参考：Stringtie的使用、用法详解
（四五二选一）
五：表达量分析——Htseq
输入文件：SAM和GTF，其中SAM需要samtools重新sort
1.下载安装
wget https://pypi.python.org/packages/46/f7/6105848893b1d280692eac4f4f3c08ed7f424cec636aeda66b50bbcf217e/HTSeq-0.7.2.tar.gz
tar zxf HTSeq-0.7.2.tar.gz
cd HTSeq-0.7.2
python setup.py install --user
vim .bashrc
PATH=/home/freshman/lbq/HTSeq-0.7.2/.local/bin:$PATH
source .bashrc
chmod 777 /home/freshman/lbq/HTSeq-0.7.2/python2/scripts/htseq-count
export PATH=$PATH:/home/freshman/lbq/HTSeq-0.7.2/*
2.计数
htseq-count -f bam -r name -s no -a 10 -t exon -i ID -m union ~/data/gly_sort.bam ~/data/ori/genes.gff3 > ~/data/counts_out.txt
# 命令参数
-f | --format default: sam 设置输入文件的格式，该参数的值可以是sam或bam。
-r | --order default: name 设置sam或bam文件的排序方式，该参数的值可以是name或pos。前者表示按read名进行排序，后者表示按比对的参考基因组位置进行排序。若测序数据是双末端测序，当输入sam/bam文件是按pos方式排序的时候，两端reads的比对结果在sam/bam文件中一般不是紧邻的两行，程序会将reads对的第一个比对结果放入内存，直到读取到另一端read的比对结果。因此，选择pos可能会导致程序使用较多的内存，它也适合于未排序的sam/bam文件。而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上，也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的，但HTSeq推荐使用name排序，且一般比对软件的默认输出结果也是按name进行排序的。
-s | --stranded default: yes 设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表示非链特异性测序；若是单端测序，yes表示read比对到了基因的正义链上；若是双末端测序，yes表示read1比对到了基因正义链上，read2比对到基因负义链上；reverse表示双末端测序情况下与yes值相反的结果。根据说明文件的理解，一般情况下双末端链特异性测序，该参数的值应该选择reverse（本人暂时没有测试该参数）。
-a | --a default: 10 忽略比对质量低于此值的比对结果。在0.5.4版本以前该参数默认值是0。
-t | --type default: exon 程序会对该指定的feature（gtf/gff文件第三列）进行表达量计算，而gtf/gff文件中其它的feature都会被忽略。
-i | --idattr default: gene_id 设置feature ID是由gtf/gff文件第9列那个标签决定的；若gtf/gff文件多行具有相同的feature ID，则它们来自同一个feature，程序会计算这些features的表达量之和赋给相应的feature ID。
-m | --mode default: union 设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知，对于原核生物，推荐使用intersection-strict模式；对于真核生物，推荐使用union模式。
-o | --samout 输出一个sam文件，该sam文件的比对结果中多了一个XF标签，表示该read比对到了某个feature上。
-q | --quiet 不输出程序运行的状态信息和警告信息。
-h | --help 输出帮助信息。
其中，gene_id这一处如果报错的话，less genes.gff3看看，它那一列的名字是不是叫其他的，改了就行。
结果如下：其中ambiguous表示该read比对到多个gene上；no_feature表示read没有比对到gene上
参考：RNA-seq分析htseq-count的使用、使用HTSeq进行有参转录组的表达量计算 | 陈连福的生信博客
六、差异表达分析——Deseq（R包）
输入文件：read count
1.安装：（此处R版本为3.6.1）
R
install.packages("BiocManager")
BiocManager::install(c("DESeq2"))
2.
网友评论