生物信息学笔记6:hisat2比对生成sam文件

在上一步,我们建立了index,可以在运行index命令的目录(我是在dnaseq\chrX_data\third\index,这个dnaseq目录是win10和linux的共享文件夹,关于该共享文件夹可能有些人可能还会有点弯路要走,我一开始是在linux下面找不到这个共享文件夹,这个后面如果有时间我会专门写一篇)下找到相应的文件如下图(里面包括了chrX.exon和chrX.ss这两个从gtf文件中提取的,以及剩下的就是index文件,有1-8,但是在比对时用chrX_tran代表):

然后,我们进入samples目录(就是我们下载了序列的目录),如下图:

注意了,上面图片示例都是在win10下截图的,只是为了更好理解,实际这一步还是在linux下操作比较方便,因为接下来要建立一个批处理文件samAA.sh是在linux下运行的。如下图:

建立方法,在linux,samples文件夹,输入命令:

touch samAA.sh

然后编辑该sh文件:

nano samAA.sh

在编辑其中输入命令:

for i in *1.fastq.gz;
do
i=${i%1.fastq.gz*};

nohup hisat2 -p 1 --dta -x /mnt/hgfs/dnaseq/chrX_data/third/index/chrX_tran -1 ${i}1.fastq.gz -2 ${i}2.fastq.gz -S /mnt/hgfs/dnaseq/chrX_data/third/sam-dta/${i}align.sam 2>/mnt/hgfs/dnaseq/chrX_data/third/sam-dta/${i}align.log &
done

保存退出后,在该文件夹下直接输入“sh samAA.sh”运行循环比对(多条序列比对)。

这个命令为什么这么写:

for语句第一行是要取出该目录下序列的文件名,就是1.fastq.gz前面的名字,比如“ERR188428_chrX_”,

do 开始循环,只要存在这样的i就继续执行,以其中循环一遍为例,比如让i=ERR188428_chrX_,

nohup开始在后台运行比对,比如把log保存在“/mnt/hgfs/dnaseq/chrX_data/third/sam-dta/ERR188428_chrX_align.log”里面,这样可以在文件里面查看是否成功,如果不成功问题提示如何。

hisat2 语句中“-p 1”是指单一进程运行比对,可以跟进电脑配置自行调整;

“--dta”是指输出的文件给stringtile继续使用;

“-S”注意这里是大写S,千万注意,这是指定位置输出比对好的sam文件,不要搞成小写s。

好了,比对结果如下图(看一下目录名,和命令行比较):

解释一下:命令行中有“/mnt/hgfs/”,后面才跟上dnaseq目录(这个是两个系统的共享目录),这个是因为linux下面的共享目录要挂载在/mnt/hgfs/下面。

因为比对时随机的,如果你用notepad++打开sam文件会发现里面的内容每次都是不一样的(但总大小一样,因为随机比对,顺序每次都不同),所以sam文件比对好以后还要进行排序,今天就到这里了。

生物信息学笔记5:用hisat2软件包建立基因组index

顾名思义,建立基因组索引,主要是提高比对的速度、效率,对剪切位点进行预测,hisat2建立基因组+转录组+SNP索引,so,为什么要建立索引:

高通量测序有成千上万条reads需要高效比对到参考基因组上,并且保证一定的准确率,答案不一定说完全正确,但一定要非常接近真实数据。需要根据参考基因组序列,经过一定算法(大部分情况是BWT或其改良算法)转换成index,把reads通过和index的比较过程进行回贴(maping到参考基因组),大幅度缩短比对maping的时间。(关于BWT算法,请参考:http://www.bio-info-trainee.com/?s=bowtie和http://www.biotrainee.com/thread-26-1-1.html,后期考虑转载过来学习以防链接失效)

步骤:

1、从注释文件里面抽取出剪切位点和外显子信息,用到hisat2文件夹下面的extract_splice_sites.py和extract_exons.py

比如,注释文件在此路径:chrX_data/genes/chrX.gtf

这个gtf文件最好自己从Ensembl上下载,或者Entrez上下载。

则命令行:

extract_splice_sites.py chrX_data/genes/chrX.gtf >chrX.ss

extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon

2、用hisat2建立index

hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran

3、当需要使用snp信息时,命令行需加入下面几行

extract_snps.py snp142Common.txt > genome.snp  #####注释:本行井号后面请勿复制,本命令主要是将snp信息文件转换成hisat2可以使用的文件

hisat2-build -p4 genome.fa --snp genome.snp --ss genome.ss --exon genome.exon genome_snp_tran

 

特别忠告:如果是用chrX练手,而且你的电脑内存大于10G,那可以尝试自己建立index,其他情况,请下载现有index,因为这步工作内存小干不了,内存大也需要时间去完成(是很多很多很多时间)。