[[ノート/ノート]]~
訪問者数 &counter(); 最終更新 &lastmod();
[[ノート/breseq/チュートリアルを試す１]] ＞ [[ノート/breseq/チュートリアルを試す２]] ＞ [[ノート/breseq/UT-Austinの進化データを試す]]~
[[ノート/breseq/Multi-Isolate１]] ＞ [[ノート/breseq/Multi-Isolate２]]
*breseq 2017-02-25 [#dafede4e]
Identification of mutations in laboratory evolved microbes from
next-generation sequencing data using breseq
3.1 Predicting mutations in a clonal sample (page 3)
$ breseq -j 4 -o Clonal_Output -r REL606.gbk SRR030257_1.fastq SRR030257_2.fastq
出力はHTMLで（見やすく）作られるので、[[この場所:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/SRR030257_output]] に置いた。
>（2017-04-08追加） 参照シーケンスをREL606とは別の「O157H7」にしてみた。~
チュートリアルでは親株のREL616との差異（変異）を見ようとして、bowtie2での参照シーケンスをREL606としてそれにリードを張り付けている。~
ここの追加では、試しに、参照シーケンスを変えてみる。NCBIが貼っているEscherichia coliの参照シーケンスのうち、適当に選んでみる。今回はO157-H7のアノテーション付きGeneBankデータを使ってみた。~
[[結果はここに置いた。:http://ginger.yy.is.sci.toho-u.ac.jp/~yamanouc/NGS2017/SRR030257-O157H7_output]]~
元の試行と極端に異なるのが、outputをまとめるための処理時間である。Summary Statisticsのタブを開くと、下の方にExecution Timeの表が載っている。O157-H7の方が全体に長くかかる傾向があるが、特にOutputのセクションが長い。
元の試行と極端に異なるのが、outputをまとめるための処理時間である。Summary Statisticsのタブを開くと、下の方にExecution Timeの表が載っている。O157-H7の方が全体に長くかかる傾向があるが、特にPolymorphic StatisticsとOutputのセクションが伸びている。（Polymorphic Statisticsが43倍、Outputが1430倍）
| step |REL606 | O157-H7 |
|Read and reference sequence file input|41 seconds|49 seconds|
|Read alignment to reference genome |23 minutes 46 seconds|36 minutes 11 seconds|
|Preprocessing alignments for candidate junction identification|2 minutes 43 seconds|6 minutes 50 seconds|
|Preliminary analysis of coverage distribution|1 minute 38 seconds|1 minute 41 seconds|
|Identifying junction candidates|14 seconds|17 seconds|
|Re-alignment to junction candidates|27 seconds|44 seconds|
|Resolving alignments with junction candidates|2 minutes 58 seconds|3 minutes 34 seconds|
|Creating BAM files|1 minute 26 seconds|1 minute 25 seconds|
|Tabulating error counts|30 seconds|49 seconds|
|Re-calibrating base error rates|1 second|2 seconds|
|Examining read alignment evidence| 5 minutes 51 seconds|15 minutes 45 seconds|
|Polymorphism statistics| 1 second| 43 seconds|
|Output | 21 seconds | 8 hours 20 minutes 28 seconds|
|Total | 40 minutes 37 seconds |9 hours 29 minutes 18 seconds|
breseqが何をやっているかを詳細に見たい。実行ログが吐き出されるので、以下に貼ってみた。これをていねいに追うと、何をどのように（＝どのようなパラメタを与えて）やっているかが見えてくる。
[yamanouc@ginger Clonal_Sample]$ breseq -j 4 -o Clonal_Output -r REL606.gbk SRR030257_1.fastq SRR030257_2.fastq
================================================================================
breseq 0.30.0 http://barricklab.org/breseq
Active Developers: Barrick JE, Deatherage DE
Contact: <jeffrey.e.barrick@gmail.com>
breseq is free software; you can redistribute it and/or modify it under the
terms the GNU General Public License as published by the Free Software
Foundation; either version 2, or (at your option) any later version.
Copyright (c) 2008-2010 Michigan State University
Copyright (c) 2011-2016 The University of Texas at Austin
If you use breseq in your research, please cite:
Deatherage, D.E., Barrick, J.E. (2014) Identification of mutations
in laboratory-evolved microbes from next-generation sequencing
data using breseq. Methods Mol. Biol. 1151: 165?188.
If you use structural variation (junction) predictions, please cite:
Barrick, J.E., Colburn, G., Deatherage D.E., Traverse, C.C.,
Strand, M.D., Borges, J.J., Knoester, D.B., Reba, A., Meyer, A.G.
(2014) Identifying structural variation in haploid microbial genomes
from short-read resequencing data using breseq. BMC Genomics 15:1039.
================================================================================
---> bowtie2 :: version 2.3.0 [/usr/local/bin/bowtie2]
---> R :: version 3.3.2 [/usr/bin/R]
+++ NOW PROCESSING Read and reference sequence file input
READ FILE::SRR030257_1
Converting/filtering FASTQ file...
Original base quality format: SANGER New format: SANGER
Original reads: 3800180 bases: 136806480
Filtered reads: 19 (?50% N) 60951 (?90% one base)
Analyzed reads: 3739210 bases: 134611560
READ FILE::SRR030257_2
Converting/filtering FASTQ file...
Original base quality format: SANGER New format: SANGER
Original reads: 3800180 bases: 136806480
Filtered reads: 369 (?50% N) 47630 (?90% one base)
Analyzed reads: 3752181 bases: 135078516
::TOTAL::
Original reads: 7600360 bases: 273612960
Analyzed reads: 7491391 bases: 269690076
[samtools] faidx Clonal_Output/data/reference.fasta
REFERENCE: REL606
LENGTH: 4629812
+++ NOW PROCESSING Read alignment to reference genome
[system] bowtie2-build -q Clonal_Output/data/reference.fasta Clonal_Output/02_reference_alignment/reference
[system] bowtie2 -t -p 4 --local -L 18 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-quals -k 2000 -i S,1,0.25 --score-min L,1,0.9 --reorder
-x Clonal_Output/02_reference_alignment/reference -U Clonal_Output/01_sequence_conversion/SRR030257_1.converted.fastq
-S Clonal_Output/02_reference_alignment/SRR030257_1.stage1.sam
--un Clonal_Output/02_reference_alignment/SRR030257_1.stage1.unmatched.fastq
Time loading reference: 00:00:00
Time loading forward index: 00:00:00
Time loading mirror index: 00:00:00
Multiseed full-index search: 00:00:33
3739210 reads; of these:
3739210 (100.00%) were unpaired; of these:
600703 (16.06%) aligned 0 times
3040853 (81.32%) aligned exactly 1 time
97654 (2.61%) aligned >1 times
83.94% overall alignment rate
Time searching: 00:00:33
Overall time: 00:00:33
[system] bowtie2 -t -p 4 --local -L 18 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-quals -k 2000 -i S,1,0.25 --score-min L,1,0.9 --reorder
-x Clonal_Output/02_reference_alignment/reference -U Clonal_Output/01_sequence_conversion/SRR030257_2.converted.fastq
-S Clonal_Output/02_reference_alignment/SRR030257_2.stage1.sam
--un Clonal_Output/02_reference_alignment/SRR030257_2.stage1.unmatched.fastq
Time loading reference: 00:00:00
Time loading forward index: 00:00:00
Time loading mirror index: 00:00:00
Multiseed full-index search: 00:00:34
3752181 reads; of these:
3752181 (100.00%) were unpaired; of these:
612595 (16.33%) aligned 0 times
3042265 (81.08%) aligned exactly 1 time
97321 (2.59%) aligned >1 times
83.67% overall alignment rate
Time searching: 00:00:34
Overall time: 00:00:34
[system] bowtie2 -t -p 4 --local -L 9 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-quals -k 2000 -i S,1,0.25 --score-min L,6,0.2 --reorder
-x Clonal_Output/02_reference_alignment/reference -U Clonal_Output/02_reference_alignment/SRR030257_1.stage1.unmatched.fastq
-S Clonal_Output/02_reference_alignment/SRR030257_1.stage2.matched.sam
Time loading reference: 00:00:00
Time loading forward index: 00:00:00
Time loading mirror index: 00:00:00
Multiseed full-index search: 00:10:58
600703 reads; of these:
600703 (100.00%) were unpaired; of these:
3722 (0.62%) aligned 0 times
10227 (1.70%) aligned exactly 1 time
586754 (97.68%) aligned >1 times
99.38% overall alignment rate
Time searching: 00:10:58
Overall time: 00:10:58
[system] bowtie2 -t -p 4 --local -L 9 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-quals -k 2000 -i S,1,0.25 --score-min L,6,0.2 --reorder
-x Clonal_Output/02_reference_alignment/reference -U Clonal_Output/02_reference_alignment/SRR030257_2.stage1.unmatched.fastq
-S Clonal_Output/02_reference_alignment/SRR030257_2.stage2.matched.sam
Time loading reference: 00:00:00
Time loading forward index: 00:00:00
Time loading mirror index: 00:00:00
Multiseed full-index search: 00:11:03
612595 reads; of these:
612595 (100.00%) were unpaired; of these:
5810 (0.95%) aligned 0 times
13420 (2.19%) aligned exactly 1 time
593365 (96.86%) aligned >1 times
99.05% overall alignment rate
Time searching: 00:11:03
Overall time: 00:11:03
[system] rm -f Clonal_Output/02_reference_alignment/reference*
[system] rm -f Clonal_Output/02_reference_alignment/SRR030257_1.stage1.unmatched.fastq
[system] rm -f Clonal_Output/02_reference_alignment/SRR030257_1.stage1.sam
[system] rm -f Clonal_Output/02_reference_alignment/SRR030257_2.stage1.unmatched.fastq
[system] rm -f Clonal_Output/02_reference_alignment/SRR030257_2.stage1.sam
[system] rm -f Clonal_Output/02_reference_alignment/SRR030257_1.stage2.matched.sam
[system] rm -f Clonal_Output/02_reference_alignment/SRR030257_2.stage2.matched.sam
+++ NOW PROCESSING Preprocessing alignments for candidate junction identification
Preprocessing alignments.
READ FILE::SRR030257_1
ALIGNED READ:100000
ALIGNED READ:200000
ALIGNED READ:300000
:
ALIGNED READ:3500000
ALIGNED READ:3600000
ALIGNED READ:3700000
READ FILE::SRR030257_2
ALIGNED READ:3800000
ALIGNED READ:3900000
ALIGNED READ:4000000
:
ALIGNED READ:7200000
ALIGNED READ:7300000
ALIGNED READ:7400000
Summary...
Aligned reads: 7481861
Read alignments: 19713263
Alignments split on indels: 21
Reads with alignments split on indels: 21
Split alignments: 6090236
Reads with split alignments: 601157
+++ NOW PROCESSING Preliminary analysis of coverage distribution
[samtools] import Clonal_Output/data/reference.fasta.fai Clonal_Output/03_candidate_junctions/best.sam Clonal_Output/03_candidate_junctions/best.unsorted.bam
[samtools] sort --threads 4 -o Clonal_Output/03_candidate_junctions/best.bam
-T Clonal_Output/03_candidate_junctions/best.bam Clonal_Output/03_candidate_junctions/best.unsorted.bam
[samtools] index Clonal_Output/03_candidate_junctions/best.bam Clonal_Output/03_candidate_junctions/best.bam.bai
REFERENCE: REL606
LENGTH: 4629812
POSITION:10000
POSITION:20000
POSITION:30000
:
POSITION:4610000
POSITION:4620000
[system] R --vanilla < /usr/local/bin/../share/breseq//coverage_distribution.r >
Clonal_Output/03_candidate_junctions/0.unique_only_coverage_distribution.tab.r.log
distribution_file=Clonal_Output/03_candidate_junctions/0.unique_only_coverage_distribution.tab
plot_file=Clonal_Output/03_candidate_junctions/0.coverage.pdf deletion_propagation_pr_cutoff=2.32374e-05
[system] rm -f Clonal_Output/03_candidate_junctions/best.unsorted.bam
[system] rm -f Clonal_Output/03_candidate_junctions/best.sam
[system] rm -f Clonal_Output/03_candidate_junctions/best.bam
[system] rm -f Clonal_Output/03_candidate_junctions/best.bam.bai
[system] rm -f Clonal_Output/03_candidate_junctions/0.coverage.pdf
[system] rm -f Clonal_Output/03_candidate_junctions/0.unique_only_coverage_distribution.tab
[system] rm -f Clonal_Output/03_candidate_junctions/0.unique_only_coverage_distribution.tab.r.log
+++ NOW PROCESSING Identifying junction candidates
READ FILE::SRR030257_1
ALIGNED READ:10000 CANDIDATE JUNCTIONS:2852
ALIGNED READ:20000 CANDIDATE JUNCTIONS:5715
ALIGNED READ:30000 CANDIDATE JUNCTIONS:8570
:
ALIGNED READ:270000 CANDIDATE JUNCTIONS:83719
ALIGNED READ:280000 CANDIDATE JUNCTIONS:87232
Passed alignment pairs examined: 100001
WARNING: Reached limit of 100000 passed alignment pairs.
Specify a greater value for --junction-alignment-pair-limit for more thorough junction prediction.
Taking top candidate junctions...
Minimum number to keep: 100
Maximum number to keep: 5000
Maximum length to keep: 462981 bases
Initial: Number = 76, Cumulative Length = 5723 bases
Testing Pos Hash Score = 11, Num = 1, Length = 75
Testing Pos Hash Score = 10, Num = 1, Length = 74
Testing Pos Hash Score = 9, Num = 2, Length = 147
Testing Pos Hash Score = 8, Num = 4, Length = 297
Testing Pos Hash Score = 7, Num = 2, Length = 149
Testing Pos Hash Score = 6, Num = 6, Length = 436
Testing Pos Hash Score = 5, Num = 3, Length = 220
Testing Pos Hash Score = 4, Num = 10, Length = 734
Testing Pos Hash Score = 3, Num = 10, Length = 742
Testing Pos Hash Score = 2, Num = 37, Length = 2849
Accepted: Number = 76, Pos Hash Score >= 2, Cumulative Length = 5723 bases
[samtools] faidx Clonal_Output/03_candidate_junctions/candidate_junction.fasta
[system] rm -f Clonal_Output/03_candidate_junctions/SRR030257_1.split.sam
[system] rm -f Clonal_Output/03_candidate_junctions/SRR030257_2.split.sam
+++ NOW PROCESSING Re-alignment to junction candidates
[system] bowtie2-build -q Clonal_Output/03_candidate_junctions/candidate_junction.fasta Clonal_Output/04_candidate_junction_alignment/candidate_junction
[system] bowtie2 -t -p 4 --local -L 10 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-quals -k 2000 -i S,1,0.25 --score-min L,1,0.70 --reorder
-x Clonal_Output/04_candidate_junction_alignment/candidate_junction -U Clonal_Output/01_sequence_conversion/SRR030257_1.converted.fastq
-S Clonal_Output/04_candidate_junction_alignment/SRR030257_1.candidate_junction.sam
Time loading reference: 00:00:00
Time loading forward index: 00:00:00
Time loading mirror index: 00:00:00
Multiseed full-index search: 00:00:13
3739210 reads; of these:
3739210 (100.00%) were unpaired; of these:
3735689 (99.91%) aligned 0 times
920 (0.02%) aligned exactly 1 time
2601 (0.07%) aligned >1 times
0.09% overall alignment rate
Time searching: 00:00:13
Overall time: 00:00:13
[system] bowtie2 -t -p 4 --local -L 10 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 --ignore-quals -k 2000 -i S,1,0.25 --score-min L,1,0.70 --reorder
-x Clonal_Output/04_candidate_junction_alignment/candidate_junction -U Clonal_Output/01_sequence_conversion/SRR030257_2.converted.fastq
-S Clonal_Output/04_candidate_junction_alignment/SRR030257_2.candidate_junction.sam
Time loading reference: 00:00:00
Time loading forward index: 00:00:00
Time loading mirror index: 00:00:00
Multiseed full-index search: 00:00:13
3752181 reads; of these:
3752181 (100.00%) were unpaired; of these:
3748811 (99.91%) aligned 0 times
840 (0.02%) aligned exactly 1 time
2530 (0.07%) aligned >1 times
0.09% overall alignment rate
Time searching: 00:00:13
Overall time: 00:00:13
[system] rm -f Clonal_Output/04_candidate_junction_alignment/candidate_junction*
+++ NOW PROCESSING Resolving alignments with junction candidates
READ FILE:SRR030257_1
READS:10000
READS:20000
:
READS:3720000
READS:3730000
READ FILE:SRR030257_2
READS:3740000
READS:3750000
:
READS:7480000
READS:7490000
[system] rm -f Clonal_Output/01_sequence_conversion/SRR030257_1.converted.fastq
[system] rm -f Clonal_Output/01_sequence_conversion/SRR030257_2.converted.fastq
[system] rm -f Clonal_Output/02_reference_alignment/SRR030257_1.reference.sam
[system] rm -f Clonal_Output/02_reference_alignment/SRR030257_2.reference.sam
[system] rm -f Clonal_Output/04_candidate_junction_alignment/SRR030257_1.candidate_junction.sam*
[system] rm -f Clonal_Output/04_candidate_junction_alignment/SRR030257_2.candidate_junction.sam*
+++ NOW PROCESSING Creating BAM files
[samtools] import Clonal_Output/03_candidate_junctions/candidate_junction.fasta.fai Clonal_Output/05_alignment_correction/junction.sam
Clonal_Output/06_bam/junction.unsorted.bam
[samtools] sort --threads 4 -o Clonal_Output/06_bam/junction.bam -T Clonal_Output/06_bam/junction.bam Clonal_Output/06_bam/junction.unsorted.bam
[system] rm -f Clonal_Output/06_bam/junction.unsorted.bam
[samtools] index Clonal_Output/06_bam/junction.bam Clonal_Output/06_bam/junction.bam.bai
[samtools] import Clonal_Output/data/reference.fasta.fai Clonal_Output/05_alignment_correction/reference.sam Clonal_Output/06_bam/reference.unsorted.bam
[samtools] sort --threads 4 -o Clonal_Output/data/reference.bam -T Clonal_Output/data/reference.bam Clonal_Output/06_bam/reference.unsorted.bam
[system] rm -f Clonal_Output/06_bam/reference.unsorted.bam
[samtools] index Clonal_Output/data/reference.bam Clonal_Output/data/reference.bam.bai
[system] rm -f Clonal_Output/05_alignment_correction/reference.sam
[system] rm -f Clonal_Output/05_alignment_correction/junction.sam
+++ NOW PROCESSING Tabulating error counts
REFERENCE: REL606
LENGTH: 4629812
POSITION:10000
POSITION:20000
POSITION:30000
:
POSITION:4610000
POSITION:4620000
+++ NOW PROCESSING Re-calibrating base error rates
[system] R --vanilla < /usr/local/bin/../share/breseq//coverage_distribution.r >
Clonal_Output/07_error_calibration/0.unique_only_coverage_distribution.tab.r.log
distribution_file=Clonal_Output/07_error_calibration/0.unique_only_coverage_distribution.tab
plot_file=Clonal_Output/output/calibration/0.unique_coverage.pdf
deletion_propagation_pr_cutoff=2.32374e-05
[system] R --vanilla in_file=Clonal_Output/07_error_calibration/base_qual_error_prob.SRR030257_1.tab
out_file=Clonal_Output/output/calibration/SRR030257_1.error_rates.pdf <
/usr/local/bin/../share/breseq//plot_error_rate.r >
Clonal_Output/07_error_calibration/SRR030257_1.plot_error_rate.log
[system] R --vanilla in_file=Clonal_Output/07_error_calibration/base_qual_error_prob.SRR030257_2.tab
out_file=Clonal_Output/output/calibration/SRR030257_2.error_rates.pdf <
/usr/local/bin/../share/breseq//plot_error_rate.r >
Clonal_Output/07_error_calibration/SRR030257_2.plot_error_rate.log
+++ NOW PROCESSING Examining read alignment evidence
REFERENCE: REL606
LENGTH: 4629812
POSITION:10000
POSITION:20000
:
POSITION:4610000
POSITION:4620000
+++ NOW PROCESSING Polymorphism statistics
[system] R --vanilla total_length=4629812 in_file=Clonal_Output/08_mutation_identification/polymorphism_statistics_input.tab
out_file=Clonal_Output/08_mutation_identification/polymorphism_statistics_output.tab qual_file=Clonal_Output/08_mutation_identification/error_counts.tab <
/usr/local/bin/../share/breseq//polymorphism_statistics.r >
Clonal_Output/08_mutation_identification/polymorphism_statistics_output.log
+++ NOW PROCESSING Output
Creating merged genome diff evidence file...
Predicting mutations from evidence...
Annotating mutations...
Drawing coverage plots...
Creating coverage plot for region: REL606:1-4629812
Creating coverage plot for region: REL606:546958-555936
Creating coverage plot for region: REL606:2031674-2054972
Creating coverage plot for region: REL606:3289962-3289977
Creating coverage plot for region: REL606:3893578-3901930
Creating read alignment for region: REL606:546957-546957
Creating read alignment for region: REL606:555937-555937
Creating read alignment for region: REL606:2031673-2031673
Creating read alignment for region: REL606:2054973-2054973
Creating read alignment for region: REL606:3289961-3289961
Creating read alignment for region: REL606:3289978-3289978
Creating read alignment for region: REL606:3893577-3893577
Creating read alignment for region: REL606:3901931-3901931
Creating read alignment for region: REL606:161041.0-161041.0
Creating read alignment for region: REL606:380188.0-380188.0
Creating read alignment for region: REL606:430835.0-430835.0
Creating read alignment for region: REL606:475288.1-475288.1
Creating read alignment for region: REL606:649391.0-649391.0
Creating read alignment for region: REL606:668787.0-668787.0
Creating read alignment for region: REL606:683496.0-683496.0
Creating read alignment for region: REL606:910316.0-910316.0
Creating read alignment for region: REL606:1004251.0-1004251.0
Creating read alignment for region: REL606:1248380.0-1248380.0
Creating read alignment for region: REL606:1286699.0-1286699.0
Creating read alignment for region: REL606:1329516.0-1329516.0
Creating read alignment for region: REL606:1976879.0-1976879.0
Creating read alignment for region: REL606:2082685.0-2082685.0
Creating read alignment for region: REL606:2732014.0-2732014.0
Creating read alignment for region: REL606:3045069.0-3045069.0
Creating read alignment for region: REL606:3155168.0-3155168.0
Creating read alignment for region: REL606:3200643.0-3200643.0
Creating read alignment for region: REL606:3248957.0-3248957.0
Creating read alignment for region: REL606:3288092.0-3288092.0
Creating read alignment for region: REL606:3339158.0-3339158.0
Creating read alignment for region: REL606:3369432.0-3369432.0
Creating read alignment for region: REL606:3370027.0-3370027.0
Creating read alignment for region: REL606:3483047.0-3483047.0
Creating read alignment for region: REL606:3490112.0-3490112.0
Creating read alignment for region: REL606:3762741.0-3762741.0
Creating read alignment for region: REL606:3875625.1-3875625.1
Creating read alignment for region: REL606:3893550.1-3893550.1
Creating read alignment for region: REL606:4100655.0-4100655.0
Creating read alignment for region: REL606:4126699.0-4126699.0
Creating read alignment for region: REL606:4201911.0-4201911.0
Creating read alignment for region: REL606:4616396.0-4616396.0
Creating read alignment for region: REL606:87812.0-87812.0
Creating read alignment for region: REL606:93946.0-93946.0
Creating read alignment for region: REL606:139812.0-139812.0
Creating read alignment for region: REL606:402056.0-402056.0
Creating read alignment for region: REL606:986054.0-986054.0
Creating read alignment for region: REL606:1318108.0-1318108.0
Creating read alignment for region: REL606:1509936.0-1509936.0
Creating read alignment for region: REL606:1812720.0-1812720.0
Creating read alignment for region: REL606:1854388.0-1854388.0
Creating read alignment for region: REL606:2167569.0-2167569.0
Creating read alignment for region: REL606:2698093.0-2698093.0
Creating read alignment for region: REL606:3357000.0-3357000.0
Creating read alignment for region: REL606:3357002.0-3357002.0
Creating read alignment for region: REL606:3357010.0-3357010.0
Creating read alignment for region: REL606:3357012.0-3357012.0
Creating read alignment for region: REL606:3357016.0-3357016.0
Creating read alignment for region: REL606:4080986.0-4080986.0
Creating read alignment for region: REL606:4462991.0-4462991.0
Creating read alignment for region: REL606:4462995.0-4462995.0
Creating read alignment for region: REL606:4462997.0-4462997.0
Creating read alignment for region: REL606__1__1__REL606__4629812__-1__0____36__36__0__0:36-37
Creating read alignment for region: REL606:1-1
Creating read alignment for region: REL606:4629812-4629812
Creating read alignment for region: REL606__16972__1__REL606__588495__1__0____36__36__0__1:36-37
Creating read alignment for region: REL606:16972-16972
Creating read alignment for region: REL606:588495-588495
Creating read alignment for region: REL606__16976__-1__REL606__590471__-1__2____36__36__0__1:36-39
Creating read alignment for region: REL606:16974-16974
Creating read alignment for region: REL606:590471-590471
Creating read alignment for region: REL606__72407__1__REL606__3650758__1__7____36__36__0__0:36-44
Creating read alignment for region: REL606:72407-72407
Creating read alignment for region: REL606:3650765-3650765
Creating read alignment for region: REL606__391964__-1__REL606__3875885__-1__-6__TTACTG__36__36__0__0:36-43
Creating read alignment for region: REL606:391964-391964
Creating read alignment for region: REL606:3875885-3875885
Creating read alignment for region: REL606__547699__-1__REL606__555924__1__0____36__36__1__1:36-37
Creating read alignment for region: REL606:547699-547699
Creating read alignment for region: REL606:555924-555924
Creating read alignment for region: REL606__589884__1__REL606__590471__-1__2____36__36__1__1:36-39
Creating read alignment for region: REL606:589886-589886
Creating read alignment for region: REL606:590471-590471
Creating read alignment for region: REL606__664688__1__REL606__1270658__1__2____36__36__1__0:36-39
Creating read alignment for region: REL606:664688-664688
Creating read alignment for region: REL606:1270660-1270660
Creating read alignment for region: REL606__809428__1__REL606__4141021__-1__5____36__36__0__0:36-42
Creating read alignment for region: REL606:809428-809428
Creating read alignment for region: REL606:4141016-4141016
Creating read alignment for region: REL606__1154778__1__REL606__4215353__-1__11____36__36__0__0:36-48
Creating read alignment for region: REL606:1154778-1154778
Creating read alignment for region: REL606:4215342-4215342
Creating read alignment for region: REL606__1733645__1__REL606__2775877__-1__2____36__36__0__1:36-39
Creating read alignment for region: REL606:1733647-1733647
Creating read alignment for region: REL606:2775877-2775877
Creating read alignment for region: REL606__1733649__-1__REL606__2774435__1__0____36__36__0__1:36-37
Creating read alignment for region: REL606:1733649-1733649
Creating read alignment for region: REL606:2774435-2774435
Creating read alignment for region: REL606__1802331__1__REL606__3944627__1__7____36__36__0__0:36-44
Creating read alignment for region: REL606:1802331-1802331
Creating read alignment for region: REL606:3944634-3944634
Creating read alignment for region: REL606__1821523__1__REL606__2775877__-1__2____36__36__0__1:36-39
Creating read alignment for region: REL606:1821525-1821525
Creating read alignment for region: REL606:2775877-2775877
Creating read alignment for region: REL606__1821528__-1__REL606__2774435__1__1____36__36__0__1:36-38
Creating read alignment for region: REL606:1821527-1821527
Creating read alignment for region: REL606:2774435-2774435
Creating read alignment for region: REL606__2103888__1__REL606__2103918__-1__11____36__36__0__0:36-48
Creating read alignment for region: REL606:2103888-2103888
Creating read alignment for region: REL606:2103907-2103907
Creating read alignment for region: REL606__2429311__1__REL606__2792997__-1__8____36__36__0__0:36-45
Creating read alignment for region: REL606:2429311-2429311
Creating read alignment for region: REL606:2792989-2792989
Creating read alignment for region: REL606__2772854__1__REL606__4524530__-1__3____36__36__1__0:36-40
Creating read alignment for region: REL606:2772854-2772854
Creating read alignment for region: REL606:4524527-4524527
Creating read alignment for region: REL606__2774196__-1__REL606__4524519__1__3____36__36__1__0:36-40
Creating read alignment for region: REL606:2774196-2774196
Creating read alignment for region: REL606:4524522-4524522
Creating read alignment for region: REL606__2775877__-1__REL606__3015769__1__2____36__36__1__0:36-39
Creating read alignment for region: REL606:2775877-2775877
Creating read alignment for region: REL606:3015771-3015771
Creating read alignment for region: REL606__3015774__-1__REL606__3894996__-1__0____36__36__0__1:36-37
Creating read alignment for region: REL606:3015774-3015774
Creating read alignment for region: REL606:3894996-3894996
Creating read alignment for region: REL606__3289961__-1__REL606__3289978__1__0____36__36__0__0:36-37
Creating read alignment for region: REL606:3289961-3289961
Creating read alignment for region: REL606:3289978-3289978
Creating read alignment for region: REL606__3630716__-1__REL606__3652533__-1__0____36__36__1__1:36-37
Creating read alignment for region: REL606:3630716-3630716
Creating read alignment for region: REL606:3652533-3652533
Creating read alignment for region: REL606__3651091__1__REL606__3741531__-1__0____36__36__1__1:36-37
Creating read alignment for region: REL606:3651091-3651091
Creating read alignment for region: REL606:3741531-3741531
Creating read alignment for region: REL606__3652533__-1__REL606__3901929__1__2____36__36__1__0:36-39
Creating read alignment for region: REL606:3652533-3652533
Creating read alignment for region: REL606:3901931-3901931
Creating read alignment for region: REL606__3749572__-1__REL606__4157637__-1__6____36__36__0__0:36-43
Creating read alignment for region: REL606:3749572-3749572
Creating read alignment for region: REL606:4157631-4157631
Creating read alignment for region: REL606__3813062__-1__REL606__4157637__-1__6____36__36__0__0:36-43
Creating read alignment for region: REL606:3813062-3813062
Creating read alignment for region: REL606:4157631-4157631
Creating read alignment for region: REL606__3944627__1__REL606__4157637__-1__7____36__36__0__0:36-44
Creating read alignment for region: REL606:3944627-3944627
Creating read alignment for region: REL606:4157630-4157630
Creating index HTML table...
[system] rm -f Clonal_Output/01_sequence_conversion/*.trims
[system] rm -f Clonal_Output/03_candidate_junctions/candidate_junction.fasta
[system] rm -f Clonal_Output/03_candidate_junctions/candidate_junction.fasta.fai
[system] rm -f Clonal_Output/05_alignment_correction/jc_evidence.gd
[system] rm -f Clonal_Output/06_bam/junction.bam
[system] rm -f Clonal_Output/06_bam/junction.bam.bai
+++ SUCCESSFULLY COMPLETED
**出力結果 [#n5e322ee]
ディレクトリ Colonal Outputを見る。出力はWebで見やすいようになっている。
ここでは、サーバーからディレクトリ Colonal Output全体をローカルWinにコピーして
その中で見て回る。
index.htmlからはじまる。
&ref(BRESEQ __ Mutation Predictions.pdf);
&ref(BRESEQ __ Marginal Predictions.pdf);
&ref(BRESEQ __ Summary Statistics.pdf);
** ログファイルから、処理を追ってみる [#xb1517e5]
&ref(2017-03-14_breseq_cmdlineの解説.txt);
// 01 sequence conversion
ファイル SRR030257_1
ファイル SRR030257_2
|
+------------------+
| シーケンス変換 | 入力をfastaやGFF3に変換する、quality scoreをSangerにする
+------------------+
|
+-------------------+
| SAMのfaidxを作る| samtools_faidx(reference.fasta)
+-------------------+
|
+--------------------+
| trimファイルを作る| calculate_trims(..)
+--------------------+
|
// 02 reference alignment
|
BUILDING MAPPING INDEX
+------------------------------------------------------------------------------+
| bowtie2-build -q reference.fasta 02_reference_alignment/reference |
+------------------------------------------------------------------------------+
|
STAGE-1 ALIGNMENT
+------------------------------------------------------------------------------+
| bowtie2 -t -p 4 --local -L 18 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 |
| --ignore-quals -k 2000 -i S,1,0.25 --score-min L,1,0.9 --reorder |
| -x 02_reference_alignment/reference |
| -U 01_sequence_conversion/SRR030257_1.converted.fastq |
| -S 02_reference_alignment/SRR030257_1.stage1.sam |
| --un 02_reference_alignment/SRR030257_1.stage1.unmatched.fastq |
+------------------------------------------------------------------------------+
| 83.94% overall alignment rate
+------------------------------------------------------------------------------+
| bowtie2 -t -p 4 --local -L 18 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 |
| --ignore-quals -k 2000 -i S,1,0.25 --score-min L,1,0.9 --reorder |
| -x 02_reference_alignment/reference |
| -U 01_sequence_conversion/SRR030257_2.converted.fastq |
| -S 02_reference_alignment/SRR030257_2.stage1.sam |
| --un 02_reference_alignment/SRR030257_2.stage1.unmatched.fastq |
+-------------------------------------------------------------------------------+
| 83.67% overall alignment rate
|
STAGE-2 ALGINMENT
+-------------------------------------------------------------------------------+
| bowtie2 -t -p 4 --local -L 9 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 |
| --ignore-quals -k 2000 -i S,1,0.25 --score-min L,6,0.2 --reorder |
| -x 02_reference_alignment/reference |
| -U 02_reference_alignment/SRR030257_1.stage1.unmatched.fastq |
| -S 02_reference_alignment/SRR030257_1.stage2.matched.sam +
+--------------------------------------------------------------------------------+
| 99.38% overall alignment rate
+--------------------------------------------------------------------------------+
| bowtie2 -t -p 4 --local -L 9 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 |
| --ignore-quals -k 2000 -i S,1,0.25 --score-min L,6,0.2 --reorder |
| -x 02_reference_alignment/reference |
| -U 02_reference_alignment/SRR030257_2.stage1.unmatched.fastq |
| -S 02_reference_alignment/SRR030257_2.stage2.matched.sam |
+--------------------------------------------------------------------------------+
| 99.05% overall alignment rate
|
MERGE SAM FILES
+---------------------------------------------------------------------------------+
| PreprocessAlignments::merge_sort_sam_files( |
| stage1_reference_sam_file_name, stage2_reference_sam_file_name, |
| reference_sam_file_name ); |
+---------------------------------------------------------------------------------+
|
// 03_candidate_junctions
|
Preprocessing alignments for candidate junction identification
+----------------------------------------------------------------------------------+
| PreprocessAlignments::preprocess_alignments( |
| settings, summary, ref_seq_info); |
+----------------------------------------------------------------------------------+
| Aligned reads: 7481861
| Read alignments: 19713263
| Alignments split on indels: 21
| Reads with alignments split on indels: 21
| Split alignments: 6090236
| Reads with split alignments: 601157
|
Preliminary analysis of coverage distribution
+-----------------------------------------------------------------------------------+
| samtools_import(data/reference.fasta.fai, |
| 03_candidate_junctions/best.sam, |
| 03_candidate_junctions/best.unsorted.bam |
+-----------------------------------------------------------------------------------+
|
+-----------------------------------------------------------------------------------+
| samtools_sort( --threads 4 -o 03_candidate_junctions/best.bam, |
| -T 03_candidate_junctions/best.bam, |
| 03_candidate_junctions/best.unsorted.bam |
+-----------------------------------------------------------------------------------+
|
+----------------------------------------------------------------------------------------+
| samtools_index( |
| 03_candidate_junctions/best.bam 03_candidate_junctions/best.bam.bai ) |
+----------------------------------------------------------------------------------------+
|
+-----------------------------------------------------------------------------------+
| CoverageDistribution::analyze_unique_coverage_distributions(settings, |
| summary, ref_seq_info, settings.coverage_junction_plot_file_name, |
| settings.coverage_junction_distribution_file_name, |
| settings.coverage_junction_done_file_name ); |
+--+ この中でRを呼んで処理 |
| R --vanilla |
| < /usr/local/bin/../share/breseq//coverage_distribution.r |
| > 03_candidate_junctions/0.unique_only_coverage_distribution.tab.r.log |
| distribution_file= |
| 03_candidate_junctions/0.unique_only_coverage_distribution.tab |
| plot_file=03_candidate_junctions/0.coverage.pdf |
| deletion_propagation_pr_cutoff=2.32374e-05 |
+-------------------------------------------------------------------------------+
|
Identifying junction candidates
+-----------------------------------------------------------------------------------+
| CandidateJunctions::identify_candidate_junctions( |
| settings, summary, ref_seq_info); |
+-----------------------------------------------------------------------------------+
| Passed alignment pairs examined: 100001
| WARNING: Reached limit of 100000 passed alignment pairs.
| Specify a greater value for --junction-alignment-pair-limit
| for more thorough junction prediction.
| Taking top candidate junctions...
| Minimum number to keep: 100
| Maximum number to keep: 5000
| Maximum length to keep: 462981 bases
| Initial: Number = 76, Cumulative Length = 5723 bases
| Accepted: Number = 76, Pos Hash Score >= 2,
| Cumulative Length = 5723 bases
|
+-----------------------------------------------------------------------------------+
| samtools_faidx 03_candidate_junctions/candidate_junction.fasta |
+-----------------------------------------------------------------------------------+
|
Re-alignment to junction candidates
+-----------------------------------------------------------------------------------+
| bowtie2-build -q 03_candidate_junctions/candidate_junction.fasta |
| 04_candidate_junction_alignment/candidate_junction |
+-----------------------------------------------------------------------------------+
|
+-----------------------------------------------------------------------------------+
| bowtie2 -t -p 4 --local -L 10 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 |
| --ignore-quals -k 2000 -i S,1,0.25 --score-min L,1,0.70 --reorder |
| -x 04_candidate_junction_alignment/candidate_junction |
| -U 01_sequence_conversion/SRR030257_1.converted.fastq |
| -S 04_candidate_junction_alignment/SRR030257_1.candidate_junction.sam |
+-----------------------------------------------------------------------------------+
| 0.09% overall alignment rate
+-----------------------------------------------------------------------------------+
| bowtie2 -t -p 4 --local -L 10 --ma 1 --mp 3 --np 0 --rdg 2,3 --rfg 2,3 |
| --ignore-quals -k 2000 -i S,1,0.25 --score-min L,1,0.70 --reorder |
| -x 04_candidate_junction_alignment/candidate_junction |
| -U 01_sequence_conversion/SRR030257_2.converted.fastq |
| -S 04_candidate_junction_alignment/SRR030257_2.candidate_junction.sam
+------------------------------------------------------------------------------------+
| 0.09% overall alignment rate
// 05 alignment_correction <-- line 1609
// * Resolve matches to new junction candidates
|
Resolving alignments with junction candidates
+------------------------------------------------------------------------------------+
| resolve_alignments( settings, summary, ref_seq_info, junction_prediction, |
| settings.read_files ); |
+------------------------------------------------------------------------------------+
|
// 05 bam
// * Create BAM read alignment database files
|
Creating BAM files
+-------------------------------------------------------------------------------------+
| samtools_import( 03_candidate_junctions/candidate_junction.fasta.fai |
| 05_alignment_correction/junction.sam 06_bam/junction.unsorted.bam ) |
+--------------------------------------------------------------------------------------+
|
+--------------------------------------------------------------------------------------+
| samtools_sort( --threads 4 -o 06_bam/junction.bam |
| -T 06_bam/junction.bam 06_bam/junction.unsorted.bam ) |
+--------------------------------------------------------------------------------------+
|
+--------------------------------------------------------------------------------------+
| samtools_index( 06_bam/junction.bam 06_bam/junction.bam.bai ) |
+--------------------------------------------------------------------------------------+
|
+--------------------------------------------------------------------------------------+
| samtools_import ( data/reference.fasta.fai |
| 05_alignment_correction/reference.sam 06_bam/reference.unsorted.bam )|
+--------------------------------------------------------------------------------------+
|
+--------------------------------------------------------------------------------------+
| samtools_sort( --threads 4 -o data/reference.bam |
| -T data/reference.bam /06_bam/reference.unsorted.bam ) |
+--------------------------------------------------------------------------------------+
|
+--------------------------------------------------------------------------------------+
| samtools_index( data/reference.bam data/reference.bam.bai ) |
+--------------------------------------------------------------------------------------+
|
Tabulating error counts
+-------------------------------------------------------------------------------------+
| error_count( settings, summary, reference_bam_file_name, |
| reference_fasta_file_name, settings.error_calibration_path, |
| settings.read_files.base_names(), coverage=true, errors=true, |
| preprocess=false, settings.base_quality_cutoff, // minimum quality score|
| read_set=num_read_files,obs_base,ref_base,quality=num_qual ); |
+-------------------------------------------------------------------------------------+
|
// Calculate error rates
|
Re-calibrating base error rates
+------------------------------------------------------------------------------------+
| CoverageDistribution::analyze_unique_coverage_distributions ( |
| settings, summary, ref_seq_info, |
| settings.unique_only_coverage_plot_file_name, |
| settings.unique_only_coverage_distribution_file_name, |
| "NULL" // means to never delete intermediates ); |
+--+ この中でfitを呼ぶ ------------------------------------------------------------+
| vector<string> lines = dist.fit(settings, |
| unique_only_coverage_distribution_file_name, |
| unique_only_coverage_plot_file_name, |
| deletion_propagation_pr_cutoff ); |
+--+ この中でRを呼ぶ ---------------------------------------------------------+
| R --vanilla |
| < /usr/local/bin/../share/breseq//coverage_distribution.r |
| > 07_error_calibration/0.unique_only_coverage_distribution.tab.r.log |
| distribution_file=07_error_calibration/0.unique_only_coverage_distribution.tab |
| plot_file=output/calibration/0.unique_coverage.pdf |
| deletion_propagation_pr_cutoff=2.32374e-05 |
+--------------------------------------------------------------------------- +
|
+----------------------------------------------------------------------------+
| R --vanilla |
| in_file=07_error_calibration/base_qual_error_prob.SRR030257_1.tab|
| out_file=output/calibration/SRR030257_1.error_rates.pdf |
| < /usr/local/bin/../share/breseq//plot_error_rate.r |
| > 07_error_calibration/SRR030257_1.plot_error_rate.log |
+----------------------------------------------------------------------------+
|
+----------------------------------------------------------------------------+
| R --vanilla |
| in_file=07_error_calibration/base_qual_error_prob.SRR030257_2.tab|
| out_file=output/calibration/SRR030257_2.error_rates.pdf |
| < /usr/local/bin/../share/breseq//plot_error_rate.r |
| > 07_error_calibration/SRR030257_2.plot_error_rate.log |
+---------------------------------------------------------------------------+
|
Examining read alignment evidence
+-------------------------------------------------------------------------------------+
| identify_mutations( settings, summary, reference_bam_file_name, |
| reference_fasta_file_name, ra_mc_genome_diff_file_name, |
| deletion_propagation_cutoffs, deletion_seed_cutoffs, |
| settings.mutation_log10_e_value_cutoff, // mutation_cutoff |
| settings.polymorphism_log10_e_value_cutoff, // polymorphism_cutoff |
| settings.polymorphism_frequency_cutoff, //polymorphism_frequency_cutoff|
| settings.polymorphism_precision_decimal, |
| settings.polymorphism_precision_places, |
| settings.print_mutation_identification_per_position_file //per_position_file );|
+--------------------------------------------------------------------------------------+
|
Polymorphism statistics
+--------------------------------------------------------------------------------------+
| ref_seq_info.polymorphism_statistics(settings, summary) |
+--+ この中でpolymorphism_statisticsを呼ぶ ------------------------------------+
| cReferenceSequences::polymorphism_statistics |
+--+ この中でRを呼ぶ --------------------------------------------------------+
| R --vanilla total_length=4629812 |
| in_file=08_mutation_identification/polymorphism_statistics_input.tab|
| out_file=08_mutation_identification/polymorphism_statistics_output.tab|
| qual_file=08_mutation_identification/error_counts.tab |
| < /usr/local/bin/../share/breseq//polymorphism_statistics.r |
| > /08_mutation_identification/polymorphism_statistics_output.log |
+------------------------------------------------------------------------------+
|
Output
+---------------------------------------------------------------------------+
| Creating merged genome diff evidence file... |
| Predicting mutations from evidence... |
| Annotating mutations... |
| Drawing coverage plots... |
| Creating index HTML table... |
+----------------------------------------------------------------------------+
|
SUCCESSFULLY COMPLETED
![[Toho University]](image/CLogo_Kihon_process.jpg "[Toho University]")
![[YYLAB]](image/yamanouc2.png "[YYLAB]")
