pysam - 多种格式基因组数据(sam/bam/vcf/bcf/cram/…)读写与处理模块(python)

在开发基因组相关流程或工具时,经常需要读取、处理和创建bam、vcf、bcf文件。目前已经有一些主流的处理此类格式文件的工具,如samtools、picard、vcftools、bcftools,但此类工具集成的大多是标准功能,在编程时如果直接调用的话往往显得不够灵活。

本文介绍的是一个处理基因组数据的python模块,它打包了htslib-1.3、samtools-1.3 和 bcftools-1.3的核心功能,能在编程时非常灵活的处理bam和bcf文件。

以下主要介绍pysam的安装和使用方法:

1. 安装

如果Linux上安装了pip,可以一键安装,在集群上的话,需要登录安装节点进行安装。

pip3 install pysam

检查是否安装成功

import pysam

2.读取bam文件(pysam.AlignmentFile)

bam是sam的二进制文件,因其占用空间少,所以都会使用bam进行存储和操作。

要读取bam文件,必须先创建一个AlignmentFile对象.

path_in = ‘./test.bam‘
samfile = pysam.AlignmentFile(path_in, "rb")

之后就可以逐行读取和处理bam文件了(顺序读取),以下打印出了bam的一行.

for line in samfile:
    print(line)
    break

但顺序读取还不够灵活,我们有时需要随机读取(提示:sam不能随机读取),pysam的fetch方法提供了随机读取功能.

直接使用fetch会报错

ValueError: fetch called on bamfile without index

提示我们需要建立(.bai)索引

samtools index corrected.bam

fetch返回的是一个迭代器(iterator),可以迭代读取内容.

for read in samfile.fetch(‘chr6‘, 28478220, 28478222):
...     print(read)

fetch方法的API如下,chr6为参考序列,后面数字分别为读取的起始和终止位置.

fetch(self, reference=None, start=None, end=None, region=None, tid=None, until_eof=False, multiple_iterators=False)

3.读取vcf/bcf文件(pysam.VariantFile)

读取方法同上,只是使用的是VariantFile方法:

gvcf = "./MHC.unified.g.vcf.gz"
vcf_in = pysam.VariantFile(gvcf)

若想随机读取,仍然需要建立索引:

首先使用bgzip压缩vcf

bgzip -c MHC.g.vcf > MHC.g.vcf.gz

然后用bcftools建立索引

bcftools index -c MHC.g.vcf.gz

使用fetch读取

for rec in vcf_in.fetch(‘chr6‘, 28577796, 28577896):
...     print(rec)
...     break

4.创建并写入到新的bam或vcf文件

pysam的核心功能是可以随心所欲的读取数据,处理之后,写入到一个新建的bam或bcf文件里.

我们可以完全自定义一些内容,然后写入到一个新的bam文件里,如下:

header = { ‘HD‘: {‘VN‘: ‘1.0‘},
            ‘SQ‘: [{‘LN‘: 1575, ‘SN‘: ‘chr1‘},
                   {‘LN‘: 1584, ‘SN‘: ‘chr2‘}] }

with pysam.AlignmentFile(tmpfilename, "wb", header=header) as outf:
    a = pysam.AlignedSegment()
    a.query_name = "read_28833_29006_6945"
    a.query_sequence="AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG"
    a.flag = 99
    a.reference_id = 0
    a.reference_start = 32
    a.mapping_quality = 20
    a.cigar = ((0,10), (2,1), (0,25))
    a.next_reference_id = 0
    a.next_reference_start=199
    a.template_length=167
    a.query_qualities = pysam.qualitystring_to_array("<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<")
    a.tags = (("NM", 1),
              ("RG", "L1"))
    outf.write(a)

同理,我们也可以读取一个已有的bam文件,逐个修改以上的属性,然后存储到一个新的bam文件里.这里不再举例.

上面设置header可能有点麻烦,容易出错,但我们可以复制一个已有bam文件的header到一个新的bam文件里.

outf = pysam.AlignmentFile(path_out, "wb", template=samfile)

以上template参数指定了模板bam文件.

5. 关闭文件

outf.close()

总结:

pysam模块非常实用,有了pysam模块,我们就可以非常灵活的操纵bam/bcf文件,而不必依赖于samtools或bcftools. pysam可以随机读取bam/bcf文件,也可以将处理后的内容自定义输出到bam/bcf文件.

以上只介绍了pysam最常见的功能,更多pysam功能请参照:http://pysam.readthedocs.io/en/latest/index.html

时间: 2024-10-28 10:40:46

pysam - 多种格式基因组数据(sam/bam/vcf/bcf/cram/…)读写与处理模块(python)的相关文章

Edit Distance编辑距离(NM tag)- sam/bam格式解读进阶

sam格式很精炼,几乎包含了比对的所有信息,我们平常用到的信息很少,但特殊情况下,我们会用到一些较为生僻的信息,关于这些信息sam官方文档的介绍比较精简,直接看估计很难看懂. 今天要介绍的是如何通过bam文件统计比对的indel和mismatch信息 首先要介绍一个非常重要的概念--编辑距离 定义:从字符串a变到字符串b,所需要的最少的操作步骤(插入,删除,更改)为两个字符串之间的编辑距离. 这也是sam文档中对NM这个tag的定义. 编辑距离是对两个字符串相似度的度量(参见文章:Edit Di

sam/bam格式

一)Sam (Sequence Alignment/Map) ------------------------------------------------- 1) SAM 文件产生背景 随着Illumina/Solexa, AB/SOLiD and Roche/454测序技术不断的进步,各种比对工具产生,被用来高效的将reads比对到参考基因组.因为这些比对工具产生不同格式的文件,导致下游分析比较困难,因此一个通用的格式可以提供一个很好的接口用于链接比对与下游分析(组装,变异等,基因分型等)

SAM/BAM文件处理

当测序得到的fastq文件map到基因组之后,我们通常会得到一个sam或者bam为扩展名的文件.SAM的全称是sequence alignment/map format.而BAM就是SAM的二进制文件(B取自binary). 那么SAM文件的格式是什么样子的呢?如果你想真实地了解SAM文件,可以查看它的说明文档.SAM由头文件和map结果组成.头文件由一行行以@起始的注释构成.而map结果是类似下面的东西: HWI-ST1001:137:C12FPACXX:7:1115:14131:66670

多种格式多种字体版本的条形码控件USPS Postnet &amp; Intelligent Mail Barcode Font Package

IDAutomation的USPS Postnet & Intelligent Mail Barcode Font Package在六个不同的格式中包含了若干不同的字体版本.其中包括字体工具,宏和源代码,以帮助将字体集成到应用程序中 具体功能: 提供的字体工具- 提供的这些字体工具可用于协助应用程序的集成.这些字体工具可以为条形码字体自动地格式化开始,结束及校验字符.所含的可用字体工具如下: .NET 框架程序集DLL (100%托管代码) C++ 头文件 C# .NET类库 Visual Ba

Android Volley获取json格式的数据

为了让Android能够快速地访问网络和解析通用的数据格式Google专门推出了Volley库,用于Android系统的网络传输.volley库可以方便地获取远程服务器的图片.字符串.json对象和json对象数组等.当然,java本身也有获取json对象的方法,然而为了更好地适应移动互联网,google专门为其做了特殊的优化,因而应该尽可能地使用Volley库. Volley官方文档:https://developer.android.com/training/volley/index.htm

mybatis 处理数组类型及使用Json格式保存数据 JsonTypeHandler and ArrayTypeHandler

mybatis 处理数组类型及使用Json格式保存数据 JsonTypeHandler and ArrayTypeHandler mybatis 比 ibatis 改进了很多,特别是支持了注解,支持了plugin inteceptor,也给开发者带来了更多的灵活性,相比其他ORM,我还是挺喜欢mybatis的. 闲言碎语不要讲,今天研究了下mybatis的typeHandler: 先看这样一张表(postgresql) create table user ( id serial not null

导入CSV格式的数据

导入CSV格式的数据 (参见http://dev.mysql.com/doc/refman/5.6/en/load-data.html) 1.数据库表(st_pptn_r) CREATE TABLE st_pptn_r ( STCD varchar(8) DEFAULT NULL, TM datetime DEFAULT NULL, DRP decimal(5,1) DEFAULT NULL, INTV decimal(5,2) DEFAULT NULL, PDR decimal(5,2) DE

第二章 导入数据到SAS | 格式规范数据读取

目录 2.1 导入数据的方法 2.2 利用导入向导读入 2.3 格式规范数据读取 2.3.1 指定原始数据位置(infile) 2.3.2 读取空格分隔原始数据(列表输入) 2.3.3 读取按列排列原始数据(列输入) 2.3.4 读取非标准格式的原始数据(格式化输入) 2.3.5 混合的输入样式(列表输入+列输入+格式输入) 2.1 导入数据的方法 将数据导入SAS的方法有很多,但可以归纳为四个基本类别,其中方法2.3是需要掌握的重点. 直接将数据输入SAS数据集 通过VIEWTABLE窗口(打

NFC技术:读写非NDEF格式的数据

1 //向nfc标签读写MifareUltraligh格式的数据 2 public class MainActivity extends Activity { 3 private CheckBox mwriteData; 4 private NfcAdapter mNfcAdapter; 5 private PendingIntent mPendingIntent; 6 7 @Override 8 protected void onCreate(Bundle savedInstanceState