生物信息学笔记7:对sam文件重新排序

通过hisat2比对得到sam文件,如果重复多次用同一个命令比对同一条序列,每次得到的sam文件其实是大小一样但内容顺序不同的文件,所以我们需要用软件进行重新排序。

在sam文件所在目录建立bam1.sh文件:

touch bam1.sh

nano bam1.sh

在文本编辑器中输入:

for i in *.sam;
do
i=${i%_align.sam*};
nohup samtools sort -@ 1 -o ../bam-dta/${i}.bam ${i}_align.sam 2>/mnt/hgfs/dnaseq/chrX_data/third/bam-dta/${i}bam.log &
done

保存后,运行该sh文件:

sh bam1.sh

运行完毕后在bam-dta目录下生成排序好的bam文件,如下图:

生物信息学笔记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,因为这步工作内存小干不了,内存大也需要时间去完成(是很多很多很多时间)。

生物信息学笔记4:hisat2的下载安装及环境配置

一、需要分析的序列,比如在samples文件夹里:

ERR188044_chrX_1.fastq.gz

ERR188044_chrX_2.fastq.gz

ERR188104_chrX_1.fastq.gz

ERR188104_chrX_2.fastq.gz

......

我就不一一列出来,这里列出来主要就是想说同一个序列有1和2,就是双端测序的意思

 

二、ubuntu18.04下面hisat2的下载和安装,这个差点忘记写。

先找个目录,比如建立个biosoft的目录,专门用于放生物信息相关软件。

1、命令行先建立后进入该目录:

mkdir ~/biosoft

cd ~/biosoft

2、进入目录后下载hisat2软件:

wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip

3、解压缩下载的zip压缩包(解压后获得hisat2-2.1.0的文件夹):

unzip hisat2-2.1.0-Linux_x86_64.zip

4、改文件目录名(我没有改名字,这步骤可以忽略,也可以把文件夹名字改短点):

mv hisat2-2.1.0 hisat2

5、加环境变量:

很多地方用vi后者vim命令,我用不习惯,我觉得用nano比较方便:

nano ~/.bashrc

然后出现了如下界面,移动光标到末尾,加入“PATH=$PATH:~/biosoft/hisat2-2.1.0:........”:

输入后,保存按ctrl+o,回车,然后ctrl+x退出。

请不要直接用export命令,该命令只能在当前登录窗口有效,断开连接后再次用同一个用户登录,也会重置环境。

保存好后,还需要让环境生效:

source ~/.bashrc

6、环境配置好后可以在任意文件目录输入:hisat2,然后出现如下信息说明配置成功:

 

 

本来想把hisat2的下载安装使用全放在一起,后来转念一想,分开的话以后搜索起来会更方便点。

生物信息学笔记3:练手序列的准备

本文参考了简书:RNA-seq数据分析---方法学文章的实战练习

我的启蒙文章,现在的生物信息学和以前的比真的是牛逼太多了,以前搞个blast、phylip等就好像差不多了,而实际上那些都算是比较简单的基本操作了。真正要搞懂的不单是碱基的差异,还有表达基因的差异、蛋白功能的差异......

简书中此文已经罗列了为什么用hisat2和stringtile,本文不赘述,反正就是这两个软件的算法可以用最少的时间和最低的资源占用,获得更好更优质的比对结果。

一、练手序列。最好入手的应该还是人类X染色体的序列文件的比对分析,因为X染色有男女的差异、人种的差异,找12个样本(人种各6个,不同人种中男女各3个),算是有一定的生物学重复吧。

序列下载可以有多个途径(建议在win10下用迅雷下载,不然你会疯的,当然和不同网络有关系):

1、linux命令行,很牛逼的样子,可是如果网络不在一个频道,这逼装的再牛也没用,下载速度比蜗牛还不如啊,重新连接次数多了还会自动停止。

命令行如下(这是约翰霍普金斯大学服务器上的数据,练手用的序列):

nohup wget ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.gz 2>download.log &
tar zxvf chrX_data.tar.gz

2、迅雷下载,这个不用细说了吧,如果这个不会,建议还是不要学生物信息了。

二、练手序列找好,还要找个参考基因序列和对应参考注释。这个最好自己去下载Ensembl版本,不要用hisat2官方提供的,因为后面涉及到的go分析等等都会用到基因名或编号,这个都需要与网络数据库配套,如果下载一个序列的编号方式不符合这些数据库的标准,那分析很难进行下去,有些编号是转换不了的。Ensembl版本全基因组的注释文件下载

三、目的就是把练手序列map到参考序列上(这个是公认的当前最完整和准确的序列,并不断完善中),参考注释是对应于参考序列的、用来标注出相应序列对应的基因的文件。通过把练手序列map到参考序列,分析出这些练手序列里基因被检测到的次数(因为练手序列通过高通量测序后的结果,这是一个碎片化的随机检测过程,而测序的样本是dna转录成mrna后,反转录的cdna,也就是相当于一个表达文库,那么比然会有基因序列的表达量的差异,这个是我的理解过程,可能会有失偏颇或者谬以千里,欢迎批评指正,共同进步),然后通过分析基因表达的次数(这个可能需要一定的统计学纠正等等),统计出一个接近于真实情况的基因表达情况,进而分析基因表达的差异,同时分析基因和蛋白功能,意图寻找不同细胞的差异所在(比如癌细胞和正常细胞)。

四、另外一个目的就是寻找未知isoforms,就是为了寻找不同基因用的,也就是当以上差异并不能很有效的说明问题关键或者需要更多的解释的时候,可以考虑寻找未知isoforms,很容易做无用功,但是一旦找到未知的关键isoforms,很容易产生重大突破。

生物信息学笔记2:主机及环境配置

主机配置:e3-1230v2+24g内存。原来用的是e3-1230v2+8g内存的主机,以为这样自己学习捣鼓捣鼓够了,可是跑起来明显内存不够,于是加内存到24g,这样小一点的序列,比如X染色体这种的是能跑起来了

系统配置:主机win10+VMware(Ubuntu 18.04 live sever amd64),用惯了win10,你让主机我装个linux太浪费了,然后分配10g内存给ubuntu系统,设置主机和虚拟机的共享文件夹,比如e:\Vhost\dnaseq,这里尽量别用大写,因为linux系统是区分大小写的,命令行输入比较麻烦。共享文件夹的建立是为了方便在win10界面下操作linux下生成的文件,这样比命令行输入效率更高。

生物信息学笔记1

一直想写个生物信息学学习的笔记或总结,今天总算是开个头。

环境配置:ubuntu14.04 lts

软件:hisat2、stringtile、R语言(包括各种package,在后面用到时候会写)、GSEA及其R语言包、R语言的FactoMiner(我感觉聚类分析有点没啥用)

思路:测序结果用hisat2比对,bam排序,stringtile计算counts、丰度等,用R语言进行各种分析(DESeq2-go-kegg、GSEA、Facto等)

很懒,写的可能只有我自己能看懂,今天第一篇,就为开个头。