Samtools截取基因组序列(Docker篇)

Samtools截取基因组序列(Docker篇)

Samtools是一个用于操作sam和bam文件的工具软件,能够对比对文件进行二进制查看、格式转换、排序及合并等,结合sam格式中的flag、tag等信息,还可以完成比对结果的统计汇总,是处理sam和bam文件(例如:转录组Tophat分析软件输出的比对结果为.bam文件,而重测序中BWA、bowtie等比对软件则主要输出为.sam文件)不可或缺的神器!

Samtools软件的安装相当麻烦,安装教程可参考:samtools你安装成功了吗?那有没有方法可以省去繁琐的安装步骤而直接使用samtools呢?答案是有的!Docker工具可以完美的解决这个问题!

Docker是生信初学者必不可少的工具,具体介绍参见文章:生物信息软件安装解决方案-docker虚拟化技术Docker工具下载安装请参考:Docker 工具安装(windows Toolbox版本)详解版!Docker 工具安装(windows 桌面版)详解版!,这里便不多介绍了。

下载镜像并运行

安装好Docker后,搜索我们omicsclass提供的docker镜像,其中便有安装好samtools软件的镜像。

$ docker search omicsclass
NAME                           DESCRIPTION                                     STARS               OFFICIAL            AUTOMATED
omicsclass/gene-family         gene-family analysis docker image               2
omicsclass/isoseq3             isoseq3 v3.3.0 build by omicsclass              0
omicsclass/bwa                 BWA v0.7.17 build by omicsclass                 0
omicsclass/rnaseq              RNA-seq analysis docker image build by omics…   0
omicsclass/samtools            samtools v1.10 build by omicsclass              0
omicsclass/biocontainer-base   Biocontainers base Image centos7                0
omicsclass/blast-plus          blast+ v2.9.0                                   0
omicsclass/blastall            legacy blastall v2.2.26                         0
omicsclass/sratoolkit          SRAtoolkit v2.10.3 and aspera v3.9.9.177872     0

下载镜像。

$ docker pull omicsclass/samtools
Using default tag: latest
latest: Pulling from omicsclass/samtools
ab5ef0e58194: Already exists
417469905807: Already exists
ed09842cc19f: Already exists
f860268ff83f: Already exists
f87dd41136a6: Already exists
90091b4f5d91: Already exists
6485f44fc594: Pulling fs layer
latest: Pulling from omicsclass/samtools
ab5ef0e58194: Already exists
417469905807: Already exists
ed09842cc19f: Already exists
f860268ff83f: Already exists
f87dd41136a6: Already exists
90091b4f5d91: Already exists
6485f44fc594: Pull complete
Digest: sha256:e641dd5b9f60d8f9d01f0d109eff72d15836d0d59a753e3e35677b1200adc4a1
Status: Downloaded newer image for omicsclass/samtools:latest

下载完成后可以检查一下。

$ docker images
REPOSITORY              TAG                 IMAGE ID            CREATED             SIZE
omicsclass/samtools     latest              9373e18781bf        8 days ago          2.04GB
omicsclass/blast-plus   latest              0220cac51a6e        8 days ago          2.55GB

之后在电脑的D盘创建samtools文件夹,并粘贴需要的染色体fasta文件。如果docker是windows Toolbox版本,需要挂载D盘到虚拟机,具体操作可参考 Docker 工具安装(windows Toolbox版本)详解版!

进入虚拟机。

$ docker run --rm -v /d/samtools:/work -it omicsclass/samtools:latest   #这里我使用的Docker是windows Toolbox版本
######################################################
#      欢迎使用组学大讲堂提供的docker镜像              #
#      问题交流请访问:www.omicsclass.com             #
######################################################

            Linux新手建议学习课程:
      --> https://www.omicsclass.com/article/702

         搭建实验室生信分析平台与docker使用详情见课程:
      --> https://www.omicsclass.com/article/1181
[root@0e9f42f25cc1  10:16:46 /work]#

查看工作目录。

ll
total 117M
-rwxrwxrwx 1 1000 ftp 117M Apr 23 10:19 Arabidopsis_thaliana.TAIR10.31.dna.toplevel.fa

samtools提取序列

samtools faidx 能够对fasta序列建立一个后缀为.fai 的文件,根据这个.fai 文件和原始的fasta文件,能够快速的提取任意区域的序列。用法:

samtools faidx Arabidopsis_thaliana.TAIR10.31.dna.toplevel.fa

该命令对输入的fasta序列有一定要求:对于每条序列,除了最后一行外,其他行的长度必须相同。

>one 
ATGCATGCATGCATGCATGCATGCATGCAT 
GCATGCATGCATGCATGCATGCATGCATGC 
ATGCAT 
>two another chromosome 
ATGCATGCATGCAT 
GCATGCATGCATGC 

最后生成的.fai文件如下, 共5列,\t分隔;

Pt      154478  4       60      61
Mt      366924  157061  60      61
4       18585056        530104  60      61
2       19698289        19424914        60      61
3       23459830        39451511        60      61
5       26975502        63302342        60      61
1       30427671        90727439        60      61

第一列 NAME:序列的名称,只保留“>”后,第一个空白之前的内容;

第二列 LENGTH:序列的长度,单位为bp;
第三列 OFFSET:第一个碱基的偏移量,从0开始计数,换行符也统计进行;
第四列 LINEBASES:除了最后一行外,其他代表序列的行的碱基数,单位为bp;
第五列 LINEWIDTH : 行宽,除了最后一行外,其他代表序列的行的长度,包括换行符,在windows系统中换行符为\r\n, 要在序列长度的基础上加2;

最后指定要提取的序列或序列在染色体上的位置,即可提取出需要的序列:

samtools faidx Arabidopsis_thaliana.TAIR10.31.dna.toplevel.fa Pt> Pt.fa             #提取叶绿体序列
samtools faidx Arabidopsis_thaliana.TAIR10.31.dna.toplevel.fa Pt:100-1000> Pt.fa    #提取叶绿体第100到1000碱基的序列


好了,今天先介绍到这里,想了解更多生物信息学相关知识,请持续关注本公众号!祝您科研愉快!

  • 发表于 2020-04-28 13:49
  • 阅读 ( 6003 )
  • 分类:软件工具

0 条评论

请先 登录 后评论
安生水
安生水

347 篇文章

作家榜 »

  1. omicsgene 698 文章
  2. 安生水 347 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. xun 82 文章
  6. 红橙子 78 文章
  7. rzx 74 文章
  8. CORNERSTONE 72 文章