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
１、アライメントクオリティ（AS）が94以上を抽出し、bam出力
filtered.bam
samsift -i tests/test.bam -f 'AS>94' -o filtered.bam
ASはSAMの最後のOptional tags（Optional field）にある。
ASとMAPQの違いはSEQanswersの回答を参照（リンク）。
２、アンマップのリードだけ出力
samsift -i tests/test.bam -f 'FLAG & 0x04' -o filtered.bam
SAMフォーマットPDFより抜粋（SAM format v1.6 PDF）。FLAGの0x04（10進数の4）はアンマップを意味する。
３、"ACCAGAGGAT"を含むリードだけ抽出
samsift -i tests/test.bam -f 'SEQ.find("ACCAGAGGAT")!=-1'
４、AとTのみ含むリードだけ抽出（それ以外の文字が配列中に見つかれば除かれる）
samsift -i tests/test.bam -f 're.match(r"^[AT]*$", SEQ)'
５、25% ランダムサンプリング（randomモジュールを使っている）
samsift -i tests/test.bam -f 'random.random()<0.25'
６、乱数生成器を固定して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'
引用
