生信步骤|转录组测序上游分析：hisat2+samtools+stringtie

学术程稻属

已于 2022-02-12 12:15:58 修改

阅读量8.7k

点赞数 10

分类专栏：生信步骤文章标签： linux

于 2022-01-12 12:53:49 首次发布

本文链接：https://blog.csdn.net/m0_64240043/article/details/122450444

版权

生信步骤专栏收录该内容

7 篇文章

订阅专栏

转录组分析在当下研究功能基因组领域十分常用。相关软件组合种类也十分丰富，本文采用了hisat2+samtools+stringtie策略从转录组数据中挖掘差异表达基因。在这里小编整理了一下此套组合的执行流程，以供日后查阅；同时分享在平台，希望能帮助到更多初学者，如有谬误也请各路大佬批评指正。

先从整体上看一下软件们所执行的功能：
hisat2：建立参考基因组索引，reads的比对
samtools：sam2bam的转化
stringtie：估算转录本表达量

所用的数据结构如下，此处用了酵母的转录组数据和参考基因组：
双端测序数据已经经过fastqc过滤，具体过滤流程本文没有涉及。示例数据仅供参考。执行时要关注文件格式转化，以及各种格式下包含的生物信息。

@biocloud:~/1223/NGS2022$ tree
.
├── gene_data.csv# 差异表达基因样例，相当于最终输出的标准答案。
├── genome
│   ├── yeast.fa# 参考基因组文件
│   ├── yeast.gff# 参考基因组注释文件
│   └── yeast_transcriptome.fa # Kallisto 所需的转录组索引文件，实现基于 pseudo alignment 的转录本定量时才会用到，本文不涉及该文件。
├── reads	# 双端测序clean data，即已经经历过reads过滤。通常下机的时候公司都会同时给出rawdata和cleandata，直接用后者就好。
│   ├── s1_y_1.fq.gz
│   ├── s1_y_2.fq.gz
│   ├── s2_y_1.fq.gz
│   └── s2_y_2.fq.gz
├── script	 # 脚本文件，需要时加入绝对路径引用该脚本即可。
│   ├── edgeR.R
│   ├── prepDE.py
│   ├── prepDE.py3
│   └── run.sh
└── src		# 可能用到的软件安装包，服务器如果安装过该软件则无需再次安装。
    ├── fastqc_v0.11.9.zip
    ├── hisat2-2.2.1-Linux_x86_64.zip
    ├── hisat2-2.2.1-OSX_x86_64.zip
    ├── kallisto_linux-v0.46.1.tar.gz
    ├── kallisto_mac-v0.46.1.tar.gz
    ├── samtools-1.14.tar.bz2
    ├── stringtie-2.2.0.Linux_x86_64.tar.gz
    └── stringtie-2.2.0.OSX_x86_64.tar.gz

1.hisat2构建参考基因组索引

# 进入存储有参考基因组yeast.fa的目录genome 
$ cd ./genome 
# 将参考基因组yeast.fa建立索引并重命名为genome.fa。注意hisat2与-build之间不要加空格。
$ hisat2-build yeast.fa genome.fa

2.hisat2将同一个处理下的双端测序reads比对到参考基因组

# 重新新建文件夹alignment并在其中执行比对操作
$ cd ..
$ mkdir alignment
$ cd alignment

# hisat2将同一个处理下的测序得到的双端测序reads比对到建立好索引的参考基因组。
# hisat2 参数-1和-2后面跟着双端测序数据（.gz格式）路径。比对后的结果保存为s1.sam。
$ hisat2 -p 6 -x ../genome/genome.fa -1 ../reads/s1_y_1.fq.gz -2 ../reads/s1_y_2.fq.gz -S ./s1.sam

3.samtools将sam文件进行二进制转化，排序以及添加索引。

# samtools view将刚刚得到的s1.sam转化为s1.bam
$ samtools view -bS s1.sam -o s1.bam //bam是sam的二进制文件，二进制转化后占用存储空间更小。

# samtools sort将s1.bam排序，产生s1.sorted.bam。后者会自动加上.bam后缀，无需命令行里添加。
$ samtools sort s1.bam s1.sorted 

# samtools index将排序后的bam文件添加索引
$ samtools index s1.sorted.bam

4.stringtie根据gtf注释和排序后bam比对结果估算转录本表达量

# 输入刚刚排序好的s1.sorted.bam文件，酵母基因组注释文件yeast.gff。
# 运行stringtie得到gtf文件(s1_out.gtf)和基因表达列表(s1_genes.list)。
$ stringtie ./s1.sorted.bam -G ../genome/yeast.gff -e -p 2 -o ./s1_out.gtf -A ./s1_genes.list

用同样的方法得到另一处理下双端测序的比对结果：s2_out.gtf 和 s2_genes.list。步骤与上述完全一致，此处不再逐一重复。每个处理下的一对双端测序结果会产生一个s1_out.gtf文件。按照上述步骤处理完另一对双端测序后应该得到两个gtf文件：s1_out.gtf 和 s2_out.gtf。

5.汇总两个处理下的转录组表达量

# 在新建的文件夹执行汇总表达量操作
$ cd ..
$ mkdir differential_expression
$ cd differential_expression

#手动新建并编辑一个sample_list.txt文件，用以存放处理的名称和上述得到的gtf文件路径。
$ vim sample_list.txt

sample_list.txt文件格式如下所示：两个处理的名称（s1，s2）+对应gtf文件的地址。注意第二行末尾不要有换行符号（\n）！！！处理名称（s1，s2）和gtf文件地址之间应该有一个空格。

s1 /mnt/1223/NGS2021/differential_expression/s1_out.gtf
s2 /mnt/1223/NGS2021/differential_expression/s2_out.gtf

6.使用stringtie的prepDE.py3脚本生成差异表达基因列表。

prepDE.py3脚本可以在github中下载（https://github.com/gpertea/stringtie）

# 注意stringtie差异基因汇总脚本有prepDE.py和prepDE.py3两个版本
# 前者适用于python2环境，后者python3环境。使用混了会报错。本文基于python3.9.5环境。
$ python prepDE.py3 -i sample_list.txt -g gene_count.csv -t transcript.csv

执行完上述步骤后得到的gene_count.csv即为差异表达基因列表，可以导入到R语言用edgeR / DESeq2包进行差异表达基因分析，有机会再整理后继教程。

sam文件基础

sam文件在序列分析中至关重要，无论是转录组分析还是基因组call SNP。认识sam文件所包含的信息有助于理解数据。我们依然以本文中生成sam文件的命令行举例。转录组分析核心比对步骤是使用hisat2将测序得到的reads比对到建立好索引的参考基因组：

$ hisat2 -p 6 -x ../genome/genome.fa -1 ../reads/s1_y_1.fq.gz -2 ../reads/s1_y_2.fq.gz -S ./s1.sam

-x 参考基因组 .fa文件
-1 双端测序第1段 .fq.gz格式
-2 双端测序第2段 .fq.gz格式
-S 输出文件地址+名字，输出结果为SAM格式
双端测序一般为read1.fq / read1.fq.gz / read1.fq.bz2格式，前后两端同时出现，并同时作为hisat2输入。将测序得到的fastq文件经过hisat2比对到参考基因组得到SAM文件。SAM文件就是序列比对文件，有上下两个部分，分别包括头部注释部分和比对结果部分：

头部注释部分

//头部注释部分以@开头：
@HD     VN:1.0  SO:unsorted //HD行：VN的版本以及比对排序类型。此处显示SO:unsorted表示没有排列顺序。
@SQ     SN:NC_001133.9  LN:230218//SQ行：参考序列目录。SN：参考序列名字。LN：参考序列长度
@SQ     SN:NC_001148.4  LN:948066
@SQ     SN:NC_001224.1  LN:85779
@SQ     SN:NC_001140.6  LN:562643
@SQ     SN:NC_001141.2  LN:439888
@SQ     SN:NC_001142.9  LN:745751
@SQ     SN:NC_001143.9  LN:666816
@PG     ID:hisat2       PN:hisat2       VN:2.0.5        CL:"/mnt/bai/public/bin/hisat2-2.0.5/hisat2-align-s --wrapper basic-0 -p 8 -x ../genome/genome.fa --dta -S ./s1.sam -1 /tmp/1909596.inpipe1 -2 /tmp/1909596.inpipe2"
//PG行：使用的比对程序名，此处为hisat2

比对结果部分

比对结果部分每行表示一个read和参考基因组的比对结果信息。前11列为主列，包含了大部分重要信息。

SRR5511068.3    99      NC_001147.6     1002516 60      75M     =       1002624 184     GTACTTAACATTCTTCTAATCATGTTAAAAGGTAAAACCTGGCCCATTTTACGATCGATCTGTAAAATCTTATAC     AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE/EEEEAEEEEEEEEEEEAEEEEEEEAEEEEEEEEEEA     AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:75 YS:i:0  YT:Z:CP NH:i:1
SRR5511068.3    147     NC_001147.6     1002624 60      76M     =       1002516 -184    GATAAGTAAGCAATGGTGGTAATTGCAATATTTTGCATATGTGCACGAGAAGAACTATTTTGAAGTAAGATCACTG    /EEEE<EAEEEEEE6E/EEEEEEEAEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEAAAAA    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:76 YS:i:0  YT:Z:CP NH:i:1
SRR5511068.9    83      NC_001146.8     299462  60      75M1S   =       299435  -103    ATTGGAAAAGAAAGTCGCGGCAAAGAGAAATGCCAATAAGACCGGGAATCAAAATTCTAAAAAGAAGAGTCAGAAG    <EAAA6<A<E/EEA/A</AE/EEEAA<EEAE6AEEAAEEAEEEAAAEE/EEEEEEEEEEEE//EEEEEE/EAAAAA    AS:i:-1 XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:75 YS:i:0  YT:Z:CP NH:i:1

主体部分以tab分隔从左到右依次为：

QNAME：Query Name，测序reads名称。如果是双端测序则会出现二次，两端各比对上一次。
FLAG：数值为2的次方数或者其加和，每个2的次方代表一种情况。详情可以参考CSDN大佬的博客：https://blog.csdn.net/genome_denovo/article/details/78712972
或者简书大佬：
https://www.jianshu.com/p/ab133ee9712c
RENAME：Reference Name，双端测序R1比对上的参考序列名，没比上就是*。
POS：position，双端测序R1起始比对上的位置序号。
MAPQ：Mapping Quality，质量分数，越高代表越准确。
CIGAR：比对结果，M代表完全匹配。
RNEXT：双端测序R2端比对情况，比对不上用*号，比对到同一段用=号。
MPOS：双端测序R2端比对位置。
ISIZE：文库插入长度，R2端位置-R1端位置+CIGAR处的值。
SEQ：序列信息。
QUAL：reads质量，用ASCII码表示。

最后总结一下：

转录组数据分析步骤不仅仅是比对产生差异表达基因，后继还会涉及差异表达基因的显著性分析，相关性分析，GO和KEGG分析等等。本文只是总结了转录组分析的上游步骤。而就上游步骤而言，可用的比对软件种类也非常多，本文只是借助hisat2+samtools+stringtie的经典步骤来学习转录组上游分析。其他比对软件日后如有涉及再行整理。