RNAseq005 转录组入门(5)：序列比对

作者: caoqiansheng | 来源:发表于2020-09-02 00:13 被阅读0次

1 比对的目的

1.1 差异表达基因
只需要确定不同的read计数，可以用bowtie, bwa这类比对工具，或者是salmon这类align-free工具，并且后者的速度更快。
1.2 寻找新剪切位点
但是如果需要找新的isoform，或RNA的可变剪切，看外显子使用差异的话，需要TopHat, HISAT2或者是STAR这类工具用于找到剪切位点。

2 RNA-Seq数据比对和DNA-Seq数据比对有什么差异

RNA-Seq不同于DNA-Seq，DNA在转录成mRNA的时候会把内含子部分去掉，所以mRNA反转的cDNA如果比对不到参考序列，会被分开，重新比对一次，判断中间是否有内含子。

3 工具抉择

在2016年的一篇综述A survey of best practices for RNA-seq data analysis，提到目前有三种RNA数据分析的策略。那个时候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已经被它的作者推荐改用HISAT进行替代。
最近的Nature Communication发表了一篇题为的Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis的文章—被称之为史上最全RNA-Seq数据分析流程，文章中在基于参考基因组的转录本分析中所用的工具，是TopHat,HISAT2和STAR，结论就是HISAT2找到junction正确率最高，但是在总数上却比TopHat和STAR少。从这里可以看出HISAT2的二类错误(纳伪）比较少，但是一类错误（弃真）就高起来。
就唯一比对而言，STAR是三者最佳的，主要是因为它不会像TopHat和HISAT2一样在PE比对不上的情况还强行把SE也比对到基因组上。而且在处理较长的read和较短read的不同情况，STAR的稳定性也是最佳的。
就速度而言，HISAT2比STAR和TopHat2平均快上2.5~100倍。

3 index下载

3.1 hisat2主页下载

http://daehwankimlab.github.io/hisat2/

image.png

http://daehwankimlab.github.io/hisat2/download/

image.png

3.2 http://ccb.jhu.edu/software.shtml

http://ccb.jhu.edu/software/hisat2/manual.shtml
ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/

image.png

cd referece && mkdir index && cd index
# hg19 index
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
# hg38 index
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg38.tar.gz

将上述链接复制放入迅雷，十分钟搞定下载

4.比对alignment

4.1 说明书：https://daehwankimlab.github.io/hisat2/manual/

4.2 用法

hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | –sra-acc <SRA accession number>} [-S <hit>]

主要参数：
-x <hisat2-idx>: 参考基因组索引文件的前缀
如下图，hisat2 hg19index解压后如下，-x参数只需要index的前缀，也就是genome，所以直接在路径后加入genome即可，如果直接使用路径，或者使用正则表达式genome.*,或是cat合并成一个genome都是不可取的
hisat2 hg19index

-1 <m1>: 双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
-2 <m2>: 双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
-U <r>: 单端数据文件。若有多组数据，使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。
–sra-acc <SRA accession number>: 输入SRA登录号，比如SRR353653，SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型，进行比对。
-S <hit> : 指定输出的SAM文件。

输入选项：
-q:输入文件为FASTQ格式。FASTQ格式为默认参数。
-qseq:输入文件为QSEQ格式。
-f:输入文件为FASTA格式。
-r:输入文件中，每一行代表一条序列，没有序列名和测序质量等。选择此项时，–ignore-quals参数也会被选择。
-c:此参数后是直接比对的序列，而不是包含序列的文件名。序列间用逗号隔开。选择此项时，–ignore-quals参数也会被选择。
-s/–skip <int>:跳过输入文件中前条序列进行比对。
-u/–qupto <int>:只使用输入文件中前条序列进行比对，默认是没有限制。
-5/–trim5 <int>:比对前去除每条序列5’端个碱基
-3/–trim3 <int>:比对前去除每条序列3’端个碱基
–phred33:输入的FASTQ文件碱基质量值编码标准为phred33，phred33为默认参数。

文章中提到：Samples 9-15 are mRNA-seq to determine effect of AKAP95 knockdown in human 293 cells (9-11) or mouse ES cells (12-15).

# -t参数可以显示时间
# -x参数指定index位置：/mnt/e/bioinfo/reference/gtf/hisat2_index/hg19/，需要在路径后添加index的前缀genome，否则会报错
# -1，-2（双向测序）参数指定数据存放位置：/mnt/e/bioinfo/rna_seq/ENA_data/
# -S指定输出sam的路径，仅有三个是人基因组数据
for i in {56..58};do hisat2 -t -x /mnt/e/bioinfo/reference/gtf/hisat2_index/hg19/genome  -1 /mnt/e/bioinfo/rna_seq/ENA_data/SRR35899${i}_1.fastq.gz  -2  /mnt/e/bioinfo/rna_seq/ENA_data/SRR35899${i}_2.fastq.gz   -S  /mnt/e/bioinfo/rna_seq/rna_seq/aligned/SRR35899$i.sam;done

我的笔记本是8G内存，差不多是半小时跑完一个数据，以下为比对运行结果

image.png

Hisat2 bowtie2比对结果解读（Hisat2 Alignment summary）

第一部分描述的是pair-end模式下的一致比对结果
1820962 (6.31%) aligned concordantly 0 times 是什么意思？aligned concordantly就是read1和read2同时合理的比对到了基因组/转录组上。
24723672 (85.68%) aligned concordantly exactly 1 time，exactly 1 time 就是只有一种比对结果。>1 times 就是read1和read2可以同时比对到多个地方。
第二部分，pair-end模式下不一致的比对结果
1820962 pairs aligned concordantly 0 times; of these: 90371 (4.96%) aligned discordantly 1 time，concordantly 0 times就是read1和read2不能同时合理的同时比对到基因组/转录组上
aligned discordantly 1 time 最难理解，discordantly比对就是read1和read2都能比对上，但是不合理，
1.比对方向不对，pair-end测序的方向是固定的；
2.read1和read2的插入片段长度是有限的
第三部分就是对剩余reads（既不能concordantly，也不能discordantly 1 time）的单端模式的比对，没什么好说的。

5 SAMtools

Samtools是一个用于操作sam和bam格式文件的应用程序集合，具有众多的功能。它从SAM（序列比对/映射）格式导入和导出，进行排序，合并和索引，并允许快速检索任何区域中的读数。SAM和BAM格式的比对文件主要由bwa、bowtie2、tophat和hisat2等序列比对工具产生，用于记录测序reads在参考基因组上的映射信息。其中，BAM格式文件是SAM文件的的二进制格式，占据内存较小且运算速度更快。

5.1 语法

Program: samtools (Tools for alignments in the SAM format)
Version: 1.10 (using htslib 1.10)
Usage:   samtools <command> [options]
Commands:
  -- Indexing 索引
     dict           创建一个序列字典文件
     faidx          建立Fasta索引文件，以.fai后缀结尾
     fqidx          建立Fastq索引文件
     index          建立sam索引文件，需对bam文件进行排序后才能构建索引
  -- Editing 编辑
     calmd          重新计算MD/NM标签和'='基数
     fixmate        为以名称排序定位的alignment填入配对坐标
     reheader       替换BAM头注释
     targetcut      切割fosmid区域（仅适用于fosmid池）
     addreplacerg   添加或替换RG标签
     markdup        重复标记
  -- File operations 文件操作
     collate        根据read name信息将bam文件进行随机打散和分组
     cat            将多个bam文件合并为一个bam文件，无需对sam文件进行排序
     merge          将多个已经sort的bam文件进行合并
     mpileup        对参考基因组每个位点做碱基堆积，用于call SNP和INDEL。主要是生成BCF、VCF文件或者pileup一个或多个bam文件
     sort           对比对后的bam文件进行排序，默认按碱基位置坐标进行排序
     split          按读取组分割文件
     quickcheck     快速检查SAM/BAM/CRAM是否完整
     fastq          将BAM转换为FASTQ
     fasta          将BAM转换为FASTA
  -- Statistics 统计
     bedcov         统计给定region区间内reads的深度，每个碱基的深度叠加在一起
     coverage       统计每个碱基位点的测序深度及百分比
     depth          统计每个碱基位点的测序深度，注意计算前必须对bam文件排序并构建索引
     flagstat       统计bam文件中read的比对情况
     idxstats       统计bam索引文件里的比对信息
     phase          杂合子项
     stats          对bam文件做详细统计, 统计结果可用mics/plot-bamstats作图
  -- Viewing 查看
     flags          查看不同flag值的含义
     tview          直观地显示reads比对到参考基因组的情况，类似于基因组浏览器
     view           SAM<->BAM<->CRAM转换
     depad          将填充的BAM转换为未填充的BAM

5.2 Samtools常用命令

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

转换：samtools view -S SRR3589912.sam -b > SRR3589912.bam（-S是最新版的samtools为了兼容以前的版本写的）
排序：samtools sort SRR3589912.bam -o SRR3589912_sorted.bam
索引：samtools index SRR3589912_sorted.bam
比对后的sam文件应该先进行格式转换，接着是排序，最后根据排序的文件建立索引文件

5.2.1 格式转换 samtools view

5.2.1.1 用法

Usage: samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]

Options:
  -b             输出BAM
  -C             输出CRAM (requires -T)
  -1             使用快速BAM压缩(implies -b)
  -u             未压缩的BAM输出(implies -b)
  -h             在SAM输出中包括头注释
  -H             仅打印SAM标头（无比对信息）
  -c             仅打印匹配记录的数量
  -o FILE        输出文件名[stdout]
  -U FILE        输出未过滤的reads到文件 [null]
  -t FILE        列出参考序列的名称和长度[null]
  -X             包含通用的索引文件
  -L FILE        仅包括与该BED FILE重叠的Reads [null]
  -r STR         仅包含在STR组中的Reads [null]
  -R FILE        仅包含FILE [null]中列出的读取组的Reads
  -d STR:STR     仅包含带有标签STR和关联值STR的Reads [null]
  -D STR:FILE    仅包含带有标签STR以及FILE [null]中列出的相关值的Reads
  -q INT         仅包含映射质量> = INT [0]的Reads
  -l STR         仅包括对库STR的Reads[null]
  -m INT         仅包含消耗CIGAR操作数的Reads
  -f INT         仅包括存在INT中所有FLAG的Read[0]
  -F INT         仅包括在INT中不存在任何FLAGS的Read[0] [0]
  -G INT         仅排除当前FLAGs为INT的Read [0]
  -s FLOAT       子样本读取（给定INT.FRAC选项值，0。FRAC是要保留的模板/读取对的分数； INT部分设置种子）
  -M             使用多区域迭代器（提高速度，删除重复项并按文件中的顺序输出读取结果）
  -x STR         读取标记以剥离[null]
  -B             折叠CIGAR操作
  -?             帮助，包括有关区域规范的注释
  -S             忽略（自动检测输入格式）
  --no-PG  do not add a PG line
      --input-fmt-option OPT[=VAL]                      以OPTION或OPTION = VALUE的形式指定一个输入文件格式选项
  -O, --output-fmt FORMAT[,OPT[=VAL]]...                指定输出格式（SAM，BAM，CRAM）
      --output-fmt-option OPT[=VAL]                     以OPTION或OPTION = VALUE的形式指定单个输出文件格式
  -T, --reference FILE                参考序列FASTA FILE 
  -@, --threads INT                   要使用的其他线程数[0]
      --write-index                   自动索引输出文件[关闭]
      --verbosity INT                 设置详细程度

5.2.1.2 samtools的view不仅可以进行格式转换，还可以进行数据的提取

# 提取1号染色体上1234~123456区域的以对read
samtools view SRR3589957_sorted.bam chr1:1234-123456 | head

5.2.1.3 SAM文件中FLAG值的理解

FLAG列在SAM文件的第二列，这是一个很重要的列，包含了很多mapping过程中的有用信息

image.png

SAM文件的flag标签包含0x4

### 小写的f是提取，大写的F是过滤
samtools view -f4 sample.bam sample.unmapped.sam

对于PE数据，在未比对成功的测序数据可以分成3类：

仅reads1没有比对成功，该提取条件包括：
该read是read1，对应于二进制FLAG的第7位，该位取1，其十进制值为64；
该read未成功比对到参考基因组，对应于二进制FLAG的第3位，该位取1，其十进制值为4；
另一配对read成功比对到参考基因组，对应于二进制FLAG的第4位，该位取0，其十进制值为8；

 # 对于取1的位点采用提取的策略，用-f参数，值设为64+4=68
 # 对于取0的位点采取过滤的策略，用-F参数，值设为8
 samtools view -u -f 68 -F 8 alignments.bam >read1_unmap.bam

仅reads2没有比对成功，该提取条件包括：
该read是read2，对应于二进制FLAG的第8位，该位取1，其十进制值为128；
该read未成功比对到参考基因组，对应于二进制FLAG的第3位，该位取1，其十进制值为4；
另一配对read成功比对到参考基因组，对应于二进制FLAG的第4位，该位取0，其十进制值为8；

# 对于取1的位点采用提取的策略，用-f参数，值设为128+4=132
# 对于取0的位点采取过滤的策略，用-F参数，值设为8
samtools view -u -f 132 -F 8 alignments.bam >read2_unmap.bam

两端reads都没有比对成功，该提取条件包括：
该read未成功比对到参考基因组，对应于二进制FLAG的第3位，该位取1，其十进制值为4；
另一配对read未成功比对到参考基因组，对应于二进制FLAG的第4位，该位取1，其十进制值为8；

# 对于取1的位点采用提取的策略，用-f参数，值设为4+8=12
$ samtools view -u -f 12 alignments.bam >pairs_unmap.bam

5.2.2 排序 samtools sort

Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     设置输出文件压缩等级。0-9，0是不压缩，9是压缩等级最高。不设置此参数时，使用默认压缩等级
  -m INT     设置每个线程运行时的内存大小，可以使用K，M和G表示内存大小
  -n         设置按照read名称进行排序
  -t TAG     Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
  -o FILE    设置最终排序后的输出文件名
  -T PREFIX  设置临时文件的前缀 PREFIX.nnnn.bam
  --no-PG不添加PG行
      --input-fmt-option OPT [= VAL]                      以OPTION或OPTION = VALUE的形式指定一个输入文件格式选项
  -O，--output-fmt FORMAT [，OPT [= VAL]] ...             指定输出格式（SAM，BAM，CRAM）
      --output-fmt-option OPT [= VAL]                     以OPTION或OPTION = VALUE的形式指定单个输出文件格式选项
      --reference FILE                                    参考序列FASTA FILE [null]
  -@, --threads INT               设置排序和压缩是的线程数量，默认是单线程
      --verbosity INT               Set level of verbosity

对bam文件进行排序，不能对sam文件进行排序。以leftmost coordinates的方式对比对结果进行排序，或者使用-n参数以read名称进行排序。将会添加适当的@HD-SO排序顺序标头标签或者如果有必要的话，将会更新现存的一个排序顺序标头标签。sort命令的输出默认是标准输出写入，或者使用-o参数时，指定bam文件输出名。sort命令还会在内存不足时创建临时文件tmpprefix.%d.bam。

samtools的排序方式有两种（常用）

# 默认方式，按照染色体的位置进行排序
 for i in {56..58};do samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam;done
# 参数-n则是根据read名进行排序。
 for i in {56..58};do samtools sort -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam;done

image.png

bam文件是Sam 文件的二进制压缩格式，保留了与sam 完成相同的内容信息。SAM/BAM 文件可以是未排序的，但是按照坐标(coodinate)排序可以线性的监控数据处理过程。samtools可以用来转化bam/sam文件，可以merg，sort aligment,可以去除duplicate,可以call snp及indels.
BAM 文件是压缩的二进制文件，对文件内容排序之后相似的内容排在一起，使得文件压缩比提高了，因此排序之后的 BAM 文件变小了，相对应的 SAM 文件就是纯文本文件，对 SAM 文件进行排序就不会改变文件大小。而且由于 RNA-seq 中由于基因表达量的关系，RNA-seq 的数据比对结果 BAM 文件使用 samtools 进行 sort 之后文件压缩比例变化会比 DNA-seq 更甚。

5.2.3 索引 samtools index

Usage: samtools index [-bc] [-m INT] <in.bam> [out.index]
Options:
  -b       创建bai索引文件，未指定输出格式时，此参数为默认参数
  -c       创建csi索引文件，默认情况下，索引的最小间隔值为2^14，与bai格式一致
  -m INT   创建csi索引文件，最小间隔值2^INT
  -@ INT   设置线程数[无]

为了能够快速访问bam文件，可以为已经基于坐标排序后bam或者cram的文件创建索引，生成以.bai或者.crai为后缀的索引文件。,必须使用排序后的文件。另外，不能对sam文件使用此命令。如果想对sam文件建立索引，那么可以使用tabix命令创建。在使用与region参数相关的其它samtools命令时，需要创建索引

# -i 参数直接用 {}就能代替管道之前的标准输出的内容
ls *sorted.bam | xargs -i samtools index {}
# 或是ls *.bam | xargs -n1  samtools index

xargs命令可参考：https://www.jianshu.com/writer#/notebooks/45449543/notes/76201605/preview

image.png
报错：默认排序可以建立index，但是按reads名排序（sort -n）无法建立index

Reference

https://www.jianshu.com/p/681e02e7f9af
https://www.jianshu.com/p/68f6e35fa4a2
https://ming-lian.github.io/2019/02/07/Advanced-knowledge-of-SAM/

RNAseq005 转录组入门(5)：序列比对

1 比对的目的

2 RNA-Seq数据比对和DNA-Seq数据比对有什么差异

3 工具抉择

3 index下载

3.1 hisat2主页下载

3.2 http://ccb.jhu.edu/software.shtml

4.比对alignment

4.1 说明书：https://daehwankimlab.github.io/hisat2/manual/

4.2 用法

5 SAMtools

5.1 语法

5.2 Samtools常用命令

5.2.1 格式转换 samtools view

5.2.1.1 用法

5.2.1.2 samtools的view不仅可以进行格式转换，还可以进行数据的提取

5.2.1.3 SAM文件中FLAG值的理解

5.2.2 排序 samtools sort

5.2.3 索引 samtools index

Reference

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

NGS测序笔记