インストール手順は大抵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から可能です。こちらを参考にしてください。
下記の記事によくまとまっているので使用する際は参考にさせて頂きましょう。
カバレッジを取得
ポジションに関するカバレッジは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
