对比软件有很多，首先你可以收藏一下，因为我们是给你介绍的，请统一使用，但要了解它的用法。

直接去主页下载index文件，然后对比fastq格式的reads，得到sam文件。

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

顺便对bam文件做一个简单的QC，参考live 系列。

接下来的四篇文章基本都是规划工作，从这篇文章开始我们才进入了RNA-Seq数据分析的核心部分。

比较还是不比较

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

RNA-Seq 数据分析有多种类型，例如寻找差异表达基因或寻找新的替代剪接。 如果您正在寻找不同表达的基因，您只需要确定不同的读取技术即可。 我们可以使用bwa这样的比较工具，或者像这样的免对齐工具，但前者更快。

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

工具选择

2016年-的综述中提到，目前RNA数据分析有三种策略。 当时的工具也主要使用，STAR和. 其中，其作者已推荐使用HISAT来代替。

最近发表了一篇题为-RNA-seq的文章——被誉为史上最完整的RNA-Seq数据分析流程，也是我一直想做的事情，但他们所做的却超出了我的想象。 文章中根据参考基因组分析转录本所用的工具是、和STAR，推论是找到正确率最高的，总量小于和STAR。 从这里可以看出，第二类错误（接受错误）相对较小，而第一类错误（放弃正确）则较高。

就唯一的比较而言，STAR是两者中最好的，主要是因为它不会像STAR和STAR那样在PE匹配失败时强行将SE与基因组进行比较。 而且在处理较长reads和较短reads的情况下，STAR的稳定性也是最好的。

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

如果你是学习RNA-Seq数据分析，里​​面提到的两篇文档一定要读3遍以上，建议每隔一段时间复习一下。 而且就比较工具而言，基本上就选有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-Seq 下降。 比对结果将存储在/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

& 会将任务丢到后台，所以这3个比较程序会同时执行，如果CPU和显存无法承受，就将&一一去掉。 比较步骤消耗了大量的显存资源，这是由于比较工具将索引数据加载到显存中造成的。 我有64G显存，理论上可以同时处理20个PE数据。 在我的笔记本电脑配置下，同时完成此步骤大约需要2个小时。

基本参数说明

对比数据时，可以悄悄读取额外的选项，主要分为以下几个部分

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

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

输出结果

比较后dnastar序列比对，将输出以下结果。 分析显示，所有数据都是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序列。

同时我也反思了错误的原因。 我默认这个实验重复 3 次 KO 和非 KO。 虽然文章的实验设计不是这样的，但是可以看出理解实验设计是非常重要的。 。

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格式的工具主要是fú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

，orby 在使用时读取名称

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

另外说说默认按照染色体位置排序，-n参数是按照read name排序。 其实还有一个-t可以按TAG排序。

谈论观点

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

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

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

例如在搭配标志（0.1.19版本没有）中，使用-f或-F参数来提取不同匹配情况的读取。

flag是描述reads比较的标志。 有 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

也就是说，如果我想过滤刚刚配对的read，我需要使用0x10

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

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

比较质量控制 (QC)

还是在-上，提到人类基因组的比对率应该是70%到90%，并且多读（multi-reads）的次数应该少一些。 此外，对齐的外显子和对齐的链（读取一个外显子和）的覆盖范围应该一致。

常用的工具有

我们来使用RSeQC，虽然使用的是写入工具，但是自然友好dnastar序列比对，安装也很方便。

# 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 文件，然后查看。

大家好奇为什么一开始是五吗，下面会贴出答案，原文是我的地址。

()

如有侵权请联系删除！