bulk RNA-Seq （2）比对到参考基因组

欢迎关注Bioinfor 生信云！

上一篇推文讲了bulk RNA-Seq 的数据清洗部分，接下来就是将清洗后的数据比对到参考基因组上（有参）。如果没有参考基因组，那么就需要我们自己去组装基因组了。

1.mapping常用的工具

1.1BWA：包括BWA-backtrack, BWA-SW and BWA-MEM，第一个算法设计用于 Illumina 序列读取高达 100bp，而其余两个用于更长的序列，范围从 70bp 到 1Mbp。对于 70-100bp Illumina 读取，BWA-MEM 也具有比 BWA-backtrack 更好的性能。

1.2Bowtie2：它特别擅长将大约 50 到 100 个字符的读取与相对较长的（e.g. mammalian）基因组对齐。Bowtie 2 通常是比较基因组学的第一步，包括call snp、ChIP-seq、RNA-seq、BS-seq。

1.3Hisat2：Hisat2 是一种快速灵敏的比对程序，用于将下一代测序读数（DNA 和 RNA）映射到人类基因组群以及单个参考基因组，主要用于RNA-seq数据。

这里我选择hisat2做比对数据处理。

2.参考序列的准备

2.1需要准备三个文件：

基因组序列（genome.fa）

基因注释文件（genes.gtf）

蛋白序列文件（proteins.fa）

2.2常用的参考基因组下载地址

NCBI：https://www.ncbi.nlm.nih.gov/

Ensembl：https://asia.ensembl.org/index.html

EnsemblPlants：http://plants.ensembl.org/index.html

UCSC：http://genome.ucsc.edu/

phytozome：https://phytozome-next.jgi.doe.gov/

还有一些物种有专有的基因组网站

3.数据预处理

3.1基因组序列一般可以直接使用

3.2注释文件可能会下载到gff格式的，需要用gffread做转换

gffread -T -o genes.gtf genes.gff

3.3蛋白序列一般也是可以直接使用的（有时候需要注意基因的编号是否一致，以及是否有可变剪切）

4.比对到参考基因组
bulk RNA-Seq 不包含内含子，所以在重新比对回参考基因组时，会被中间的内含子隔开，这种情况叫做splice-
alignment。
4.1软件安装（linux）

conda install hisat2
conda install samtools

4.2给参考基因组构建Index

hisat2-build genome.fa \ #参考基因组
./index/genome \ #index名称
1>hisat2-build.log 2>&1 #将日志信息存在.log文件里面

4.3hisat2比对

hisat2 --new-summary \#使用新的日志格式
-p 4 \#线程数
-1 ../clean_data/CER3_1_R1.clean.fastq.gz \#质控过滤后的文件，双端测序 -1 -2 来指示
-2 ../clean_data/CER3_1_R2.clean.fastq.gz \
2>CER3_1.mapping.log | samtools view -o CER3_1.bam #比对完的结果通过管道符交给samtools转为bam文件，这一步是为了节省空间

输入：参考基因组的index文件，每个样本的测序数据
输出：bam格式的比对结果
可以使用awk，或者for循环批量处理

4.4排序及bam文件构建index

samtools sort -@ 2 -o CER3_1.s.bam CER3_1.bam #bam文件排序
samtools index CER3_1.s.bam #bam文件构建index用于IGV查看比对结果
samtools flagstat CER3_1.s.bam >CER3_1.s.metrics.txt #保存日志

5.使用IGV查看比对结果

IGV是broad institute 开发的数据可视化软件，用来可视化sam/bam文件。https://software.broadinstitute.org/software/igv/