macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

samのフィルタリングツール SAMsift

 

 

SAMsiftはKarel BřindaさんがGithubで公開されている、samを様々な条件でフィルタリングできるツール。

 

インストール

mac os10.13のPython 3.6.2 :: Anaconda 3-5.0.0 でテストした。

本体 GIthub

#Anaconda環境
conda install -c bioconda samsift

#Anaconda使ってないならpipで導入
pip install samsift

> samsift -h

$ samsift -h

 

Program: samsift (advanced filtering and tagging of SAM/BAM alignments using Python expressions)

Version: 0.2.5

Author:  Karel Brinda <kbrinda@hsph.harvard.edu>

 

Usage:   samsift [-i FILE] [-o FILE] [-f PY_EXPR] [-c PY_CODE] [-m STR] 

                    [-0 PY_CODE] [-d PY_EXPR] [-t PY_EXPR]

 

Basic options:

  -h, --help            show this help message and exit

  -v, --version         show program's version number and exit

  -i FILE               input SAM/BAM file [-]

  -o FILE               output SAM/BAM file [-]

  -f PY_EXPR            filtering expression [True]

  -c PY_CODE            code to be executed (e.g., assigning new tags) [None]

  -m STR                mode: strict (stop on first error)

                              nonstop-keep (keep alignments causing errors)

                              nonstop-remove (remove alignments causing errors) [strict]

 

Advanced options:

  -0 PY_CODE            initialization [None]

  -d PY_EXPR            debugging expression to print [None]

  -t PY_EXPR            debugging trigger [True]

 

 

ラン

ここではSAMsiftのGithubレポジトリ:test.sam(samsift/tests/tests/)を使っている。

> head -n 5 samsift/tests/tests/test.sam 

$ head -n 5 samsift/tests/test.sam 

@SQ SN:gi|561108321|ref|NC_018143.2| LN:4411709

r00001 0 gi|561108321|ref|NC_018143.2| 2798553 60 100M * 0 0 AAAGGAATGCATCACCAATCTCGTCGGCGAGGCTTCGTGCCTCCTCGCCCGGCGCGCGGGTCGCGCCCGGGTCACTCTCGCCGGCGAACCCGACATAGGC 2/454483235524234201022/1252422523223533243212222421513222236220215432131222242315235121223130142647 NM:i:2 MD:Z:51C23C24 AS:i:90 XS:i:0

r00002 0 gi|561108321|ref|NC_018143.2| 445039 60 100M * 0 0 TCATCGGCAGATGCCCCCCGAGCGCAGCTTCCACGCGGCGCCGCGACGCCGTGTGCTGGTTCGTCACCGCCCGACCGACGCGGGCCCAGTGGTCGAGCTG 11/242522323122622254122142322321121236232321205072426233223223204322022022633343113212322232-462422 NM:i:2 MD:Z:17G45A36 AS:i:90 XS:i:0

r00003 0 gi|561108321|ref|NC_018143.2| 2014256 60 100M * 0 0 GGCAGCCCCGATGGGGGACCAGATCGGCCCCTGGGCCGCCTCGCTCACCGCGAGGCTCGCGAGGAGTCCGACCAGCGCGGCGCCCACGGCCACAATCACG 2220242203322212123344242421/2-0042322131152322232426231022262222244315/353/227312462213575444122223 NM:i:2 MD:Z:19G50G29 AS:i:90 XS:i:0

r00004 0 gi|561108321|ref|NC_018143.2| 3675128 60 100M * 0 0 GTTGGGCGGGGTGGAATTCATCCATTTCATTCAGTGCCCGTTGCGAATCCCCAAGCTACCCCGACGGCGACCAGAGGATGTCGATGGGGACGGCGGCGAG 28542253205043252212122151202.212332342342273312355423123122235524002520/2421/242130/623202/22221632 NM:i:0 MD:Z:100 AS:i:100 XS:i:0

 

1、アライメントクオリティ(AS)が94以上を抽出し、bam出力

filtered.bam
samsift -i tests/test.bam -f 'AS>94' -o filtered.bam 

ASはSAMの最後のOptional tags(Optional field)にある。

f:id:kazumaxneo:20180915174705p:plain

ASとMAPQの違いはSEQanswersの回答を参照(リンク)。

 

2、アンマップのリードだけ出力

samsift -i tests/test.bam -f 'FLAG & 0x04' -o filtered.bam 

f:id:kazumaxneo:20180915174309p:plain

SAMフォーマットPDFより抜粋(SAM format  v1.6 PDF)。FLAGの0x04(10進数の4)はアンマップを意味する。

 

3、"ACCAGAGGAT"を含むリードだけ抽出

samsift -i tests/test.bam -f 'SEQ.find("ACCAGAGGAT")!=-1'

 

4、AとTのみ含むリードだけ抽出(それ以外の文字が配列中に見つかれば除かれる)

samsift -i tests/test.bam -f 're.match(r"^[AT]*$", SEQ)'

pythonのre.モジュールでの正規表現に従う。

 

5、25% ランダムサンプリング(randomモジュールを使っている)

samsift -i tests/test.bam -f 'random.random()<0.25'

 

6、乱数生成器を固定して25% ランダムサンプリング

samsift -i tests/test.bam -f 'random.random()<0.25' -0 'random.seed(42)'

ペアエンドリードのそれぞれに対して同じ値を使ったりする。

 

リストに名前があるリードだけ出力

samsift -i tests/test.bam -0 'q=open("tests/qnames.txt").read().splitlines()' -f 'QNAME in q'

 

 

タグ追加。リード長としてlnタグを追加し、平均ベースクオリティとしてabタグを追加する。

samsift -i tests/test.bam -c 'ln=len(SEQ);ab=1.0*sum(QUALa)/ln'

オリジナル NM:i:2 MD:Z:51C23C24 AS:i:90 XS:i:0

 ↓

出力    NM:i:2 MD:Z:51C23C24 AS:i:90 XS:i:0 ln:i:100 ab:f:50.62

 

 

全リードをアンマップに書き換える。

samsift -i tests/test.bam -c 'a.flag|=0x4;a.reference_start=-1;a.cigarstring="";a.reference_id=-1;a.mapping_quality=0'

 

引用

GitHub - karel-brinda/samsift: SAMsift: advanced filtering and tagging of SAM/BAM alignments using Python expressions.