生信分析 - Pysam应用

罗大黑学生信 2021-07-26

2177

Pysam应用

Pysam包是一个处理基因组数据的python模块，它打包了htslib-1.3、samtools-1.3 和 bcftools-1.3的核心功能，能在编程时非常灵活的处理bam和bcf文档，实现python处理基因组数据的无缝衔接，而不用在python进程内部调用samtools、bcftools等软件。

函数

Pysam的函数有很多，主要的读取函数有：

AlignmentFile：读取BAM/CRAM/SAM文件；
VariantFile：读取变异数据（VCF或者BCF）；
TabixFile：读取由tabix索引的文件；
FastaFile：读取fasta序列文件；
FastqFile：读取fastq测序序列文件

FastaFile：读取fasta序列文件

import pysam


## 1 构建FastaFile对象,调用后返回对象，dir(对象)查看其内置方法
refGenome="Homo_sapiens_assembly19.fasta"
refseq = pysam.FastaFile(refGenome)
>>> dir(refseq)
['__class__', '__contains__', '__delattr__', '__dir__', '__doc__', '__enter__', '__eq__', '__exit__', '__format__', '__ge__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__init__', '__le__', '__len__', '__lt__', '__ne__', '__new__', '__pyx_vtable__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', '_open', 'close', 'closed', 'fetch', 'filename', 'get_reference_length', 'is_open', 'lengths', 'nreferences', 'references']




## 2 用fetch函数随机读取序列， 查看其使用方法help(object.fetch())
### 提取整条序列
refseq.fetch("chr1")
>>> 'AGCTACTGCTAGCATACGATCTAACGTAGCTCTCTCAGCCGATATTCGCGAT'


### 提取具体位点的碱基，比如上述从左数起第5,6个碱基'AC'
####  Python风格半开区间：半开区间碱基位置编号从0开始。refseq.fetch("chr1", 4, 6)
>>>'AC'


#### Samtools风格闭区间：碱基位置编号从1开始
refseq.fetch(region="chr1:5-6")
>>>'AC'

AlignmentFile：读取BAM/CRAM/SAM文件

## 1 构建AlignmentFile对象,调用后返回对象，dir(对象)查看其内置方法
samfile = pysam.AlignmentFile("example.bam", "rb") 
>>> dir(samfile)
['__class__', '__delattr__', '__dir__', '__doc__', '__enter__', '__eq__', '__exit__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__iter__', '__le__', '__lt__', '__ne__', '__new__', '__next__', '__pyx_vtable__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', '_open', 'add_hts_options', 'category', 'check_index', 'check_truncation', 'close', 'closed', 'compression', 'count', 'count_coverage', 'description', 'duplicate_filehandle', 'fetch', 'filename', 'find_introns', 'find_introns_slow', 'format', 'get_index_statistics', 'get_reference_length', 'get_reference_name', 'get_tid', 'getrname', 'gettid', 'has_index', 'head', 'header', 'index_filename', 'is_bam', 'is_bcf', 'is_closed', 'is_cram', 'is_open', 'is_read', 'is_remote', 'is_sam', 'is_stream', 'is_valid_reference_name', 'is_valid_tid', 'is_vcf', 'is_write', 'lengths', 'mapped', 'mate', 'mode', 'nocoordinate', 'nreferences', 'parse_region', 'pileup', 'reference_filename', 'references', 'reset', 'seek', 'tell', 'text', 'threads', 'unmapped', 'version', 'write']




## 2 用fetch函数随机读取序列， 查看其使用方法help(object.fetch()),  dir(object)查看内置方法
### 提取出比对到目标区域内的全部reads。返回的是一个迭代器，可以通过for循环或者next函数从中取出reads，reads是用AlignedSegment对象表示，可以通过该对象的内置方法再对这条reads进行一些查询操作。allreads = samfile.fetch("chr17")
>>> next(allreads).to_dict()
{'seq': 'GCTTTGGAAGAACCAATCAAGAAAGGATCCTGGGTGTTTGTATTTGCAGTCAAGTCTTCCAATTCACTGCACTGTGAAGAAAACAAGCTAGCAGAACATTTTGTTTCCTCACTAAGGTGATGTTCCTG', 'ref_name': 'chr17', 'ref_pos': '41243582', 'next_ref_pos': '41243648', 'map_quality': '60', 'flag': '99', 'name': 'NB500965:27:H7JT7AFXX:1:21207:11689:5709:chr17-1-41243561-GAAGTGCATGGA:GAAGTGCATGGA', 'length': '216', 'qual': '//6/EEEEEEEA/EEEEE<EEEEEAEEEAEEEEEEE<EEEEEAEEEEEE/EEAEEEEAEE/<AEEEEEAEEEEEE<<AAEEEE6EEEEEEAAAAEEAAA/A<A6<A6<A<AEA<E<<6<AAA6AEEE/', 'next_ref_name': '=', 'tags': ['MD:Z:126', 'RG:Z:1', 'NM:i:0', 'AS:i:126', 'XS:i:0', 'Bc:Z:chr17-1-41243561-GAAGTGCATGGA', 'bc:Z:GAAGTGCATGGA', 'tr:Z:chr17-1-41243796-41'], 'cigar': '2S126M'}




## 3. 用pileup函数提取具体位点的pileup信息， 类似：samtools mpileup
pileup(self, contig=None, start=None, stop=None, region=None, reference=None, end=None, **kwargs)
类似于samtools的pileup操作，返回一个迭代器，每次迭代返回一个位点的PileupColumn对象。### 创建pileup对象
### 调用pileup 函数
### 调用内置方法


sampileup = samfile.pileup(region = chrom + ':' + pos + '-' + pos, truncate=True, max_depth=1000000, stepper='nofilter')
>>> sampileup = samfile.pileup(region = "chr17:41243582-41243582", truncate=True, max_depth=1000000, stepper='nofilter')
>>> for read in sampileup:
...     for pileupRead in read.pileups:
...             dir(pileupRead)
...             print(pileupRead)
...             break
... 
['__class__', '__delattr__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__le__', '__lt__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', 'alignment', 'indel', 'is_del', 'is_head', 'is_refskip', 'is_tail', 'level', 'query_position', 'query_position_or_next']
NB500965:27:H7JT7AFXX:1:11104:2773:7813:chr17-1-41243571-TAGACATAGAGC:TAGACATAGAGC  99  17  41243572  60  128M  17  41243647  128  CCTCATTTGTTTGGAAGAACCAATCAAGAAAGGATCCTGGGTGTTTGTATTTGCAGTCAAGTCTTCCAATTCACTGCACTGTGAAGAAAACAAGCTAGCAGAACATTTTGTTTCCTCACTAAGGTGAT  array('B', [21, 27, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 32, 36, 36, 36, 36, 14, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 14, 36, 32, 36, 36, 32, 36, 36, 36, 36, 36, 36, 36, 36, 32, 36, 36, 36, 36, 36, 36, 36, 14, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 32, 32, 36, 32, 32, 36, 36, 32, 36, 32, 36, 32, 36, 36, 36, 36, 36, 36, 36, 36, 32, 36, 36, 36, 32, 36, 36, 21, 27, 36, 36, 27, 36, 36, 32, 27, 32, 32, 14, 32, 36, 36, 27, 32, 36, 32, 32, 32, 32, 32, 36, 32, 32, 32, 36, 32, 27])  [('MD', '128'), ('RG', '1'), ('NM', 0), ('AS', 128), ('XS', 0), ('Bc', 'chr17-1-41243571-TAGACATAGAGC'), ('bc', 'TAGACATAGAGC'), ('tr', 'chr17-1-41243796-41')]  9  0  0  0  0  0  0


>>> sampileup = samfile.pileup(region = "chr17:41243582-41243582", truncate=True, max_depth=1000000, stepper='nofilter')
>>> for read in sampileup:
...     for pileupRead in read.pileups:
...             dir(pileupRead.alignment)
...             break
... 
['__class__', '__copy__', '__deepcopy__', '__delattr__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__le__', '__lt__', '__ne__', '__new__', '__pyx_vtable__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', 'aend', 'alen', 'aligned_pairs', 'bin', 'blocks', 'cigar', 'cigarstring', 'cigartuples', 'compare', 'flag', 'from_dict', 'fromstring', 'get_aligned_pairs', 'get_blocks', 'get_cigar_stats', 'get_forward_qualities', 'get_forward_sequence', 'get_overlap', 'get_reference_positions', 'get_reference_sequence', 'get_tag', 'get_tags', 'has_tag', 'header', 'infer_query_length', 'infer_read_length', 'inferred_length', 'is_duplicate', 'is_paired', 'is_proper_pair', 'is_qcfail', 'is_read1', 'is_read2', 'is_reverse', 'is_secondary', 'is_supplementary', 'is_unmapped', 'isize', 'mapping_quality', 'mapq', 'mate_is_reverse', 'mate_is_unmapped', 'mpos', 'mrnm', 'next_reference_id', 'next_reference_name', 'next_reference_start', 'opt', 'overlap', 'pnext', 'pos', 'positions', 'qend', 'qlen', 'qname', 'qqual', 'qstart', 'qual', 'query', 'query_alignment_end', 'query_alignment_length', 'query_alignment_qualities', 'query_alignment_sequence', 'query_alignment_start', 'query_length', 'query_name', 'query_qualities', 'query_sequence', 'reference_end', 'reference_id', 'reference_length', 'reference_name', 'reference_start', 'rlen', 'rname', 'rnext', 'seq', 'setTag', 'set_tag', 'set_tags', 'tags', 'template_length', 'tid', 'tlen', 'to_dict', 'to_string', 'tostring']

# read ID

qname = pileupRead.alignment.query_name

# mapping quality

bq = pileupRead.alignment.mapping_quality

# get NM tag

# 以元祖的形式返回tag

> pileupRead.alignment.tags
如：[('NM', 12), ('MD', '106'), ('MC', '73M55S'), ('AS', 88), ('XS', 19)]

# count number of INDELs in the read sequence

# 统计indel 的个数长度

cigar = pileupRead.alignment.cigar
for (op, value) in cigar:
    # 1 for insertion, 2 for deletion
    if op == 1 or op == 2:
        nIndel += value

# Number of mismatches except INDEL

# 判断碱基错配的个数

mismatch = max(0, NM - nIndel)

# paired read

# 判断read 是R1 or R2

if pileupRead.alignment.is_read1:
    pairRead = 'R1'
if pileupRead.alignment.is_read2:
    pairRead = 'R2'

# +/- strand

# 判断read 是'Forward' or 'Reverse'

'Reverse' if pileupRead.alignment.is_reverse else 'Forward'

# indel

# 后面的碱基是否处在indel上，若在insertion上则返回正数，数值为insertion长度，若在deltion上则为负数，数值为deltion长度。

pileupRead.indel

# query_sequence

# 获取比对的read的序列

pileupRead.alignment.query_sequence

# query_position

# 该碱基位点在query上的位置，0代表第一个。若处在deltion上或者是refskip则返回None。

pileupRead.query_position

# query_qualities

# 获取某位置碱基的质量

pileupRead.alignment.query_qualities[pileupRead.query_position]

cigar记录的元组格式

下面的例子就代表了128M（0 代表 M）

>>> pileupRead.alignment.cigar
[(0, 128)]

参考:

https://www.dazhuanlan.com/2019/12/08/5dec78f1a79fb/

https://zhuanlan.zhihu.com/p/297858072

https://www.dazhuanlan.com/2019/12/08/5dec78f1a79fb/

https://vimsky.com/examples/detail/python-module-pysam.html

https://m.jb51.net/article/148115.htm

http://blog.sina.cn/dpool/blog/s/blog_7948610e01013xvn.html

数据库

文章转载自罗大黑学生信，如果涉嫌侵权，请发送邮件至：contact@modb.pro进行举报，并提供相关证据，一经查实，墨天轮将立刻删除相关内容。

生信分析 - Pysam应用

评论