搜索
查看: 5457|回复: 4

[mRNA-seq] 转录组入门分析5_序列比对

[复制链接]

9

主题

15

帖子

132

积分

注册会员

Rank: 2

积分
132
发表于 2017-8-27 04:41:12 | 显示全部楼层 |阅读模式
本帖最后由 laofuzi 于 2017-8-27 04:43 编辑

1.     目的:
这次Jimmy群主布置的作业比较多,我个人将任务分解了一下,大致分为以下几个小任务:
1.1了解RNA-seq比对软件
1.2 学会下载hisat2的index文件,把fastq格式的reads比对得到sam文件
1.3 学会将sam文件转为bam文件
1.4 排序,注意N和P两种排序区别
1.5 将索引好bam文件载入IGV,截图几个基因
1.6 便对bam文件进行简单QC
2.     RNA-seq比对软件
首先说一下HISAT2,HISAT全称为Hierarchical Indexing for Spliced Alignment of Transcripts,由约翰霍普金斯大学开发。它取代Bowtie/TopHat程序,能够将RNA-Seq的读取与基因组进行快速比对。HISAT利用大量FM索引,以覆盖整个基因组。Index的目的主要使用与序列比对。由于物种的基因组序列比较长, 如果将测序序列与整个基因组进行比对,则会非常耗时。因此采用将测序序列和参考基因组的Index文件进行比对,会节省很多时间。
以人类基因组为例,它需要48,000个索引,每个索引代表~64,000 bp的基因组区域。这些小的索引结合几种比对策略,实现了RNA-Seq读取的高效比对,特别是那些跨越多个外显子的读取。尽管它利用大量索引,但HISAT只需要4.3 GB的内存。这种应用程序支持任何规模的基因组,包括那些超过40亿个碱基的。(来自plob)
另外一个软件是STAR(SplicedTranscripts Alignments to a Reference,STAR)软件,是ENCODE计划御用软件之一,可以装配短序列和长序列。
在使用上,HISAT2表现出最快的速度和最准确的拼接比对,但是没有STAR的敏感度高。HISAT2找到junction正确率最高,但是在总数上却比TopHat和STAR少。HISAT2的二类错误(纳伪)比较少,但是一类错误(弃真)就高起来。 就唯一比对而言,STAR是三者最佳的,主要是因为它不会像TopHat和HISAT2一样在PE比对不上的情况还强行把SE也比对到基因组上。而且在处理较长的read和较短read的不同情况,STAR的稳定性也是最佳的。 就速度而言,HISAT2比STAR和TopHat2平均快上2.5~100倍。(来自徐洲更)
3.     hisat2的Index的下载
Hisat2的官网地址为https://ccb.jhu.edu/software/hisat2/index.shtml。由于文章采用的是hg19、mm10参考基因组。所以,到官网的右侧Index列表中下载。
[AppleScript] 纯文本查看 复制代码
# 构建目录
mkdir /mnt/d/reference/index
mkdir /mnt/d/reference/index/hisat
#下载
wget -P /mnt/d/reference/index/hisat [url=ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz]ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz[/url]
wget -P /mnt/d/reference/index/hisat [url=ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz]ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz[/url]
#解压,-C指定解压目录
tar -zxvf /mnt/d/reference/index/hisat/hg19.tar.gz -C /mnt/d/reference/index/hisat
tar -zxvf /mnt/d/reference/index/hisat/mm10.tar.gz -C /mnt/d/reference/index/hisat
# 移除下载的压缩包
rm /mnt/d/reference/index/hisat/*.tar.gz
#查看解压文件
ll /mnt/d/reference/index/hisat/hg19
ll /mnt/d/reference/index/hisat/mm10

Tips_1:不用单独给下载的hg19和mm10的index文件再分别构建目录,解压的时候会自动创建hg19和mm10目录。

Tips_2: 解压的文件中,包含genome.*.ht2的8个文件和一个shell脚本。

Tips_3:如果所要分析的物种,没有现成的index文件,还可以自行构建index文件。在解压缩的文件里面,提供了make_hg19.sh和make_mm10.sh两个文件,里面详细记录了创建索引的过程,可以参考构建自己的index文件。

这里,微吐槽一下,Hisat2提供的index文件没有给出md5值,下载的完整性如何也不知道。


4.     hisat2命令参数

hisat2基本用法就是hisat2 [options]* -x <hisat2-idx>{-1 <m1> -2 <m2> | -U <r> } [-S <hit>],基本就是提供index的位置,PE数据或者是SE数据存放位置。最主要的四个参数是:

-x 指定基因组索引

-1 指定第一个fastq文件

-2 指定第二个fastq文件

-S 指定输出的SAM文件

hisat2的官网操作手册地址为:https://ccb.jhu.edu/software/hisat2/manual.shtml

主要参数:

-x<hisat2-idx>  

Thebasename of the index for the reference genome. The basename is the name of anyof the index files up to but not including the final .1.ht2 / etc. hisat2 looksfor the specified index first in the current directory, then in the directoryspecified in the HISAT2_INDEXES environment variable.

提供参考基因组索引文件index的位置。


-1<m1>

Comma-separatedlist of files containing mate 1s (filename usually includes _1), e.g. -1flyA_1.fq,flyB_1.fq. Sequences specified with this option must correspondfile-for-file and read-for-read with those specified in <m2>. Reads maybe a mix of different lengths. If - is specified, hisat2 will read the mate 1sfrom the "standard in" or "stdin" filehandle.

双端测序结果的第一个文件。若有多组数据,使用逗号将文件分隔。Reads的长度可以不一致。


-2<m2>

Comma-separatedlist of files containing mate 2s (filename usually includes _2),e.g. -2 flyA_2.fq,flyB_2.fq.Sequences specified with this option must correspond file-for-file andread-for-read with those specified in <m1>.Reads may be a mix of different lengths. If -is specified, hisat2 willread the mate 2s from the "standard in" or "stdin"filehandle.

双端测序结果的第二个文件。若有多组数据,使用逗号将文件分隔,并且文件顺序要和-1参数对应。Reads的长度可以不一致。


-U<r>

Comma-separatedlist of files containing unpaired reads to be aligned, e.g. lane1.fq,lane2.fq,lane3.fq,lane4.fq.Reads may be a mix of different lengths. If -is specified, hisat2 getsthe reads from the "standard in" or "stdin" filehandle.

单端数据文件。若有多组数据,使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。


–sra-acc<SRA accession number>

Comma-separatedlist of SRA accession numbers, e.g. --sra-accSRR353653, SRR353654. Information about read typesis available at http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?sp=runinfo&acc=sra-acc&retmode=xml,where sra-acc is SRA accession number. If users run HISAT2 on a computercluster, it is recommended to disable SRA-related caching (see the instructionat SRA-MANUAL).

输入SRA登录号,比如SRR353653,SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型,进行比对。


-S<hit>

Fileto write SAM alignments to. By default, alignments are written to the"standard out" or "stdout" filehandle (i.e. the console).

指定输出的SAM文件。

在emacs中编辑如下脚本:
[AppleScript] 纯文本查看 复制代码
#!/bin/bash

for i in `seq 56 58`
do
hisat2 -t -x /mnt/d/reference/index/hisat/hg19/genome -1 /mnt/d/AKAP95/fastq/SRR35899${i}_1.fastq.gz -2 /mnt/d/AKAP95/fastq/SRR35899${i}_2.fastq.gz -S  /mnt/d/AKAP95/aligned/SRR35899${i}.sam 
done

for i in `seq 59 62`
do
hisat2 -t -x /mnt/d/reference/index/hisat/mm10/genome -1 /mnt/d/AKAP95/fastq/SRR35899${i}_1.fastq.gz -2 /mnt/d/AKAP95/fastq/SRR35899${i}_2.fastq.gz -S  /mnt/d/AKAP95/aligned/SRR35899${i}.sam 
done

将脚本存于aligned目录下,命名为map.sh。

Tips_1: -t参数主要是显示输入加载index文件和比对序列所用的时间。

Tips_2: 在命令中,-x参数的命令格式为,-x <hisat2-idx>,本例中为-x/mnt/d/reference/index/hisat/mm10/genome,其中的genome是index文件的前缀(<ht2-idx>  Index filename prefix (minus trailing.X.ht2).),本人在刚开始键入命令的时候,以为genome是一个文件夹,所以也就没有加genome,也并没有在意,但是运行代码的时候报错:“Could not locatea HISAT2 index corresponding to basename”,后来问了一下谷哥才知道问题所在。

Tips_3:徐洲更在其代码中“hisat2 …SRR35899${i}.sam &”提及,&会把任务丢到后台,所以会同时执行这3个比对程序,如果CPU和内存承受不住,去掉&一个个来。本人电脑有点渣,没试过。
[AppleScript] 纯文本查看 复制代码
#运行该脚本
sh /mnt/d/AKAP95/aligned/map.sh

Tips_1: 比对结果的解读

Time loadingforward index: 00:00:06

Time loadingreference: 00:00:00

Multiseedfull-index search: 00:26:46

28856780 reads;of these:

  28856780 (100.00%) were paired; of these:   

    1838758 (6.37%) aligned concordantly 0times    # 1838758(6.37%)没有一致对齐

    24733251 (85.71%) aligned concordantlyexactly 1 time  #  24733251(85.71%)精确的一次对齐

    2284771 (7.92%) aligned concordantly >1times   #2284771 (7.92%) 有超过一次的定位

    ----

    1838758 pairs aligned concordantly 0 times;of these:   

      90903 (4.94%) aligned discordantly 1 time

    ----

    1747855 pairs aligned 0 times concordantlyor discordantly; of these:

      3495710 mates make up the pairs; ofthese:

        2034758 (58.21%) aligned 0 times

        1221302 (34.94%) aligned exactly 1 time  

        239650 (6.86%) aligned >1 times

96.47% overallalignment rate

Time searching:00:26:46

Overall time: 00:26:52

本人比较了本人的单线程跑的结果(如上)和其他人跑的结果。

#下图是rumbleD

#下图是Panda姐的结果

可以看出,三个结果略有不同,但差距不大。
本人查了一下官网手册,其中参数-p下面有一段话:
If your computer has multipleprocessors/cores, use -p. The -p option causes HISAT2 to launch a specifiednumber of parallel search threads. Each thread runs on a differentprocessor/core and all threads find alignments in parallel, increasingalignment throughput by approximately a multiple of the number of threads (though in practice, speedup is somewhat worse than linear) (加速得到的结果会比单线程略差
这三个结果差异不大,正如官网手册说的“though in practice,speedup is somewhat worse than linear”。

5.     SAM转化为BAM

SAM(sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式。目前处理SAM格式的工具主要是SAMTools。Samtools是一个用来处理BAM格式(SAM的二进制格式)的比对文件的工具箱。它能够输入和输出SAM(序列比对/图谱)格式的文件,对其进行排序、合并、建立索引等处理,还可以快速获取任意区域的reads。 这个网上资料很多,基本上就两个版本。
编辑如下脚本,另存为sam.sh,并运行
[AppleScript] 纯文本查看 复制代码
#!/bin/bash
#转到特定文件夹
cd /mnt/d/AKAP95/aligned
# 首先将比对后的sam文件转换成bam文件
# 利用的是samtools的view选项,参数-S 输入sam文件;参数-b 指定输出的文件为bam;最后重定向写入bam文件
for i in `seq 56 62`
  do samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
done
# 将所有的bam文件进行排序
for i in `seq 56 62`
  do samtools sort SRR35899${i}.bam -o SRR35899${i}.sorted.bam
done
# 将所有的排序文件建立索引,索引文件.bai后缀
for i in `seq 56 62`
  do samtools index SRR35899${i}.sorted.bam
done

sh /mnt/d/AKAP95/aligned/sam.sh

6.    两种排序方法
在先前阅读官方手册的时候,就已经知道对于samtools的sort命令在对序列进行排序的时候,存在两种排序方法。默认通过最左边的坐标排序(染色体的位置);当使用 -n 时,通过read的名字排序。
为了有一个更直观的认识,还是直接做一遍。采用徐洲更的代码。
[AppleScript] 纯文本查看 复制代码
cd /mnt/d/AKAP95/aligned
head -1000 SRR3589956.sam > test.sam
samtools view -b test.sam > test.bam
samtools view test.bam | head

#下图为未排序前
[AppleScript] 纯文本查看 复制代码
samtools sort test.bam -o test.sorted.bam
samtools view test.sorted.bam | head
#下图为采用默认排序,可以看出以染色体位置进行排序
[AppleScript] 纯文本查看 复制代码
head -1000 SRR3589956.sam > test2.sam
samtools view -b test2.sam > test2.bam
samtools view test2.bam | head
samtools sort test2.bam -n -o test2.sorted.bam
samtools view test2.sorted.bam | head

# 当使用 -n 时,通过read的名字排序

7.    序列比对质控
比对完,通过统计比对率等相关信息,判断比对结果可否通过二次质控(QC of Alignment)。如若通过,方可进行后续分析。
常用工具有
采用RSeQC进行QC。查了一下官网,需要python2.7的安装环境(Prerequisite: gcc; python2.7; numpy; R)。本人的python3.6.1版本。要使这个科学计算环境可以在计算机上运行要求就是能够多个Python版本并存,尤其是2.x和3.x的并存。可以利用 conda 的虚拟环境管理功能在 Python2 和 Python3 之间自由切换。鉴于目前的情况,需要再配置一个python2.7版本环境。
[AppleScript] 纯文本查看 复制代码
# 基于 python2.7 创建一个名为python27的环境
conda create --name python27 python=2.7
#列出所有的环境,来检查一下截至目前所安装的环境,conda将会显示所有环境的列表,当前环境会显示在一个括号内,有时也会在目前活动的环境前边加上*号。
conda info --envs
# 激活python27环境
source activate python27
#如果当前已经激活了某个Python环境,那么就可以在当前环境开始安装第三方包。
conda install rseqc
#升级rseqc
conda update rseqc
#如果不再使用当前环境,可以选择关闭当前工作环境
source deactivate
#下载bed文件,基因组覆盖率的QC需要提供bed文件,可以直接RSeQC的网站下载,或者可以用gtf转换
wget [url=http://dldcc-web.brc.bcm.edu/lilab/liguow/RSeQC/dat/hg19_RefSeq.bed.gz]http://dldcc-web.brc.bcm.edu/lil ... /hg19_RefSeq.bed.gz[/url]
#更换目录
cd /mnt/d/reference/bed/
#解压
gunzip hg19_RefSeq.bed.gz
#定位到安装目录
cd /home/huiguang/anaconda3/bin/
#使用bam_stat.py
bam_stat.py -i /mnt/d/AKAP95/aligned/SRR3589956.sorted.bam

[AppleScript] 纯文本查看 复制代码
# read_distribution.py
read_distribution.py -i /mnt/d/AKAP95/aligned/SRR3589956.sorted.bam -r /mnt/d/reference/bed/hg19_RefSeq.bed
[AppleScript] 纯文本查看 复制代码
# clipping_profile.py
clipping_profile.py -i /mnt/d/AKAP95/aligned/SRR3589956.sorted.bam -s "PE" -o out
ls -l out.clipping_profile.* |awk '{print $9}'
Tips: null device应该是说,没有图形化输出设备。因为已经将图形生成pdf,所以,这个报错没有太大问题。
[AppleScript] 纯文本查看 复制代码
#移动文件
mv out.clipping_profile.* /mnt/d/AKAP95/
结果包括一个R文件,两个pdf图文件和一个原始Excel文件。其中生成的R文件是用于做pdf图的,如果选择的是“SE”则会生成一个pdf文件。本例参数为PE,所以有两个文件。pdf图形如下:
[AppleScript] 纯文本查看 复制代码
# deletion_profile.py
deletion_profile.py -i /mnt/d/AKAP95/aligned/SRR3589956.sorted.bam -l 51 -o out
ls -l out.deletion_profile* |awk '{print $9}'
mv out.deletion_profile.* /mnt/d/AKAP95/
[AppleScript] 纯文本查看 复制代码
# geneBody_coverage.py
geneBody_coverage.py -i /mnt/d/AKAP95/aligned/SRR3589956.sorted.bam -r /mnt/d/reference/bed/hg19_RefSeq.bed -o output
ls -l output.geneBodyCoverage* |awk '{print $9}'
mv output.geneBodyCoverage* /mnt/d/AKAP95/



Tips: geneBody_coveraged可以同时输入多个bam文件,将所有的结果生成一张图。
8.     IGV截图

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:一行脚本判断你的fastq测序数据的质量编码方式
下一篇:unbuntu 安装mysql 报错E:invalid operation install,求助!!
回复

使用道具 举报

9

主题

22

帖子

168

积分

注册会员

Rank: 2

积分
168
发表于 2017-8-28 11:33:38 | 显示全部楼层
非常好,学习啦!
回复 支持 反对

使用道具 举报

3

主题

20

帖子

336

积分

中级会员

Rank: 3Rank: 3

积分
336
发表于 2017-11-23 15:32:19 | 显示全部楼层
mark下!
回复

使用道具 举报

3

主题

20

帖子

336

积分

中级会员

Rank: 3Rank: 3

积分
336
发表于 2017-11-24 15:56:22 | 显示全部楼层
hg19.tar.gz和mm10.tar.gz这两个文件下载的时候直接没有速度,醉了!
回复 支持 反对

使用道具 举报

3

主题

7

帖子

98

积分

注册会员

Rank: 2

积分
98
发表于 2018-1-20 20:37:40 | 显示全部楼层
做个记号,谢谢。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-7-24 09:17 , Processed in 0.052350 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.