Pythonでバイオインフォマティクス - BAM/SAMファイルを扱う

f:id:kimoppy126:20201025155625p:plain

インストール手順は大抵conda かpipで可能なので、これらをいちいち紹介しても仕方ない気がしてきました。今回からインストール部分は省いてご紹介します。

前回の記事ではアライメント情報を格納するためのフォーマット、BAM/SAMファイルをご紹介しました。 www.kimoton.com

今回はそんなBAM/SAMファイルをPythonで扱う例を見ていきます。
BAM/SAMファイルを扱うためのツールといえば、Samtoolsです。Samtoolsは内部でHTSlibというC言語製のAPIを使用しています。
pysamはこのHTSlibのラッパーであり、バイオインフォマティクス関連ファイル(SAM / BAM / VCF / BCF / BED / GFF / GTF / FASTA / FASTQ)の読み取りと書き込み、およびSamtoolsとBCFtoolsパッケージのCLI機能へのpython APIを提供しています。

BAMファイルをとってくる

早速、シーケンサによる標準的な出力フォーマットであるFASTQファイルをとってきましょう。今回は1000人ゲノムプロジェクトのデータを用います。
1000人ゲノムプロジェクトのデータはFTPで公開されており、 Data portalから検索することができます。

BAM/SAMファイルのフォーマットって何?って方は下記の記事を参考にして下さい。 https://www.kimoton.com/entry/sam

.bamとそのインデックスファイルである.bam.baiをとってきます。 インデックスファイルは

$ wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/exome_alignment/NA18489.chrom20.ILLUMINA.bwa.YRI.exome.20121211.bam
$ wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/exome_alignment/NA18489.chrom20.ILLUMINA.bwa.YRI.exome.20121211.bam.bai

BAM/SAMファイルの確認

BAMファイルを扱う際には、Samtoolsというツールを頻繁に使用します。

samtools view -hを使うと、バイナリファイルであるBAMファイルの中身を見ることができます(-h はヘッダーを含めるためのオプション)。

# samtools view でBAMファイルの中身を確認
$ samtools view NA18489.chrom20.ILLUMINA.bwa.YRI.exome.20121211.bam | head
SRR100025.13877927      163     20      59995   60      76M     =       60021   101     GACTCAGATCCAGAGGTGGAAGAGGAAGGAAGCTTGGAACCCTATAGAGTTGCTGAGTGCCAGGACCAGATCCTGG    =?DCE@GA@FGBHBIG@HGBEH@HHBFIIDFIGGDIICDGGHF@B=FAG7CGHFG6A<C,A:=C;A=,?<7<<=@C    X0:i:1  X1:i:0  MD:Z:0N0N0N0N0N0N70     RG:Z:SRR100025  AM:i:37 NM:i:6  SM:i:37 XN:i:6  MQ:i:60 XT:A:U  BQ:Z:GIMKLEH@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@A
SRR100025.62130839      99      20      59997   60      52M24S  =       60229   295     CTCAGATCCAGAGGTGGAAGAGGAAGGAAGCTTGGAACCCTATAGAGTTGCTGAGTGCCAGGACCAGATACTGGGC    BCEBHCBGHCIDIIAIGBDGBHICFHE?EFCDCIFCG=IIGABAECF9CD@#########################    X0:i:1  X1:i:0  XC:i:52 MD:Z:0N0N0N0N48 RG:Z:SRR100025  AM:i:37 NM:i:4  SM:i:37 XN:i:4  MQ:i:60 XT:A:U  BQ:Z:JKLGI@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SRR100025.58882446      121     20      60006   37      76M     =       60006   0       AGAGGTGGAAGAGGAAGGAAGCTTGGAACCCTATAGAGTTGCTGAGTGCCAGGACCAGATCCTGGCCCTAAACAGG    B@D<;<;<5E/E<<AFE+DGEFEAE<CAEDGACBGIGHGDJJDHHHBHGIFIEBGJFHCDIJCHFIIFACCAIDEB    X0:i:1  X1:i:0  MD:Z:76 RG:Z:SRR100025  AM:i:0  NM:i:0  SM:i:37 XT:A:U  BQ:Z:@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@A
SRR100025.58882446      181     20      60006   0       64M12S  =       60006   0       CCGGGCACCCAGCAACTCTATAGGGTTCCAAGCTTCCTTCCTCTTCCACCTCTGGATCTGAGTCACCTCTCTATTA    #############@>DAF2A?A;<>A*E9CDHH@>GFB>DF?>D?FI9FHAI7HDACIBGFGCG?GHAGAF??C==    XC:i:64 RG:Z:SRR100025
SRR100025.76297965      99      20      60007   60      55M21S  =       60023   91      GAGGTGGAAGAGGAAGGAAGCTTGGAACCCTATAGAGTTGCTGAGTGCCAGGACCAGATCCTGGCCCTAAACAGGG    @>CC@HGABCBH@BFIHBECBFCIGBABFFGAB>GBFBCJB?H>ECIE;@GF?@######################    X0:i:1  X1:i:0  XC:i:55 MD:Z:55 RG:Z:SRR100025  AM:i:37 NM:i:0  SM:i:37 MQ:i:60 XT:A:U  BQ:Z:@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SRR100025.13877927      83      20      60021   60      76M     =       59995   -101    AGGAAGCTTGGAACCCTATAGAGTTGCTGAGTGCCAGGACCAGATCCTGGCCCTAAACAGGTGGTAAGGAAGGAGA    B;<@D?DBA:AA@AFBCCCHHGIECHFDIHICGGFFIIBJKHIADIJDIIHFJBCDCHFIHCIHBDGHGCFGGDDB    X0:i:1  X1:i:0  MD:Z:76 RG:Z:SRR100025  AM:i:37 NM:i:0  SM:i:37 MQ:i:60 XT:A:U  BQ:Z:@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@A
SRR100025.76297965      147     20      60023   60      76M     =       60007   -91     GAAGCTTGGAACCCTATAGAGTTGCTGAGTGCCAGGACCAGATCCTGGCCCTAAACAGGTGGTAAGGAAGGAGAGA    C?BDFBAGED@IGFBBBGHEHDCHJCFFEC=GHFBGBIJGHBDIICE;FDD@CCBJGGGBHG@CFGFBEFEDDCC=    X0:i:1  X1:i:0  MD:Z:76 RG:Z:SRR100025  AM:i:37 NM:i:0  SM:i:37 MQ:i:60 XT:A:U  BQ:Z:@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SRR100025.29222641      163     20      60057   60      76M     =       60075   77      GGACCAGATCCTGGCCCTAAACAGGTGGTAAGGAAGGAGAGAGTGAAGGAACTGCCAGGTGACACACTCCCACCAT    =D@DE@DA@FGFIHHHHFAFFGBIIAIG??EFFB>HH>HCHCH3ICCECB>GFGIH?HH@JCB3,1/@)6A;><5:    X0:i:1  X1:i:0  MD:Z:76 RG:Z:SRR100025  AM:i:37 NM:i:0  SM:i:37 MQ:i:60 XT:A:U  BQ:Z:@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SRR100025.29222641      83      20      60075   60      16S60M  =       60057   -77     ACCAGATCCTGCCCCTAAACAGGTGGTAAGGAAGGAGAGAGTGAAGGAACTGCCAGGTGACACACTCCCACCATGG    #################C?HFG=CAF9CEE4CGFEF6FEF:=@DFFG@BHCCGHFGC;BAH>H;E+FFH@FG??B@    X0:i:1  X1:i:0  XC:i:60 MD:Z:60 RG:Z:SRR100025  AM:i:37 NM:i:0  SM:i:37 MQ:i:60 XT:A:U  BQ:Z:@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SRR100025.55446050      81      20      60077   37      3S73M   hs37d5  17693622        0       TAAACAGGTGGTAAGGAAGGAGAGAGTGAAGGAACTGCCAGGTGACACACTCCCACCATGGACCTCTGGGATCCTA    ####DDB<;AA>BBE9CBG=FCGGGHDCCFAEA@HCIBJGFGCGBJBI?I@GGIAAJBCIH?IIDICHHGBCHG>B    X0:i:1  X1:i:0  XC:i:73 MD:Z:73 RG:Z:SRR100025  AM:i:37 NM:i:0  SM:i:37 XT:A:U  BQ:Z:@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# オプション -h を付けると、ヘッダーを含めたファイルの中身を確認できる 
$ samtools view -h NA18489.chrom20.ILLUMINA.bwa.YRI.exome.20121211.bam | head
@HD     VN:1.0  SO:coordinate
@SQ     SN:1    LN:249250621    M5:1b22b98cdeb4a9304cb5d48026a85128     UR:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
@SQ     SN:2    LN:243199373    M5:a0d9851da00400dec1098a9255ac712e     UR:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
@SQ     SN:3    LN:198022430    M5:fdfd811849cc2fadebc929bb925902e5     UR:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
@SQ     SN:4    LN:191154276    M5:23dccd106897542ad87d2765d28a19a1     UR:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
@SQ     SN:5    LN:180915260    M5:0740173db9ffd264d728f32784845cd7     UR:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
@SQ     SN:6    LN:171115067    M5:1d3a93a248d92a729ee764823acbbc6b     UR:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
@SQ     SN:7    LN:159138663    M5:618366e953d6aaad97dbe4777c29375e     UR:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
@SQ     SN:8    LN:146364022    M5:96f514a9929e410c6651697bded59aec     UR:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
@SQ     SN:9    LN:141213431    M5:3e273117f15e0a400f01055d9f393768     UR:ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human

pysam でBAM/SAMファイルを扱う

pysamを使うと、C言語製のsamtoolsをPython経由で使用することができます。
インストールはconda 若しくはpipから可能です。こちらを参考にしてください。

pysam でBAM/SAMファイルを読み込むにはAlignmentFileクラスを使用します。

bam = pysam.AlignmentFile('NA18489.chrom20.ILLUMINA.bwa.YRI.exome.20121211.bam', 'rb')

生成されるAlignmentFileオブジェクトはアライメントセクションのジェネレータになっています。

> print(next(bam))
SRR100025.13877927 163    19 59994  60 76M 19 60020  76 GACTCAGATCCAGAGGTGGAAGAGGAAGGAAGCTTGGAACCCTATAGAGTTGCTGAGTGCCAGGACCAGATCCTGG    array('B', [28, 30, 35, 34, 36, 31, 38, 32, 31, 37, 38, 33, 39, 33, 40, 38, 31, 39, 38, 33, 36, 39, 31, 39, 39, 33, 37, 40, 40, 35, 37, 40, 38, 38, 35, 40, 40, 34, 35, 38, 38, 39, 37, 31, 33, 28, 37, 32, 38, 22, 34, 38, 39, 37, 38, 21, 32, 27, 34, 11, 32, 25, 28, 34, 26, 32, 28, 11, 30, 27, 22, 27, 27, 28, 31, 34])   [('X0', 1), ('X1', 0), ('MD', '0N0N0N0N0N0N70'), ('RG', 'SRR100025'), ('AM', 37), ('NM', 6), ('SM', 37), ('XN', 6), ('MQ', 60), ('XT', 'U'), ('BQ', 'GIMKLEH@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@A')]
> print(next(bam))
SRR100025.62130839 99 19 59996  60 52M24S  19 60228  52 CTCAGATCCAGAGGTGGAAGAGGAAGGAAGCTTGGAACCCTATAGAGTTGCTGAGTGCCAGGACCAGATACTGGGC    array('B', [33, 34, 36, 33, 39, 34, 33, 38, 39, 34, 40, 35, 40, 40, 32, 40, 38, 33, 35, 38, 33, 39, 40, 34, 37, 39, 36, 30, 36, 37, 34, 35, 34, 40, 37, 34, 38, 28, 40, 40, 38, 32, 33, 32, 36, 34, 37, 24, 34, 35, 31, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])    [('X0', 1), ('X1', 0), ('XC', 52), ('MD', '0N0N0N0N48'), ('RG', 'SRR100025'), ('AM', 37), ('NM', 4), ('SM', 37), ('XN', 4), ('MQ', 60), ('XT', 'U'), ('BQ', 'JKLGI@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@')]

また、ヘッダーは.header属性に辞書に近い形式で格納されています。

> headers = bam.header
> headers.items
odict_items([('HD', {'VN': '1.0', 'SO': 'coordinate'}), ('SQ', [{'SN': '1', 'LN': 249250621, 'M5': '1b22b98cdeb4a9304cb5d48026a85128', 'UR': 'ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human'}, {'SN': '2', 'LN': 243199373, 'M5': 'a0d9851da00400dec1098a9255ac712e', 'UR': 'ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human'},

ヘッダーの中身を確認

headers = bam.header
for tag, records in headers.items():
    print (tag)
    for i, record in enumerate(records):
        if type(record) == dict:
            print('\t%d' % (i + 1))
            for field, value in record.items():
                print('\t\t%s\t%s' % (field, value))
        else:
            print('\t\t%s' % record)
HD
        VN
        SO
SQ
    1
        SN  1
        LN  249250621
        M5  1b22b98cdeb4a9304cb5d48026a85128
        UR  ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
    2
        SN  2
        LN  243199373
        M5  a0d9851da00400dec1098a9255ac712e
        UR  ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
.
.
.
    86
        SN  hs37d5
        LN  35477943
        M5  5b6a4b3a81a2d3c134b7d14bf6ad39f1
        UR  ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz        AS:NCBI37       SP:Human
RG
    1
        ID  SRR100025
        LB  Solexa-51039
        SM  NA18489
        PI  220
        CN  BI
        PL  ILLUMINA
        DS  SRP004074
PG
    1
        ID  bwa_index
        PN  bwa
        VN  0.5.9-r16
        CL  bwa index -a bwtsw $reference_fasta
    2
        ID  bwa_aln_fastq
        PN  bwa
        PP  bwa_index
        VN  0.5.9-r16
        CL  bwa aln -q 15 -f $sai_file $reference_fasta $fastq_file
.
.
.
    11
        ID  bam_mark_duplicates
        PN  picard
        PP  bam_merge
        VN  1.53
        CL  java $jvm_args -jar MarkDuplicates.jar INPUT=$bam_file OUTPUT=$markdup_bam_file ASSUME_SORTED=TRUE METRICS_FILE=/dev/null VALIDATION_STRINGENCY=SILENT
    12
        ID  bam_merge.1
        PN  picard
        PP  bam_mark_duplicates
        VN  1.53
        CL  java $jvm_args -jar MergeSamFiles.jar INPUT=$bam_file(s) OUTPUT=$merged_bam VALIDATION_STRINGENCY=SILENT
CO
        $known_indels_file(s) = ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_mapping_resources/ALL.wgs.indels_mills_devine_hg19_leftAligned_collapsed_double_hit.indels.sites.vcf.gz
        $known_indels_file(s) .= ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_mapping_resources/ALL.wgs.low_coverage_vqsr.20101123.indels.sites.vcf.gz
        $known_sites_file(s) = ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_mapping_resources/ALL.wgs.dbsnp.build135.snps.sites.vcf.gz

レコードの情報を確認

あるレコードについて、リファレンスID、参照しているリファレンスの名称、リファレンスの開始・終了位置を調べます。pysam で格納される開始位置・終了位置はSAMフォーマットと異なり0始まりとなっていることに注意してください(SAMファイル内の値 -1)。

> print(rec.query_name, rec.reference_id, bam.getrname(rec.reference_id), rec.reference_start, rec.reference_end)
SRR100025.76297965 19 20 60006 60061

60061 - 60006 = 55塩基がマッピングされていることが分かります。
続いてCIGAR文字列を確認します。

> print(rec.cigarstring)
55M21S

CIGAR文字列からは、残る21塩基はクリップ(soft clip)されていることが分かります。

続いて、アライメント開始位置・アライメント終了位置、アライメントされた長さを確認します。

> print(rec.query_alignment_start, rec.query_alignment_end, rec.query_alignment_length)
0 55 55

CIGAR文字列からわかるように、開始位置は0、終了位置は55となっていることが確認できます。

更に、各塩基のPhred Qualityも格納されています。

# クエリ配列のPhred Qualityスコアを取得
> print(rec.query_qualities)
array('B', [31, 29, 34, 34, 31, 39, 38, 32, 33, 34, 33, 39, 31, 33, 37, 40, 39, 33, 36, 34, 33, 37, 34, 40, 38, 33, 32, 33, 37, 37, 38, 32, 33, 29, 38, 33, 37, 33, 34, 41, 33, 30, 39, 29, 36, 34, 40, 36, 26, 31, 38, 37, 30, 31, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])

# アライメントされた部分だけver
> print(rec.query_alignment_qualities)
array('B', [31, 29, 34, 34, 31, 39, 38, 32, 33, 34, 33, 39, 31, 33, 37, 40, 39, 33, 36, 34, 33, 37, 34, 40, 38, 33, 32, 33, 37, 37, 38, 32, 33, 29, 38, 33, 37, 33, 34, 41, 33, 30, 39, 29, 36, 34, 40, 36, 26, 31, 38, 37, 30, 31, 2])

# クエリ配列自体を取得
> print(rec.query_sequence)
GAGGTGGAAGAGGAAGGAAGCTTGGAACCCTATAGAGTTGCTGAGTGCCAGGACCAGATCCTGGCCCTAAACAGGG

アラインされた領域の可視化

リードの各位置についてアラインされたカウントをプロットします。

counts = [0] * 76
for n, rec in enumerate(bam.fetch('20', 0, 10000000)):
    for i in range(rec.query_alignment_start, rec.query_alignment_end):
        counts[i] += 1
freqs = [x / (n + 1.) for x in counts]
fig, ax = plt.subplots(figsize=(16,9))
ax.plot(range(1, 77), freqs)

プロットから、中央に凹みがあることが分かります。
ペアエンド法でシーケンスした結果のため、このような結果になると推測されます。因みにペアエンド法は両端から読むってことです(下記参照)。

参考:http://www.clst.riken.jp/files/6114/3338/5221/10_Principle_for_Illumina_sequencer.pdf

Phredスコアの分布を可視化

各ポジションのPhreadスコアについて、最大、95%タイル、50%タイル、5%タイルを取得し、積上げ折れ線グラフであらわします。

# 各ポジションのPhread スコアについて、最大、95%タイル、50%タイル、5%タイルを取得
maxs = [max(phreds[i]) for i in range(76)]
tops = [np.percentile(phreds[i], 95) for i in range(76)]
medians = [np.percentile(phreds[i], 50) for i in range(76)]
bottoms = [np.percentile(phreds[i], 5) for i in range(76)]
# プロット用に各パーセンタイルについて差分を取得
medians_fig = [x - y for x, y in zip(medians, bottoms)]
tops_fig = [x - y for x, y in zip(tops, medians)]
maxs_fig = [x - y for x, y in zip(maxs, tops)]

# 積み上げ折れ線グラフをプロット
fig, ax = plt.subplots(figsize=(16,9))
ax.stackplot(range(1, 77), (bottoms, medians_fig, tops_fig, maxs_fig))
ax.plot(range(1, 77), maxs, 'k-')

アラインされた領域を可視化した際に凹んでいた35bp-42bp辺りのスコアが顕著に低いことが分かります。

pysamstats でBAM/SAMファイルを扱う

pysamstatsは、BAM/SAMファイルについて、様々な統計値を算出するためのOSSになります。 インストールはconda 若しくはpipから可能です。こちらを参考にしてください。

下記の記事によくまとまっているので使用する際は参考にさせて頂きましょう。

kazumaxneo.hatenablog.com

カバレッジを取得

ポジションに関するカバレッジはload_coverageメソッドを用いることで取得することができます。
これを用いて、75000bp から 80000bpまでのカバレッジをプロットします。

bam = pysam.AlignmentFile('NA18489.chrom20.ILLUMINA.bwa.YRI.exome.20121211.bam', 'rb')
stats = pysamstats.load_coverage(bam, chrom='20', start=75000, end=80000)
plt.plot(stats.pos, stats.reads_all)
plt.show()

Python APIではオブジェクトの属性値として各種値が格納されていたのに対し、
CLIで取得する場合はタブ区切りで出力されます。

$ pysamstats --type coverage NA18489.chrom20.ILLUMINA.bwa.YRI.exome.20121211.bam | head
chrom   pos     reads_all       reads_pp
20      59995   1               1
20      59996   1               1
20      59997   2               2
20      59998   2               2
20      59999   2               2
20      60000   2               2
20      60001   2               2
20      60002   2               2
20      60003   2               2