从零开始生物信息学(1):WGS初探 |
您所在的位置：网站首页 › wg是什么意思医学 › 从零开始生物信息学(1):WGS初探 |
从零开始生物信息学(1):WGS初探
|
前言 最近开始自学生物信息学,由于本科和研究生都是计算机专业,对生物知识了解甚少,所以最近恶补了不少生物和全基因组测序的相关知识,并且把整个过程记录下来,这是一个漫长艰巨的过程,也希望自己能够一直坚持下来. WGS全基因组测序（WGS, Whole Genome Sequencing）是下一代测序技术，用于快速，低成本地确定生物体的完整基因组序列。基因组的深度测序对于临床研究的意义重大，解读WGS数据并了解基因组突变在健康和疾病中的重要性是精准医疗的基石. WGS分析流程能分为三大块，数据处理、检测变异和综合分析,目前主要可以分为以下几个目的: 质量控制. 也就是数据质控,在目前的高通量测序方法中,会发生低质量污染读数,会影响接下来的数据分析,因此质控对于WGS的数据质量保证十分重要. 基因组对齐.本质上就是大规模的字符串匹配问题,求解两个序列的相似度, 可以通过动态规划的方法以及大名顶顶的blast算法来解决. 变异检测. 其中先搞清楚变异的种类,目前主要包括几种: 单碱基变异，即单核苷酸多态性(SNP) 较短的连续的插入(Insertion)或者缺失(Deletion),简称 Indel,这里的较短指的是不超过50 base pair(bp)的线性长度 基因组结构性变异（SVs, Structure Variantions,这个就相对复杂,包括50bp以上的长片段序列插入或者删除（Big Indel）、串联重复（Tandem repeate）、染色体倒位（Inversion）、染色体内部或染色体之间的序列易位（Translocation)等等,可视化如下: 而所谓的变异检测就是检测说上面所述的变异类型. 当然基因测序也不是一蹴而就的,经历的漫长的发展,目前最流行的第三代测序技术(NGS, Next Generation Sequencing),之前两代测序技术对于社会的贡献也是十分巨大,这里就不展开讲述了,可以参考我另外一份笔记基本技术的发展, 这部分笔记主要在于实操,对于现有的数据进行分析. 正文这个笔记主要是记录基因数据的的处理流程以及熟悉整体的操作流程. 首先,目前大多数的基因分析的软件和包都是在linux环境下进行的,所以不了解linux的童鞋也可以先了解下linux的操作,由于本人计算机专业,所以堆linux也用的比较多,这个压力不大. 在这次笔记中主要用到以下几个工具包: samtools ,这是一个用于操作原始数据的工具,包括堆基因数据的查看,排序和合并等等,十分强大 bwa, 这是一个序列比对到参考基因组上的工具,能够快速的进行基因组比对.GATK, 强大的工具包,用于分析变异信息sratoolkit,用于将NCBI数据装换的工具包tabix, 压缩和解压工具安装samtools,网上有很多教程,需要安装一大堆依赖,其实只要在终端输入:sudo apt-get install samtools bwa,访问bwa官网,选择最新版本下载,解压,终端直接进入文件夹目录,直接make一下,就会生成一个bwa的可执行文件. GATK可以参照官网的安装方法,目前已经更新到4.0的版本,进入官网下载,解压之后,将文件目录加入系统路径即可,如果不知道怎么加可以参考这个教程 sratoolkit, 可以直接在官网下载最新的版本,或者直接执行下面命令: wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.2/sratoolkit.2.9.2-centos_linux64.tar.gz tar xzvf sratoolkit.2.9.2-centos_linux64.tar.gz 然后同样加目录加入到系统环境变量中即可 5 tabix,是一个高效压缩和解压的工具包,安装也相当简单: sudo apt-get install tabix 数据准备由于目前人类全基因组序列的长度约为3Gb,在单机电脑可能无法快速的进行WGS的数据分析,因此这里就直接采用比较简单的数据,大肠杆菌的数据作数据分析,是为了明白整个WGS数据分析的流程,以后对于人类基因组数据也可以采用相同的策略 下载数据这种基因叫做E.coli K12,基因组只有4.6Mb,数据直接可以NCBI (National Center for Biotechnology Information)这个数据库网站中找到,直接下载到本地,然后解压重命名为.fa的文件: wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz gzip -dc GCF_000005845.2_ASM584v2_genomic.fna.gz > E.coli_K12_MG1655.fa这个是一个fastq的文件,你可以执行 head E.coli_K12_MG1655.fa 就可以看到以下的内容: >NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genomeAGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGAC第一行是这个基因组的信息,第二行开始就是原始的ACGT的碱基数据了,因此十分的直观易懂. 然后采用samtools工具堆这份数据建立一个索引: samtools faidx E.coli_K12_MG1655.fa 会在相同目录中生成一个E.coli_K12_MG1655.fa.fai,这个索引文件可以高效的对原始数据进行位点的查询,是一个十分实用的功能. 有了原始数据,我们还需要测序的数据,可以直接在NCBI中下载Illumina MiSeq测序平台的SRR1770413测序数据. 首先,也是下载数据: wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR177/SRR1770413/SRR1770413.sra但是这个格式是测序平台的压缩格式,需要转化为统一的fastq格式用于后续处理.这时就用到NCBI的官方工具包sratoolkit cd ./ #进入sratookit 的目录 fastq-dump --split-files SRR1770413.sra #将SRA文件转换为fastq文件这时会在本地生成两个文件,分别为SRR1770413_1.fastq和SRR1770413_2.fastq,这是E.coli K12数据的两份read文件(应该是一份从左到右,一份从右到左) 然后对着两份文件进行压缩一下,当然不压缩也是可以的: bgzip -f SRR1770413_1.fastq bgzip -f SRR1770413_2.fastq现在给大家看看我的目录结构,分别有以下的东西: ./bio/ ├── input │ ├── E.coli_K12_MG1655.fa │ ├── E.coli_K12_MG1655.fa.fai │ ├── SRR1770413_1.fastq.gz │ └── SRR1770413_2.fastq.gz └── output测序文件比对基因的比对也是生物信息学一个很重要的应用, 本质就是把我们测序生成的数据匹配到参考基因组中, 样本长度越大,所需要的时间也就越多.这里由于不涉及到原理,直接使用强大的bwa的包来进行, 首先是构建比对序列的比对索引,相当于给每个测序文件一个序号: ./your_bwa_file/bwa index ./input/E.coli_K12_MG1655.fa这里用的是之前说过的编译过后的bwa可执行文件,然后可以看到执行过程 [bwa_index] Pack FASTA... 0.04 sec [bwa_index] Reverse the packed sequence... 0.01 sec [bwa_index] Construct BWT for the packed sequence... [bwa_index] 0.40 seconds elapse. [bwa_index] Construct BWT for the reverse packed sequence... [bwa_index] 0.39 seconds elapse. [bwa_index] Update BWT... 0.01 sec [bwa_index] Update reverse BWT... 0.01 sec [bwa_index] Construct SA from BWT and Occ... 0.13 sec [bwa_index] Construct SA from reverse BWT and Occ... 0.13 sec之后会在E.coli_K12_MG1655.fa 目录下生成一些文件: ├── E.coli_K12_MG1655.fa ├── E.coli_K12_MG1655.fa.amb ├── E.coli_K12_MG1655.fa.ann ├── E.coli_K12_MG1655.fa.bwt ├── E.coli_K12_MG1655.fa.fai ├── E.coli_K12_MG1655.fa.pac ├── E.coli_K12_MG1655.fa.rbwt ├── E.coli_K12_MG1655.fa.rpac ├── E.coli_K12_MG1655.fa.rsa ├── E.coli_K12_MG1655.fa.sa接着使用bwa进行数据比对: ./your_bwa_file/bwa mem -R '@RG\tID:foo\tPL:illumina\tSM:E.coli_K12' ./input/E.coli_K12_MG1655.fa ./input/SRR1770413_1.fastq.gz ./input/SRR1770413_2.fastq.gz | samtools view -Sb - > ./output/E_coli_K12.bam这里相当于将测序文件定位在参考基因组上进行之后再用linux的管道和samtools工具包把输出转化乘bam格式并保存到本地.如果对bwa men 命令不熟悉 也可以输入bwa mem -h来看参数设置,这里上网查了下, -R参数是必须的参数,其中PL参数也不是随便设置的,需要在指定的测序平台上选一个,目前是从ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN里面选择一个.这个需要注意一下 然后对bam文件进行排序,排序的作用就是便于后期的数据删除和下一步的匹配: samtools sort ./output/E_coli_K12.bam ./output/E_coli_K12.sort默认匹配是按照参考序列位置从小到大进行排序,但是这时候由于测序时采用了PCR扩增,高中应该学过,就是复制多个基因组序列进行测序,会有大量的重复的测序数据,因此需要去重,这时就用到了GATK的工具包用于标记和去重: gatk MarkDuplicates -I ./output/E_coli_K12.sort.bam -O ./output/E_coli_K12.sort.markdup.bam -M ./output/E_coli_K12.sort.markdup_metrics.txt最后还是用samtools对标记重复后的数据创建索引: samtools index ./output/E_coli_K12.sort.markdup.bam变异检测这里用GATK工具包完成变异检测,首先为参考序列产生一个字典文件,并且存放在input的目录下(为了保证和fastq文件在同一目录下),在接下来变异检测中需要用到 gatk CreateSequenceDictionary -R ./input/E.coli_K12_MG1655.fa -O ./input/E.coli_K12_MG1655.dict首先用HaplotypeCaller工具包生成中间的gvcf文件 gatk HaplotypeCaller -R ./input/E.coli_K12_MG1655.fa --emit-ref-confidence GVCF -I ./output/E_coli_K12.sort.markdup.bam -O ./output/E_coli_K12.g.vcf 需要几分钟的时间(感觉GATK这个工具包还是有点慢),然后进行变异检测,生成一个的文件 gatk GenotypeGVCFs -R ./input/E.coli_K12_MG1655.fa -V ./output/E_coli_K12.g.vcf -O ./output/E_coli_K12.vcf最后,同样的,对这个VCF文件进行压缩，并用tabix为它构建索引 bgzip -f ./output/E_coli_K12.vcf tabix -p vcf ./output/E_coli_K12.vcf.gz这是最终的目录: ./bio/ ├── input │ ├── E.coli_K12_MG1655.dict │ ├── E.coli_K12_MG1655.fa │ ├── E.coli_K12_MG1655.fa.amb │ ├── E.coli_K12_MG1655.fa.ann │ ├── E.coli_K12_MG1655.fa.bwt │ ├── E.coli_K12_MG1655.fa.fai │ ├── E.coli_K12_MG1655.fa.pac │ ├── E.coli_K12_MG1655.fa.rbwt │ ├── E.coli_K12_MG1655.fa.rpac │ ├── E.coli_K12_MG1655.fa.rsa │ ├── E.coli_K12_MG1655.fa.sa │ ├── SRR1770413_1.fastq.gz │ └── SRR1770413_2.fastq.gz └── output ├── E_coli_K12.bam ├── E_coli_K12.g.vcf ├── E_coli_K12.g.vcf.idx ├── E_coli_K12.sort.bam ├── E_coli_K12.sort.markdup.bam ├── E_coli_K12.sort.markdup.bam.bai ├── E_coli_K12.sort.markdup_metrics.txt ├── E_coli_K12.vcf.gz ├── E_coli_K12.vcf.gz.tbi └── E_coli_K12.vcf.idx欢迎大家关注我的知乎专栏：从零开始生物信息学 相同内容也可以关注我的微信公众号: 壹读基因： 
|
【本文地址】
今日新闻 |
推荐新闻 |