有很多比较软件。 首先，请收集它们。 因为我们是来向您介绍它们的，所以请统一使用它们并了解如何使用它们。

只要去首页下载index文件，然后对比fastq格式的reads就可以得到sam文件了。

然后转成bam文件并排序（注意N排序和P排序的区别）。 索引完成后，加载IGV，截取几个基因看看！

顺便对bam文件做一下简单的QC，参考我的基因组系列直播。

前四篇文章基本上都是准备工作。 从本文开始，我们进入RNA-Seq数据分析的核心部分。

比较还是不比较

在比较之前，我们先要明白比较的目的是什么？ RNA-Seq数据比较和DNA-Seq数据比较有什么区别？

RNA-Seq数据分析可以分为多种类型，例如寻找差异表达基因或寻找新的选择性剪接。 如果我们要寻找差异表达的基因，我们只需要确定不同的read技术即可。 我们可以使用bwa等比较工具，或者align-free工具，后者速度更快。

但如果您需要找到新的或替代的 RNA 剪接，您将需要像 STAR 这样的工具来找到剪接位点。 由于RNA-Seq与DNA-Seq不同，当DNA转录为mRNA时，内含子部分被去除。 因此，如果反向的mRNA cDNA无法与参考序列进行比较，则会将其分离并重新比对，以确定中间是否有内含子。

工具选择

在 2016 年的一篇评论 A of best for RNA-seq data 中，提到目前 RNA 数据分析有三种策略。 当时主要使用的工具是STAR和. 其中，其作者已推荐使用HISAT来代替。

近日，一篇题为《by a Broad-RNA-seq》的文章发表——被称为史上最全面的RNA-Seq数据分析流程。 这也是我一直想做的，但他们做得很好。 超出我的想象。 文章中基于参考基因组进行转录本分析所使用的工具是、和STAR。 结论是准确率最高，但总数少于STAR。 从这里可以看出，第2类错误（接受谎言）相对较少，但第1类错误（抛弃真相）数量较多。

就独特的比较而言，STAR是三者中最好的，主要是因为它在PE比较失败时不会将SE强制到基因组上，就像和一样。 而且，在处理长读和短读的不同情况时，STAR的稳定性也是最好的。

速度方面，比STAR和平均速度快2.5至100倍。

如果你是学习RNA-Seq数据分析，上面提到的两篇文档一定要阅读三遍以上，建议每隔一段时间复习一下。 但就比较工具而言，基本上就选有STAR的。

下载索引

首先，问自己一个问题，为什么比较时需要使用索引？ 强烈推荐大家阅读Jimmy写的算法原理探索和算法原理探索。 但这只是一个建议，你不需要实际看到，反正你也不会理解。

高通量测序遇到的第一个问题是，必须在合理的时间内将数千甚至数亿的reads映射到参考基因组上，并且错误率必须在可接受的范围内。 为了提高比对速度，需要通过BWT算法将参考基因组序列转换为索引，而我们比对的序列实际上是索引的子集。 当然，转录组比较还需要考虑选择性剪接，因此比较复杂。

因此，我们不是将读数直接粘贴到基因组上，而是将读数与索引进行比较。 人类索引一般都是现成的。 我建议您下载现成的。 我曾经尝试使用服务器自己创建索引，所花费的时间让我怀疑自己的人生。

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

如果你觉得你的电脑配置还可以，或者没有现成的索引，可以通过方法创建

# 其实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

我的是i7-7700处理器，内存64G。 运行资源效率如下：

形式比较

基本用法是 []* -x {-1 -2 | -U } [-S]，基本上提供了索引、PE数据或SE数据存储位置的位置。 然而，其他可选参数是一个很好的进步方式。 如果您是新手，则使用默认参数即可。

因为RNA-Seq数据来自~，6个样本的PE数据，也就是有6个循环，可以写脚本或者直接在命令行运行

下面命令运行的目录是/mnt/f/Data/，我的参考基因组的索引数据存放在/mnt/f/Data//index/hg19/，RNA-seq数据存放在/ mnt/f/Data/ 在 RNA 测序下。 比较结果将存储在/mnt/f/Data/RNA-Seq/中

mkdir -p RNA-Seq/aligned
for i in `seq 57 62`
do
    hisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 SRR35899${i}_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam &
done

& 会将任务扔到后台，因此这三个比较程序会同时执行。 如果CPU和内存无法承受，请将&一一删除。 比较步骤消耗大量内存资源。 这是由于比较工具必须将索引数据放入内存而导致的。 我有64G内存，理论上可以同时处理20个PE数据。 在我的电脑配置下dnastar序列比对，同时完成这一步大约需要2个小时。

基本参数说明

在数据比较过程中，您可以有额外的安静读取选项，主要分为以下几个部分：

注：soft表示比对的reads中只有部分与参考序列匹配，有部分不匹配。 也就是说，100bp的read会匹配前20bp或者后20bp，或者最后20bp的比较效果不太好。

因此，比较结果受比较选项、评分选项、可变剪切选项和PE选项的影响。 在我有生之年，我应该写一篇文章来介绍这些选项对结果的影响。

输出结果

比较后，将输出以下结果。 解读一下，所有数据都是100%，96.68%的配对数据不是一次性比较的，1.23%的数据比较是唯一比较，2.09%是多重比较。 那么96.68%的人根本没有数据可以比较。 如果不按顺序，就会有0.05%的比较。 如果再用单端数据进行比较，则95.20%的数据没有进行比较，3.60%的数据进行了一次比较，1.20%的数据进行了多次比较。 零和零的总和加起来就是8%的对比！ ！ ！

这个整体对比率让我对生活产生了质疑。 怎么可能？ 我查看了输出记录，发现有几个结果的比对率都在90%以上。 我想了想，发现这个实验似乎是在人类和老鼠身上都做过的。 于是又去GEO查了一下记录，突然发现自己差点翻车了。

图9-15是人293细胞(9-11)或小鼠ES细胞(12-15)中的mRNA序列。

同时我也反思了错误的原因。 我默认 KO 和非 KO 各重复 3 次。 事实上，文章的实验设计并不是这样的。 可见，理解实验设计非常重要。 。

mkdir -p RNA-Seq/aligned
for i in `seq 56 58`
do
    hisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 SRR35899${i}_2.fastq.gz -S RNA-Seq/SRR35899${i}.sam &
done

以上是修改后的代码

三轴

SAM（/）数据格式是目前高通量测序中存储比对数据的标准格式。 当然也可以用来存储未对齐的数据。因此，SAM的格式说明

目前处理SAM格式的工具主要是Heng Li编写的。 除了C语言版本外，还有Java、Pysam、lisp的cl-sam等版本。 主要功能如下：

最常用的三种工具是格式转换、排序和索引。 高级教程是通过阅读文档来提高的。

for i in `seq 56 58`
do
    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
done

笔记

Jimmy说要仔细判断sam的两种排序方式的区别，于是我截取了前100行数据，分别排序，然后查看结果。

head -1000 SRR3589957.sam > test.sam
samtools view -b  test.sam > test.bam
samtools view test.bam | head

默认排序

samtools sort test.bam default
samtools view default.bam | head

按 排序，或使用 -n 时按读取名称排序

samtools sort -n test.bam sort_left
samtools view sort_left.bam | head

也就是说，默认是根据染色体位置排序dnastar序列比对，而-n参数则是根据read name排序。 当然还有一个-t根据TAG排序。

谈谈看法

的view是一个非常实用的子命令。 除了前面的格式转换之外，它还可以进行数据提取和提取。

例如，提取1号染色体1234区域的比对读数

samtools view SRR3589957_sorted.bam chr1:1234-123456 | head

例如，使用flag（0.1.19版本中不可用）和 ，使用-f或-F参数来提取不同匹配情况的reads。

Flag是一种描述读取比较情况的标记。 共有12种，可以组合使用。

0x1    PAIRED    paired-end (or multiple-segment) sequencing technology
0x2    PROPER_PAIR    each segment properly aligned according to the aligner
0x4    UNMAP    segment unmapped
0x8    MUNMAP    next segment in the template unmapped
0x10    REVERSE    SEQ is reverse complemented
0x20    MREVERSE    SEQ of the next segment in the template is reverse complemented
0x40    READ1    the first segment in the template
0x80    READ2    the last segment in the template
0x100    SECONDARY    secondary alignment
0x200    QCFAIL    not passing quality controls
0x400    DUP    PCR or optical duplicate
0x800    SUPPLEMENTARY    supplementary alignment

你可以先查看一下整体情况

samtools flagstat SRR3589957_sorted.bam

换句话说，如果我想过滤完全匹配的reads，我需要使用0x10

samtools view -b -f 0x10 SRR3589957_sorted.bam chr1:1234-123456  > flag.bam
samtools flagstat flag.bam

我有生之年应该写一篇文章好好介绍一下。

比较质量控制 (QC)

还是在A of best for RNA-seq data中，提到人类基因组的比对率应该是70%到90%，并且multi-reads的次数要少。 此外，外显子上的读取和链（外显子上的读取和 ）的覆盖范围应该一致。

常用的工具包括

让我们使用 RSeQC。 毕竟使用自己写的工具本来就很友好，安装也很方便。

# Python2.7环境下
pip install RSeQC

总共有以下几个文件。 根据它们的名字就可以知道它们的功能是什么。

首先对bam文件进行统计分析。 从结果来看，满足70~90的对比率要求。

bam_stat.py -i SRR3589956_sorted.bam

基因组覆盖度的QC需要提供bed文件，可以直接从RSeQC网站下载，也可以使用gtf转换

read_distribution.py -i RNA-Seq/aligned/SRR3589956_sorted.bam -r reference/hg19_RefSeq.bed

IGV视图

加载参考序列、注释和 BAM 文件并查看。

你好奇为什么一开始是五吗？ 下面揭晓答案。 原文是我的地址。

()

如有侵权请联系删除！