~ubuntu-branches/ubuntu/utopic/idba/utopic-proposed

« back to all changes in this revision

Viewing changes to src/tools/idba_tran_test.cpp

Committer: Package Import Robot
Author(s): Andreas Tille
Date: 2014-02-13 13:57:55 UTC
Revision ID: package-import@ubuntu.com-20140213135755-c97jpl9wl5yzfwfw

Tags: upstream-1.1.1

Import upstream version 1.1.1

files added:

AUTHORS

ChangeLog

Makefile.am

Makefile.in

NEWS

README

aclocal.m4

bin/Makefile.am

bin/Makefile.in

config.h.in

configure

configure.ac

depcomp

install-sh

lib/Makefile.am

lib/Makefile.in

missing

script

script/Makefile.am

script/Makefile.in

script/run-unittest.py

script/scan.py

script/validate_blat

script/validate_blat_parallel

src/assembly

src/assembly/assembly_utility.cpp

src/assembly/assembly_utility.h

src/assembly/local_assembler.cpp

src/assembly/local_assembler.h

src/basic

src/basic/atomic_integer.h

src/basic/bit_operation.h

src/basic/functional.h

src/basic/hash.h

src/basic/histgram.h

src/basic/kmer.h

src/basic/math.cpp

src/basic/math.h

src/basic/pool.h

src/container

src/container/hash_map.h

src/container/hash_set.h

src/container/hash_table.h

src/container/managed_list.h

src/graph

src/graph/bit_edges.h

src/graph/contig_builder.h

src/graph/contig_graph.cpp

src/graph/contig_graph.h

src/graph/contig_graph_branch_group.cpp

src/graph/contig_graph_branch_group.h

src/graph/contig_graph_path.h

src/graph/contig_graph_vertex.h

src/graph/contig_info.cpp

src/graph/contig_info.h

src/graph/hash_graph.cpp

src/graph/hash_graph.h

src/graph/hash_graph_branch_group.cpp

src/graph/hash_graph_branch_group.h

src/graph/hash_graph_path.h

src/graph/hash_graph_vertex.h

src/graph/scaffold_graph.cpp

src/graph/scaffold_graph.h

src/graph/scaffold_graph_path.h

src/graph/scaffold_graph_vertex.h

src/graph/vertex_status.h

src/misc

src/misc/blat_record.cpp

src/misc/blat_record.h

src/misc/hash_aligner.cpp

src/misc/hash_aligner.h

src/misc/log.cpp

src/misc/log.h

src/misc/options_description.cpp

src/misc/options_description.h

src/misc/utils.cpp

src/misc/utils.h

src/release

src/release/filterfa.cpp

src/release/fq2fa.cpp

src/release/idba.cpp

src/release/idba_hybrid.cpp

src/release/idba_tran.cpp

src/release/idba_ud.cpp

src/release/parallel_blat.cpp

src/release/raw_n50.cpp

src/release/sim_reads.cpp

src/release/sim_reads_tran.cpp

src/release/split_scaffold.cpp

src/release/validate_contigs_blat.cpp

src/sequence

src/sequence/compact_sequence.cpp

src/sequence/compact_sequence.h

src/sequence/read_library.h

src/sequence/sequence.cpp

src/sequence/sequence.h

src/sequence/sequence_io.cpp

src/sequence/sequence_io.h

src/sequence/sequence_reader.cpp

src/sequence/sequence_reader.h

src/sequence/sequence_writer.cpp

src/sequence/sequence_writer.h

src/sequence/short_sequence.cpp

src/sequence/short_sequence.h

src/tools

src/tools/fa2fq.cpp

src/tools/filter_blat.cpp

src/tools/filter_contigs.cpp

src/tools/idba_tran_test.cpp

src/tools/parallel_rna_blat.cpp

src/tools/print_graph.cpp

src/tools/sample_reads.cpp

src/tools/scaffold.cpp

src/tools/shuffle_reads.cpp

src/tools/sort_psl.cpp

src/tools/sort_reads.cpp

src/tools/split_fa.cpp

src/tools/split_fq.cpp

src/tools/test.cpp

src/tools/validate_component.cpp

src/tools/validate_contigs_mummer.cpp

src/tools/validate_reads_blat.cpp

src/tools/validate_rna.cpp

test

test/Makefile.am

test/Makefile.in

Show diffs side-by-side

added added

removed removed

src/tools/idba_tran_test.cpp

/**

* @file idba_tran.cpp

* @brief An iterative de Bruijn graph assembler for transcriptome sequencing reads.

* @author Yu Peng (ypeng@cs.hku.hk)

* @version 1.0.13

* @date 2012-09-14

#include <cmath>

#include <cstdio>

#include <deque>

#include <fstream>

#include <iostream>

#include <sstream>

#include <vector>

#include <sys/stat.h>

#include <unistd.h>

#include "assembly/assembly_utility.h"

#include "assembly/local_assembler.h"

#include "basic/bit_operation.h"

#include "basic/histgram.h"

#include "graph/contig_graph.h"

#include "graph/hash_graph.h"

#include "misc/hash_aligner.h"

#include "misc/log.h"

#include "misc/options_description.h"

#include "misc/utils.h"

#include "sequence/sequence.h"

#include "sequence/sequence_io.h"

#include "sequence/short_sequence.h"

using namespace std;

struct IDBAOption

{

string directory;

string read_file;

string long_read_file;

int mink;

int maxk;

int step;

int inner_mink;

int inner_step;

int prefix_length;

int min_count;

int min_support;

int min_contig;

double similar;

int max_mismatch;

int seed_kmer_size;

int num_threads;

int min_pairs;

int max_gap;

bool is_no_local;

bool is_use_all;

bool is_no_coverage;

bool is_no_correct;

bool is_pre_correction;

bool is_no_internal;

string reference;

IDBAOption()

{

directory = "out";

mink = 20;

maxk = 60;

step = 10;

inner_mink = 10;

inner_step = 5;

prefix_length = 3;

min_count = 2;

min_support = 1;

min_contig = 200;

similar = 0.95;

max_mismatch = 3;

seed_kmer_size = 30;

num_threads = 0;

min_pairs = 3;

max_gap = 0;

is_no_local = false;

is_use_all = false;

is_no_coverage = false;

is_no_correct = false;

is_pre_correction = false;

}

string log_file()

{ return directory + "/log"; }

string kmer_file()

{ return directory + "/kmer"; }

string align_file(int kmer_size)

{ return directory + FormatString("/align-%d", kmer_size); }

string align_local_file(int kmer_size)

100

{ return directory + FormatString("/align-local-%d", kmer_size); }

101

102

string graph_file(int kmer_size)

103

{ return directory + FormatString("/graph-%d.fa", kmer_size); }

104

105

string contig_file(int kmer_size)

106

{ return directory + FormatString("/contig-%d.fa", kmer_size); }

107

108

string contig_info_file(int kmer_size)

109

{ return directory + FormatString("/contig-info-%d.fa", kmer_size); }

110

111

string local_contig_file(int kmer_size)

112

{ return directory + FormatString("/local-contig-%d.fa", kmer_size); }

113

114

string contig_file()

115

{ return directory + "/contig.fa"; }

116

117

string scaffold_file()

118

{ return directory + "/scaffold.fa"; }

119

120

string ref_contig_file()

121

{ return directory + "/ref_contig.fa"; }

122

123

string transcript_file(int kmer_size)

124

{ return directory + FormatString("/transcript-%d.fa", kmer_size); }

125

126

string label_transcript_file(int kmer_size)

127

{ return directory + FormatString("/label-transcript-%d.fa", kmer_size); }

128

129

string component_file(int kmer_size)

130

{ return directory + FormatString("/component-%d", kmer_size); }

131

};

132

133

AssemblyInfo assembly_info;

134

IDBAOption option;

135

double median = 0;

136

double sd = 0;

137

int read_length = 0;

138

int max_isoforms = 3;

139

int max_component_size = 30;

140

141

void BuildHashGraph(int kmer_size);

142

void Assemble(HashGraph &hash_graph);

143

void AlignReads(const string &contig_file, const string &align_file);

144

void CorrectReads(int kmer_size);

145

void LocalAssembly(int kmer_size, int new_kmer_size);

146

void Iterate(int kmer_size, int new_kmer_size);

147

void FindIsoforms(ContigGraph &contig_graph, deque<Sequence> &transcripts);

148

149

bool CompareCoverage(const ContigGraphVertexAdaptor &x, const ContigGraphVertexAdaptor &y)

150

{

151

return x.coverage() > y.coverage();

152

}

153

154

bool CompareCoveragePath(const ContigGraphPath &x, const ContigGraphPath &y)

155

{

156

return CompareCoverage(x.back(), y.back());

157

}

158

159

double Similarity(ContigGraphPath &x, ContigGraphPath &y)

160

{

161

deque<ContigGraphVertexAdaptor> a;

162

deque<ContigGraphVertexAdaptor> b;

163

164

for (unsigned i = 0; i < x.num_nodes(); ++i)

165

a.push_back(x[i]);

166

for (unsigned i = 0; i < y.num_nodes(); ++i)

167

b.push_back(y[i]);

168

169

sort(a.begin(), a.end());

170

sort(b.begin(), b.end());

171

172

unsigned i = 0;

173

unsigned j = 0;

174

double common = 0;

175

while (i < a.size() && j < b.size())

176

{

177

if (a[i] < b[j])

178

++i;

179

else if (a[i] > b[j])

180

++j;

181

else

182

{

183

common += a[i].contig_size();

184

++i;

185

++j;

186

}

187

}

188

189

return common / max(x.size(), y.size());

190

}

191

192

void FindIsoforms(ContigGraphPath &path, deque<ContigGraphPath> &isoforms, ContigGraph &contig_graph, deque<ContigGraphPath> &component_isoforms)

193

{

194

if ((int)isoforms.size() >= max_isoforms)

195

return;

196

197

int kmer_size = contig_graph.kmer_size();

198

199

ContigGraphVertexAdaptor current = path.back();

200

if (current.out_edges().size() != 0)

201

{

202

deque<ContigGraphVertexAdaptor> candidates;

203

for (int x = 0; x < 4; ++x)

204

{

205

if (current.out_edges()[x])

206

{

207

ContigGraphVertexAdaptor next = contig_graph.GetNeighbor(current, x);

208

209

if (next.status().IsUsed())

210

continue;

211

212

candidates.push_back(next);

213

}

214

}

215

216

sort(candidates.begin(), candidates.end(), CompareCoverage);

217

218

for (unsigned i = 0; i < candidates.size(); ++i)

219

{

220

ContigGraphVertexAdaptor next = candidates[i];

221

path.Append(next, - kmer_size + 1);

222

next.status().SetUsedFlag();

223

FindIsoforms(path, isoforms, contig_graph, component_isoforms);

224

next.status().ResetUsedFlag();

225

path.Pop();

226

}

227

}

228

else

229

{

230

bool flag = true;

231

232

if (flag && path.size() > 300)

233

{

234

isoforms.push_back(path);

235

}

236

}

237

}

238

239

bool Compare(int x, int y)

240

{

241

deque<ShortSequence> &reads = assembly_info.reads;

242

if (reads[x] != reads[y])

243

return reads[x] < reads[y];

244

else

245

return reads[x+1] < reads[y+1];

246

}

247

248

void SortReads(AssemblyInfo &assembly_info)

249

{

250

deque<ShortSequence> &reads = assembly_info.reads;

251

deque<int> aux;

252

for (unsigned i = 0; i < reads.size(); i += 2)

253

{

254

if (reads[i+1] < reads[i])

255

swap(reads[i], reads[i+1]);

256

aux.push_back(i);

257

}

258

259

sort(aux.begin(), aux.end(), Compare);

260

261

deque<ShortSequence> new_reads;

262

int last = -1;

263

for (int i = 0; i < (int64_t)aux.size(); ++i)

264

{

265

int id = aux[i];

266

if (last == -1 || reads[id] != reads[last] || reads[id+1] != reads[last+1])

267

{

268

new_reads.push_back(reads[id]);

269

new_reads.push_back(reads[id+1]);

270

}

271

272

last = id;

273

}

274

275

reads.swap(new_reads);

276

}

277

278

int main(int argc, char *argv[])

279

{

280

OptionsDescription desc;

281

282

desc.AddOption("out", "o", option.directory, "output directory");

283

desc.AddOption("read", "r", option.read_file, FormatString("fasta read file (<=%d)", ShortSequence::max_size()));

284

desc.AddOption("long_read", "l", option.long_read_file, FormatString("fasta long read file (>%d)", ShortSequence::max_size()));

285

desc.AddOption("mink", "", option.mink, FormatString("minimum k value (<=%d)", Kmer::max_size()));

286

desc.AddOption("maxk", "", option.maxk, FormatString("maximum k value (<=%d)", Kmer::max_size()));

287

desc.AddOption("step", "", option.step, "increment of k-mer of each iteration");

288

desc.AddOption("inner_mink", "", option.inner_mink, "inner minimum k value");

289

desc.AddOption("inner_step", "", option.inner_step, "inner increment of k-mer");

290

desc.AddOption("prefix", "", option.prefix_length, "prefix length used to build sub k-mer table");

291

desc.AddOption("min_count", "", option.min_count, "minimum multiplicity for filtering k-mer when building the graph");

292

desc.AddOption("min_support", "", option.min_support, "minimum supoort in each iteration");

293

desc.AddOption("num_threads", "", option.num_threads, "number of threads");

294

desc.AddOption("seed_kmer", "", option.seed_kmer_size, "seed kmer size for alignment");

295

desc.AddOption("min_contig", "", option.min_contig, "minimum size of contig");

296

desc.AddOption("similar", "", option.similar, "similarity for alignment");

297

desc.AddOption("max_mismatch", "", option.max_mismatch, "max mismatch of error correction");

298

//desc.AddOption("min_pairs", "", option.min_pairs, "minimum number of pairs");

299

//desc.AddOption("max_gap", "", option.max_gap, "maximum gap in reference");

300

desc.AddOption("no_local", "", option.is_no_local, "do not use local assembly");

301

desc.AddOption("no_coverage", "", option.is_no_coverage, "do not iterate on coverage");

302

desc.AddOption("no_correct", "", option.is_no_correct, "do not do correction");

303

desc.AddOption("pre_correction", "", option.is_pre_correction, "perform pre-correction before assembly");

304

desc.AddOption("max_isoforms", "", max_isoforms, "maximum number of isoforms");

305

desc.AddOption("max_component_size", "", max_component_size, "maximum size of components");

306

307

try

308

{

309

desc.Parse(argc, argv);

310

311

if (option.read_file == "" && option.long_read_file == "")

312

throw logic_error("not enough parameters");

313

314

if (option.maxk < option.mink)

315

throw invalid_argument("mink is larger than maxk");

316

317

if (option.maxk > (int)Kmer::max_size())

318

throw invalid_argument("maxk is too large");

319

}

320

catch (exception &e)

321

{

322

cerr << e.what() << endl;

323

cerr << "IDBA-Tran - Iterative de Bruijn Graph Assembler for next-generation transcriptome sequencing data." << endl;

324

cerr << "Usage: idba_tran -r read.fa -o output_dir" << endl;

325

cerr << "Allowed Options: " << endl;

326

cerr << desc << endl;

327

exit(1);

328

}

329

330

MakeDir(option.directory);

331

332

LogThread log_thread(option.log_file());

333

334

string begin_file = option.directory + "/begin";

335

fclose(OpenFile(begin_file, "wb"));

336

337

if (option.num_threads == 0)

338

option.num_threads = omp_get_max_threads();

339

else

340

omp_set_num_threads(option.num_threads);

341

cout << "number of threads " << option.num_threads << endl;

342

343

ReadInput(option.read_file, option.long_read_file, assembly_info);

344

cout << "reads " << assembly_info.reads.size() << endl;

345

cout << "long reads " << assembly_info.long_reads.size() << endl;

346

347

SortReads(assembly_info);

348

cout << "reads " << assembly_info.reads.size() << endl;

349

cout << "long reads " << assembly_info.long_reads.size() << endl;

350

351

read_length = assembly_info.read_length();

352

cout << "read_length " << read_length << endl;

353

354

if (option.is_pre_correction)

355

{

356

int kmer_size = (option.maxk + option.mink)/2;

357

cout << "kmer " << kmer_size << endl;

358

BuildHashGraph(kmer_size);

359

AlignReads(option.contig_file(kmer_size), option.align_file(kmer_size));

360

CorrectReads(kmer_size);

361

assembly_info.ClearStatus();

362

}

363

364

int old_kmer_size = 0;

365

int kmer_size = option.mink;

366

while (true)

367

{

368

cout << "kmer " << kmer_size << endl;

369

370

if (kmer_size == option.mink)

371

BuildHashGraph(kmer_size);

372

else

373

Iterate(old_kmer_size, kmer_size);

374

375

assembly_info.ClearStatus();

376

AlignReads(option.contig_file(kmer_size), option.align_file(kmer_size));

377

378

CorrectReads(kmer_size);

379

assembly_info.ClearStatus();

380

381

old_kmer_size = kmer_size;

382

kmer_size = min(option.maxk, kmer_size + option.step);

383

LocalAssembly(old_kmer_size, kmer_size);

384

385

if (old_kmer_size == option.maxk)

386

break;

387

}

388

389

kmer_size = option.maxk;

390

391

deque<Sequence> contigs;

392

deque<string> names;

393

ReadSequence(option.contig_file(kmer_size), contigs, names);

394

FastaWriter writer(option.contig_file());

395

for (unsigned i = 0; i < contigs.size(); ++i)

396

{

397

if ((int)contigs[i].size() >= option.min_contig)

398

writer.Write(contigs[i], names[i]);

399

}

400

401

string end_file = option.directory + "/end";

402

fclose(OpenFile(end_file, "wb"));

403

404

fflush(stdout);

405

406

return 0;

407

}

408

409

void BuildHashGraph(int kmer_size)

410

{

411

BuildKmerFile(assembly_info, kmer_size, option.min_count, option.prefix_length, option.kmer_file());

412

413

HashGraph hash_graph(kmer_size);

414

ReadKmerFile(option.kmer_file(), hash_graph);

415

416

hash_graph.RefreshEdges();

417

418

if (!option.is_no_internal)

419

InsertInternalKmers(assembly_info, hash_graph, option.min_count);

420

421

if (option.reference != "")

422

{

423

deque<Sequence> ref_contigs;

424

ReadSequence(option.ref_contig_file(), ref_contigs);

425

#pragma omp parallel for

426

for (int64_t i = 0; i < (int64_t)ref_contigs.size(); ++i)

427

hash_graph.InsertUncountKmers(ref_contigs[i]);

428

hash_graph.RefreshEdges();

429

}

430

431

Assemble(hash_graph);

432

}

433

434

void Assemble(HashGraph &hash_graph)

435

{

436

cout << "kmers " << hash_graph.num_vertices() << " "<< hash_graph.num_edges() << endl;

437

438

int kmer_size = hash_graph.kmer_size();

439

double min_cover = max(1, (kmer_size == option.mink ? option.min_count : option.min_support));

440

441

Histgram<int> hist = hash_graph.coverage_histgram();

442

//double expected_coverage = hist.mean();

443

444

deque<Sequence> contigs;

445

deque<ContigInfo> contig_infos;

446

hash_graph.Assemble(contigs, contig_infos);

447

hash_graph.clear();

448

449

{

450

HashGraph tmp_hash_graph;

451

tmp_hash_graph.swap(hash_graph);

452

}

453

454

ContigGraph contig_graph(kmer_size, contigs, contig_infos);

455

contigs.clear();

456

contig_infos.clear();

457

458

if (!option.is_no_coverage)

459

{

460

contig_graph.RemoveStandAlone(kmer_size);

461

462

int bubble = contig_graph.RemoveBubble();

463

cout << "merge bubble " << bubble << endl;

464

465

contig_graph.RemoveLocalLowCoverage(min_cover, option.min_contig, 0.1);

466

}

467

468

contig_graph.SortVertices();

469

contig_graph.GetContigs(contigs, contig_infos);

470

WriteSequence(option.graph_file(kmer_size), contigs);

471

contigs.clear();

472

contig_infos.clear();

473

474

if (!option.is_no_coverage)

475

{

476

double ratio = 0.25;

477

478

deque<Sequence> multi_contigs;

479

deque<ContigInfo> multi_contig_infos;

480

contig_graph.GetContigs(multi_contigs, multi_contig_infos);

481

PrintN50(multi_contigs);

482

483

contig_graph.Trim(10);

484

contig_graph.MergeSimilarPath();

485

contig_graph.GetContigs(multi_contigs, multi_contig_infos);

486

487

contig_graph.InitializeTable();

488

contig_graph.IterateComponentCoverage2(option.min_contig, ratio, min_cover, 1e100, 1.1, max_component_size);

489

contig_graph.GetContigs(multi_contigs, multi_contig_infos);

490

491

contig_graph.Trim(10);

492

contig_graph.Prune(kmer_size);

493

contig_graph.GetContigs(multi_contigs, multi_contig_infos);

494

495

contig_graph.MergeSimilarPath();

496

}

497

498

deque<Sequence> multi_contigs;

499

deque<ContigInfo> multi_contig_infos;

500

contig_graph.SortVertices();

501

contig_graph.GetContigs(multi_contigs, multi_contig_infos);

502

PrintN50(multi_contigs);

503

WriteSequence(option.contig_file(kmer_size), multi_contigs);

504

WriteContigInfo(option.contig_info_file(kmer_size), multi_contig_infos);

505

506

deque<Sequence> transcripts;

507

508

FindIsoforms(contig_graph, transcripts);

509

510

int index = 0;

511

for (unsigned i = 0; i < transcripts.size(); ++i)

512

{

513

if (transcripts[i].size() >= 300)

514

transcripts[index++] = transcripts[i];

515

}

516

transcripts.resize(index);

517

518

PrintN50(transcripts);

519

WriteSequence(option.transcript_file(kmer_size), transcripts, FormatString("transcript-%d", kmer_size));

520

}

521

522

void AlignReads(const string &contig_file, const string &align_file)

523

{

524

deque<Sequence> contigs;

525

ReadSequence(contig_file, contigs);

526

527

HashAligner hash_aligner(option.seed_kmer_size, option.min_contig, 2);

528

hash_aligner.Initialize(contigs);

529

530

int64_t num_aligned_reads = AlignReads(assembly_info, hash_aligner, option.similar, align_file, true);

531

cout << "aligned " << num_aligned_reads << " reads" << endl;

532

}

533

534

void CorrectReads(int kmer_size)

535

{

536

if (option.is_no_correct)

537

return;

538

539

deque<Sequence> contigs;

540

deque<string> names;

541

deque<ContigInfo> contig_infos;

542

ReadSequence(option.contig_file(kmer_size), contigs, names);

543

CorrectReads(assembly_info, contigs, contig_infos, option.align_file(kmer_size), option.max_mismatch);

544

//WriteSequence(option.contig_file(kmer_size), contigs);

545

WriteContig(option.contig_file(kmer_size), contigs, contig_infos, FormatString("contig-%d", kmer_size));

546

}

547

548

void LocalAssembly(int kmer_size, int new_kmer_size)

549

{

550

EstimateDistance(option.align_file(kmer_size), median, sd);

551

if (median < 0 || median != median || sd != sd || sd > 2*median)

552

{

553

cout << "invalid insert distance" << endl;

554

deque<Sequence> local_contigs;

555

WriteSequence(option.local_contig_file(kmer_size), local_contigs, FormatString("local_contig_%d", kmer_size));

556

return;

557

}

558

559

deque<ShortSequence> &reads = assembly_info.reads;

560

561

deque<Sequence> contigs;

562

ReadSequence(option.contig_file(kmer_size), contigs);

563

564

LocalAssembler local_assembler;

565

local_assembler.Initialize(assembly_info, contigs);

566

local_assembler.set_num_threads(option.num_threads);

567

local_assembler.set_mink(option.inner_mink);

568

local_assembler.set_maxk(new_kmer_size);

569

local_assembler.set_step(option.inner_step);

570

local_assembler.set_min_contig(option.min_contig);

571

local_assembler.set_insert_distance(median, sd);

572

573

FILE *falign = OpenFile(option.align_file(kmer_size), "rb");

574

int buffer_size = (1 << 20) * option.num_threads;

575

for (int64_t offset = 0; offset < (int64_t)reads.size(); offset += buffer_size)

576

{

577

int64_t size = min((int64_t)buffer_size, (int64_t)(reads.size() - offset));

578

vector<HashAlignerRecord> all_records(size);

579

580

ReadHashAlignerRecordBlock(falign, all_records);

581

#pragma omp parallel for

582

for (int i = 0; i < size; ++i)

583

{

584

HashAlignerRecord &record = all_records[i];

585

586

if (record.match_length != 0)

587

local_assembler.AddReadByHashAlignerRecord(record, offset + i);

588

}

589

}

590

fclose(falign);

591

592

deque<Sequence> local_contigs;

593

594

if (!option.is_no_local)

595

local_assembler.Assemble(local_contigs);

596

597

int num_seed_contigs = 0;

598

for (unsigned i = 0; i < contigs.size(); ++i)

599

{

600

if ((int)contigs[i].size() > option.min_contig)

601

++num_seed_contigs;

602

}

603

604

cout << "seed contigs " << num_seed_contigs << " local contigs " << local_contigs.size() << endl;

605

WriteSequence(option.local_contig_file(kmer_size), local_contigs, FormatString("local_contig_%d", kmer_size));

606

}

607

608

void Iterate(int kmer_size, int new_kmer_size)

609

{

610

deque<Sequence> contigs;

611

ReadSequence(option.contig_file(kmer_size), contigs);

612

613

deque<Sequence> local_contigs;

614

ReadSequence(option.local_contig_file(kmer_size), local_contigs);

615

616

deque<Sequence> multi_contigs;

617

ReadSequence(option.graph_file(kmer_size), multi_contigs);

618

619

deque<Sequence> transcripts;

620

ReadSequence(option.transcript_file(kmer_size), transcripts);

621

622

uint64_t sum = 0;

623

for (unsigned i = 0; i < contigs.size(); ++i)

624

sum += contigs[i].size();

625

HashGraph hash_graph(kmer_size);

626

hash_graph.reserve(sum);

627

628

deque<Sequence> old_contigs;

629

old_contigs.insert(old_contigs.end(), contigs.begin(), contigs.end());

630

old_contigs.insert(old_contigs.end(), local_contigs.begin(), local_contigs.end());

631

old_contigs.insert(old_contigs.end(), multi_contigs.begin(), multi_contigs.end());

632

old_contigs.insert(old_contigs.end(), transcripts.begin(), transcripts.end());

633

contigs.clear();

634

local_contigs.clear();

635

multi_contigs.clear();

636

637

IterateHashGraph(assembly_info, new_kmer_size, option.min_support, hash_graph, old_contigs);

638

kmer_size = new_kmer_size;

639

old_contigs.clear();

640

641

if (kmer_size < read_length)

642

hash_graph.RefreshEdges();

643

else

644

hash_graph.AddAllEdges();

645

646

Assemble(hash_graph);

647

}

648

649

void FindIsoforms(ContigGraph &contig_graph, deque<Sequence> &transcripts)

650

{

651

transcripts.clear();

652

653

deque<deque<ContigGraphVertexAdaptor> > components;

654

deque<string> component_strings;

655

contig_graph.GetComponents(components, component_strings);

656

657

int num[100] = {0};

658

for (unsigned i = 0; i < components.size(); ++i)

659

{

660

if ((int)components[i].size() < max_component_size)

661

++num[components[i].size()];

662

else

663

{

664

int total = 0;

665

for (unsigned j = 0; j < components[i].size(); ++j)

666

total += components[i][j].contig_size();

667

}

668

}

669

670

int kmer_size = contig_graph.kmer_size();

671

int output_component = 0;

672

673

//FastaWriter transcript_writer(option.label_transcript_file(kmer_size));

674

FastaWriter component_writer(option.component_file(kmer_size));

675

for (unsigned i = 0; i < components.size(); ++i)

676

{

677

for (unsigned j = 0; j < components[i].size(); ++j)

678

{

679

Sequence seq = components[i][j].contig();

680

component_writer.Write(seq, FormatString("component_%d_%d", i, j));

681

}

682

}

683

684

//int tran_index = 0;

685

for (unsigned i = 0; i < components.size(); ++i)

686

{

687

int num_begin = 0;

688

int num_end = 0;

689

for (unsigned j = 0; j < components[i].size(); ++j)

690

{

691

if (components[i][j].in_edges().size() == 0)

692

++num_begin;

693

if (components[i][j].out_edges().size() == 0)

694

++num_end;

695

}

696

697

if (num_end < num_begin)

698

{

699

for (unsigned j = 0; j < components[i].size(); ++j)

700

components[i][j].ReverseComplement();

701

swap(num_begin, num_end);

702

}

703

704

if ((int)components[i].size() <= max_component_size)// * 100)

705

{

706

{

707

++output_component;

708

}

709

710

deque<ContigGraphPath> all_isoforms;

711

for (unsigned j = 0; j < components[i].size(); ++j)

712

{

713

if (components[i][j].in_edges().size() == 0)

714

{

715

ContigGraphPath path;

716

path.Append(components[i][j], 0);

717

components[i][j].status().SetUsedFlag();

718

deque<ContigGraphPath> isoforms;

719

FindIsoforms(path, isoforms, contig_graph, all_isoforms);

720

components[i][j].status().ResetUsedFlag();

721

722

all_isoforms.insert(all_isoforms.end(), isoforms.begin(), isoforms.end());

723

724

for (unsigned k = 0; k < isoforms.size(); ++k)

725

{

726

Sequence contig;

727

ContigInfo contig_info;

728

isoforms[k].Assemble(contig, contig_info);

729

730

transcripts.push_back(contig);

731

732

//transcript_writer.Write(contig, FormatString("transcript_%d_%d", i, tran_index++));

733

}

734

}

735

}

736

737

for (unsigned j = 0; j < all_isoforms.size(); ++j)

738

{

739

for (unsigned k = 0; k < all_isoforms[j].num_nodes(); ++k)

740

all_isoforms[j][k].status().SetUsedFlag();

741

}

742

743

for (unsigned j = 0; j < components[i].size(); ++j)

744

{

745

if (!components[i][j].status().IsUsed())

746

{

747

transcripts.push_back(components[i][j].contig());

748

749

//transcript_writer.Write(components[i][j].contig(), FormatString("transcript_%d_%d", i, tran_index++));

750

}

751

}

752

}

753

else

754

{

755

for (unsigned j = 0; j < components[i].size(); ++j)

756

{

757

transcripts.push_back(components[i][j].contig());

758

759

//transcript_writer.Write(components[i][j].contig(), FormatString("transcript_%d_%d", i, tran_index++));

760

}

761

}

762

}

763

764

contig_graph.ClearStatus();

765

}

766

Older »