一. 简介
在基因组数据分析过程中,我们可以将fastq,fasta,ONT等不同格式的数据mapping到参考基因组上,然后通过IGV进行可视化,结合GFF3文件手工对感兴趣的区段进行校正,mapping,工具很多,这里主要介绍李恒开发的BWA和Minimap2。
二. 软件安装
- 首先在服务器安装miniconda后,建python2.7的环境,修改~/.bashrc文件方便进入
[abc@Server Bwa]$ vi ~/.bashrc
alias py27='. activate;conda activate py27' #.bashrc添加一行
[abc@Server Bwa]$ py27 #每次登录只需输入py27即可进入python2.7环境
(py27) [abc@Server Bwa]$ #前面显示环境py27
- 然后,在bioconda网页分别搜BWA、Minimap2和samtools进行软件安装,
(py27) [abc@Server Bwa]$ conda install -c bioconda bwa #安装bwa
(py27) [abc@Server Bwa]$ conda install -c bioconda minimap2 #安装minimap2
(py27) [abc@Server Bwa]$ conda install -c bioconda samtools #安装samtools
三. Mapping数据到参考基因组
- 二代数据,建议用bwa mem:
(py27) [abc@Server Bwa]$ nohup bwa index ref.fas 2>index.log & #建立索引
(py27) [abc@Server Bwa]$ nohup bwa mem ref.fas a1.fq a2.fq 1>a.sam 2>a.sam.log & #mapping fastq数据a1.fq和a2.fq(paired-end)到参考基因组ref.fas
(py27) [abc@Server Bwa]$ nohup bwa mem ref.fas b.fq 1>b.sam 2>b.sam.log & #mapping fastq数据b.fq(单末端测序文件)到参考基因组ref.fas
(py27) [abc@Server Bwa]$ nohup bwa mem ref.fas c.fasta 1>c.sam 2>c.sam.log & #mapping fasta数据c.fasta(可以是基因组,cds)到参考基因组ref.fas
(py27) [abc@Server Bwa]$ nohup bwa mem -x pacbio/ont2d ref.fas a1.fq a2.fq 1>a.sam 2>a.sam.log & #对于超长读长的reads,比如PacBio和Nanopore测序仪产生的reads,该命令容易报错,估计和bwa版本有关
(py27) [abc@Server Bwa]$ nohup samtools view -bS a.sam 1> a.bam 2>a.bam.log& #将sam文件转为bam文件
(py27) [abc@Server Bwa]$ nohup grep 'LG02' a.sam >a_LG02.sam& #如果sam文件过大可以只取其中一个染色体数据进行可视化
(py27) [abc@Server Bwa]$ nohup samtools sort a.bam a.sort 2>a.sort.bam.log& #对bam文件进行sort,生成a.sort.bam文件
(py27) [abc@Server Bwa]$ nohup samtools merge a.sort.merge.bam a1.sort.bam a2.sort.bam a3.sort.bam ...& #可以对多个sort过的bam文件进行合并
(py27) [abc@Server Bwa]$ nohup samtools index abc.sort.merge.bam abc.sort.bam.merge.bai & #对sort.bam文件建立索引,生成对应bai文件
可以将bam文件和对应的bai文件拷到台式机,利用IGV进行可视化。
- 三代数据,建议用minimap2:
(py27) [abc@Server Bwa]$ nohup minimap2 -d target.mmi target.fas & #建立索引
(py27) [abc@Server Bwa]$ nohup minimap2 -a target.mmi a.fas >a.sam 2>a.sam.log & # mapping fasta数据到参考基因组
(py27) [abc@Server Bwa]$ nohup minimap2 -a a.fas b.fas >b.sam 2>b.sam.log & # 不建索引直接mapping b.fas到a.fas
(py27) [abc@Server Bwa]$ nohup minimap2 -x map-ont -d target.mmi a.fas >a.sam 2>a.sam.log& # mapping ont数据a.fas到参考基因target
(py27) [abc@Server Bwa]$ nohup minimap2 -ax map-ont ref.fas ont.fq.gz >ont.sam 2>ont.sam.log& # mapping ont数据ont.fq.gz到参考基因ref.fas
生成的sam文件参照前面samtools相关命令,转成sort.bam和sort.bai文件,拷到台式机上面利用IGV进行可视化。
四. 利用IGV进行可视化
下载安装http://software.broadinstitute.org/software/igv/download
打开IGV,Genomes->Load Genome from File加载基因组数据,File->Load from File加载bam文件和GFF3注释文件,bam文件需要有对应的bai文件存在于同一文件夹中:
注意:在加载GFF3文件的时候,需先去掉GFF3文件中包含gene和mRNA的行,不然会将mRNA和gene识别为UTR,利用窄长方形表示,看不到内含子的直线
(py27) [abc@Server Bwa]$ more abc.gff3|grep -v 'gene'|grep -v 'mRNA' >abc_1.gff3
可以通过选择加右键对每个Track的参数进行调整,如coverage的范围,序列显示形式等,下面是不同数据加载的效果图: